1...? .. ..mfiv . ' 3“ .” J6... ”H“ . {Mum 2mm}. 1.. lw .cs ‘ l1. V . it)... i. 53:: , J. . m ....a a. .1)... .. .. ...}...Jflfi. 3%..“ .. $759.1 ‘0 521;. .2 ...». r .I .21 . ..l 9’ .E £92153. 31...} . 2’4. ..QJ)J.H . ... .‘xy; s l 1:5”; STATEU Z IIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIIIIIIIIIIIII 31293 01688 5406 This is to certify that the dissertation entitled A CONTINUUM-BASED SHELL ELEMENT FOR LAMINATED COMPOSITES UNDER LARGE DEFORMATION presented by Chienhom Lee has been accepted towards fulfillment of the requirements for Ph.D. . Mechanics degreein Mas? Major professor/ Date S~/2?/?8 [WSU is an Affirmative Action/ Equal Opportunity Institution 0-12771 I LIBRARY Michigan State Unlverslty PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINE return on or before date due. MTE DUE DATE DUE DATE DUE “M J 82 1/98 6101mm.“ A CONTINUUM-BASED SHELL ELEMENT FOR LAMINATED COMPOSITES UNDER LARGE DEFORMATION By Chienhom Lee A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Materials Science and Mechanics 1998 ABSTRACT A CONTINUUM-BASED SHELL ELEMENT FOR LAMINATED COMPOSITES UNDER LARGE DEFORMATION By Chienhom Lee Driving by lack of accuracy of existing finite element programs in analyzing laminated composites under large deformation, a continuum-based shell element is proposed in this study. The objective is to develop an accurate but inexpensive (in terms of computer time) shell element that can solve large scale engineering problems. The new shell element is based on the Generalized Zigzag Theory to better describe transverse shear stresses and kinky in- plane displacements through the laminate thickness. It also uses the rate-of-deformation tensor and the 'h'uesdell rate of Cauchy stress, in an updated Lagrangian sense, to describe kinematic and kinetic relations for a structure under large deformation. The accuracy of the proposed shell element is demonstrated by comparing its numerical results with several well-recognized investigations based on theoretical and experimental approaches. To Yufang, Grace and Diane iii ACKNOWLEDGMENTS Dear Heavenly Father, thank You for Your grace, which always amazes me. Without Your help, I would never have been able to finish this dissertation. Thank You for putting all the helping hands around me when I needed them: my most sincere appreciation to my advisor, Dr. Dahsin Liu for his wholehearted support and friendship, also for his persistence in keeping hold of academic excellence with love and kindness; my truthful gratitude to my manager, Mr. Chuan Seng Lee for his patience, encouragements, and constant prayers; my special thanks to my company TRW VSSI for the support of my time and equipment; my deepest appreciation to one of my committee members and colleagues, Dr. Ying-Kuo Lee for his mentoring and his inspiring enthusiasm in finding truths in Continuum Mechanics; my gratitude to my genuine friend, Dr. Chun-Ying Lee and my colleague, Dr. Huai-Yang Chiang for their encouragement and advice; my special thankfulness to my committee members, Drs. Nicholas Altiero, Thomas Fence, and Zhengfang Zhou for their help in reviewing and approving this dissertation; my appreciation to Dr. Joop Nagtegaal and Mr. Lance Hill of HKS Inc. for their help in handling the software; my sincere appreciation to those who pray for me and my family, their names are too many to be listed; finally, gratefulness from my heart to my parents who always reflect Your love on me. My Lord, please allow me to have one more request: so many of those mentioned above do not know You. Please show Your mercy to them so that they can also enjoy their lives with You eternally. I know You love them much more than I do. In the glorious name of Jesus Christ, Amen. iv TABLE OF CONTENTS LIST OF TABLES viii LIST OF FIGURES ix 1 INTRODUCTION 1 1.1 Motivation .................................... 1 1.2 Laminate Theories ................................ 2 1.3 Nonlinear Analysis . ............................... 7 1.4 Formulation for Large Defamation ....................... 8 1.5 Problem— Solving Methodology .......................... 9 1.6 Organization of the Thesis ............................ 10 2 LARGE DEFORMATION ANALYSIS 14 2.1 Introduction .................................... 14 2.2 Descriptions of Kinematic Relations ...................... 15 2.3 Rate of Deformation ............................... 16 2.4 Descriptions of Kinetic Relations ........................ 17 2.5 Measure of Objective Stress Rate ........................ 18 2.6 Rate Form of Static Equilibrium Equations .................. 22 3 GOVERNING EQUATIONS 3.1 Introduction .................................... 3.2 Differential Equations .............................. 3.3 Variational Formulation ............................. 3.4 Variational Equations .............................. 3.5 Linearization of Variational Equations ..................... 3.6 Cauchy Stress Update .............................. 3.7 Newton-Raphson Iteration ............................ 4 INCREMENTAL DISPLACEMENT FIELD 4.1 Introduction .................................... 4.2 Displacement Field ................................ 4.3 Displacement Continuity Conditions ...................... 4.4 Interlaminar Shear Stress Continuity Conditions ............... 4.5 Free Surface Shear Stress Conditions ...................... 5 FINITE ELEMENT FORMULATION 5.1 Introduction .................................... 5.2 Incremental Strains ................................ 5.3 Incremental Displacement Gradients ...................... 5.4 Finite Element Description of Displacements ................. 5.5 Description of Geometry ............................. 5.6 Stiffness Matrices and Force Vectors ...................... 6 NUMERICAL STUDIES 6.1 Introduction .................................... vi 26 26 27 28 29 29 33 34 35 35 36 36 38 4O 44 44 45 47 50 53 56 62 62 6.2 Cross-Ply Laminate Under Cylindrical Bending ................ 63 6.3 Bending of Rectangular Laminate ........................ 74 6.4 Shear Locking ................................... 80 6.5 Patch Tests .................................... 81 6.6 Large Deflection Effects in Unsymmetric Cross-ply Composite Laminates . 84 6.7 Comparison of Laminates Under Large Deformation with Experimental Data 93 7 CONCLUSIONS AND RECOMMENDATIONS 98 7.1 Conclusions .................................... 98 7.2 Recommendations ................................ 100 A CONSISTENT LINEARIZATION 102 B SOME SPECIAL MATRICES FOR FINITE ELEMENT FORMULA- TION 105 C TRUESDELL INCREMENTAL CONSTITUTIVE EQUATION 117 D ABAQUS/Standard CONVERGENCY CRITERIA 119 D.1 The Solution of Nonlinear Problems ...................... 119 D2 Steps, Increments, and Iterations ........................ 120 D.3 Convergency ................................... 121 DA Automatic Incrementation Control ....................... 122 BIBLIOGRAPHY 125 vii 1.1 6.1 6.2 6.3 6.4 LIST OF TABLES Summary of various laminate theories ..................... 12 Normalized, central transverse-deflection of a simply supported, isotropic, square laminate under uniform transverse-pressure .............. 77 Normalized deflection and stresses of a simply supported, orthotropic, square laminate under uniform transverse-pressure .................. 78 Normalized, central transverse-deflection of a simply supported, rectangular [0/ 90/0] laminate under sinusoidal transverse-pressure ............ 79 Results of membrane patch test ......................... 83 viii 1.1 2.1 2.2 4.1 5.1 5.2 6.1 6.2 6.3 6.4 LIST OF FIGURES List of considerations and flow chart of the thesis organization ....... 13 Illustration of relative motion between two particles in space ........ 24 Illustration of time and material derivatives of Cauchy stress under a rigid body rotation ................................... 25 Coordinate systems of laminate shell ...................... 43 Shape functions .................................. 60 Illustration of element and natural coordinate systems for a typical element 61 In-plane displacement distribution of a 8:4, [0/ 90/0] laminate under cylin- drical bending .................................. 66 In-plane normal stress distribution of a 8:4, [0 / 90/ 0] laminate under cylin- drical bending .................................. 67 'Ii'ansverse shear stress distribution of a 8:4, [0/ 90/0] laminate under cylin- drical bending .................................. 68 'II'ansverse shear stress distribution of a S210, [0/ 90/0] laminate under cylin- drical bending .................................. 69 ix 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17 Comparison of U101 and the Generalized Zigzag theoretical solutions for in-plane displacement distribution ....................... 70 Comparison of U101 and the Generalized Zigzag theoretical solutions for in—plane normal stress distribution ....................... 71 Comparison of U101 and the Generalized Zigzag theoretical solutions for transverse shear stress distribution ....................... 72 Comparison of U101 and the Generalized Zigzag theoretical solutions for transverse deflection ............................... 73 Small load vs. transverse deflection curve of a S=225, [0/90] unsymmetric laminate ...................................... 87 Small load vs. in-plane force curve of a S=225, [0/ 90] unsymmetric laminate 88 Large load vs. transverse deflection curve of a S=225, [0/90] unsymmetric laminate ...................................... 89 Large load vs. transverse deflection curve of a S=11.25, [0/90] unsymmetric laminate ...................................... 90 Large load vs. in-plane force curve of a S=11.25, [0/90] unsymmetric laminate 91 Length-to-thickness ratio vs. transverse deflection curve of a small loading, [0/90] unsymmetric laminate .......................... 92 Normalized, central transverse-deflection of a [60].: laminate under uniform pressure ...................................... 95 Normalized, central transverse-deflection of a [0/90]2 laminate under uniform pressure ...................................... 96 Normalized, central transverse-deflection of a {—604 /604] laminate under uni- form pressure ................................... 97 D.1 First and second iteration in an increment ................... 124 xi Chapter 1 INTRODUCTION 1.1 Motivation With their superior strength-toweight ratio, stifl'ness tailored-ability, formability, and envi- ronmental stability, advanced polymer-matrix composite materials have become a key ingre- dient in the design of future automotive and aircraft components. The complex structural responses and design versatility associated with such materials provide a strong motivation for developing a numerical design tool (besides experimental trial-and-errors). Many finite element programs have been widely used in the automotive and aerospace industries to investigate structural problems in the context of large deformation. Yet their lack of accuracy in analyzing laminated composite structures has been well recognized. The reason is simply because most programs are designed for isotropic materials. On the other hand, although many laminate theories have been postulated to improve the accuracy of laminated composite analysis, very few of them have been adopted and introduced to the finite-element community. The primary reason is Ahat they are too sophisticated to implement. Obviously, in the finite element simulation of laminated composites, there exists a gap between what has been developed and advanced and what has been actually used. The motivation behind this thesis, therefore, is to bridge such a gap. In order to achieve this goal, a shell element is to be formulated to fulfill the accuracy requirements for both displacement field and stress state and yet retain computational efficiency. The method of achieving the goal is to first conduct thorough reviews on all existing large deformation theories used in finite element formulation as well as the latest development of laminate theories for composite analysis. Compromises between accuracy and efficiency are then taken in developing a new shell element that has degrees of freedom that are low enough for reasonable efficiency in computation and yet high enough for acceptable accuracy in results. Figure 1.1 outlines what need to be considered in order to formulate a laminated com- posite shell element for large deformation analysis. The following sections provide literature reviews of existing laminate theories and nonlinear finite element analysis. Also included are descriptions of the problem-solving methodology and organization of the thesis. 1.2 Laminate Theories When a laminated composite structure is subjected to out-of-plane loading, accurate trans- verse stresses are very important in structural design and failure analysis . If a laminated composite structure is moderately thick or thick, or if it consists of a matrix with low shear modulus, the transverse shear deformation may not be negligible. In addition, it should be noted that delamination is a primary damage mode in laminated composites. Both central delamination and edge delamination occur in laminated composite structures quite often. In general, central delamination may occur as a result of impact loading [41] and edge delamination may be attributed to free edge effect [56]. Since the interlaminar stresses are responsible for delamination, a correct prediction of transverse stresses is critical to laminated composite analysis. Various laminate theories for composites analysis have been proposed in the past and are briefly summarized as follows. (1) First-order Shear Deformation Theory The First-order Shear Deformation Theory, developed by Reissner [65] and Mindlin [51] independently, is known as the Reissner-Mindlin (RM) Theory. It relaxed the Kirchhoff- Love hypothesis that required normals to the mid-plane to remain normal throughout de- formation. By including two additional rotational degrees of freedom, the normals are then free to rotate with respect to the mid-plane during deformation. This type of deforma- tion implies constant transverse shear stresses through the shell thickness. Consequently, shear correction factors are required for the equilibrium process. The RM Theory provides accurate displacements and stresses for thin and moderately thick, isotropic structures. However, as can be seen later, the theory leads to unsatisfactory displacements and stresses for laminated composites. Nevertheless, the finite element formulation based on RM The ory is still the most widely used in commercial software for investigating structures made of conventional metals and composites. The reason is believed to be its high efficiency in computation resulting from using only five degrees of freedom. (2) High-order Shear Deformation Theories Many refined shear deformation theories have been presented to improve the predic- tion of displacements and stresses (especially transverse stresses) for laminated composites. Literature reviews regarding these theories can be found in the book by Palazotto and Dennis [55] and the dissertation by Li [37]. This category of laminate theories is based on an assumed displacement field that is of high-order polynomial functions of the thickness coordinate. Accordingly, both the in—plane and transverse displacements are smooth and continuous through the laminate thickness. In reality, however, the abrupt change of ma- terial properties across the laminate interfaces usually results in kinky distributions of the in—plane displacements. Another deficiency of the High-order Shear Deformation Theories is the prediction of double-valued transverse stresses on the laminate interfaces. This results from the single-valued strains on the interface being multiplied by difl'erent material prop- erties in different layers. This deficiency originates in the assumption of continuous in-plane displacements, resulting in continuous strain distribution through the laminate thickness. Although this unsatisfactory result can be avoided by a recovering technique based on equi- librium equations, a theory that can give correct displacement field and stress state based on constitutive equations is strongly preferred. (3) Layerwise Theories In order to resolve the deficiencies of the High-order Shear Deformation Theories, it may be intuitive to describe each composite laminate as an assembly of individual layers. Some quasi-threedimensional techniques based on so-called Layerwise Theories were pro- posed [6, 33, 37, 46]. These theories treated each layer individually and imposed one or more continuity conditions on the laminate interfaces to preserve the kinky displacement distributions and continuous shear stress states through the laminate thickness. As a result, the total number of degrees of freedom were reduced. Li and Liu [38] showed that the afore- mentioned High-order Shear Deformation Theories were merely simplified cases of the Gen- eralized Layerwise Theory that they proposed. Despite the high accuracy of displacements and stresses obtained from the Layerwise Theories, a large number of degrees-of-freedom proportional to the total number of layers in the laminate was needed. Thus, the theories are computationally inefficient. This disadvantage is especially true when the layer num- ber of composite laminates becomes overwhelmingly large. In view of the advantages and disadvantages of the High-order Shear Deformation Theories and the Layerwise Theories, a theory that would be a compromise between the numerical accuracy and computational efficiency is highly desired. (4) Zigzag Theories A group of theories called Zigzag Theories, [17, 37, 72, 73], uses a presumed displacement field for each layer and utilizes interlaminar continuity conditions (displacement, stresses, or both) to assemble individual layers. The name ”zigzag” is due to their capability of rep- resenting the kinky distributions of in-plane displacements through the laminate thickness when the laminated composite is subjected to bending. Similar to the Layerwise Theory, the Zigzag Theories do not need any shear correction factor; and consequently all the transverse shear strains and stresses can be calculated based on the constitutive equations. The Generalized Zigzag Theory (ZIGZAG) presented by Li and Liu [38] indicated that all other Zigzag Theories were special cases of theirs. As pointed out by Liu and Li [42], the third-order Generalized Zigzag Theory is a simplified case of the third-order Generalized Layerwise Theory because the layer-dependent variables were only designated to the wrath- order and the first-order terms (as opposed to all four terms in the Generalized Layerwise Theory). Therefore, the total number of degrees of freedom of the Generalized Zigzag Theory was layer-number dependent. It was then reduced to be layer-number independent after using the continuity conditions of displacement and interlaminar shear stress. Hence, the major advantage of the Generalized Zigzag Theory over the Generalized Layerwise Theory is that the number of degrees of freedom is independent of the number of layers. Thus, it gives higher computational efficiency. However, it is also because there is no layer-dependent variables that the Generalized Zigzag Theory is less accurate than the Generalized Layerwise Theory. In comparison with the previously mentioned High-order Shear Deformation Theories, the Generalized Zigzag Theory gives correct kinky in-plane displacement fields and trans- verse stress states. However, due to the fact that the assumed displacement fields for the Generalized Zigzag Theory require transverse deflection to be Cl continuous in finite ele- ment formulation, additional degrees of freedom must be introduced. The same situation also happens to any of the high-order shear deformation theories. The Generalized Zigzag Theory is therefore more complicated in formulation and less efficient than the First-order Shear Deformation Theory. (5) Quasi-layerwise Theories and Others Following the same method described in the Generalized Zigzag Theory, Li [37] devel- oped some Quasi-layerwise Theories with the use of only two layer-dependent variables in different orders. Although the accuracy in displacements and stresses can be improved more or less in comparison with the Generalized Zigzag Theory, the Quasi-layerwise Theo- ries sufl'er numerical deficiency because they are very sensitive to the selection of coordinate systems. Global-local Superposition and Double-superposition Theories are two other ideas presented by Li and Liu [39]. The theories utilize the thickness coordinate of a local layer in combination with the laminate thickness coordinate. As a result, the total number of degrees of freedom of these theories is independent of the total number of layers in the lam- inate. Although some of these theories can provide higher accuracy, they are less eflicient than the Generalized Zigzag Theory because a larger number of degrees of freedom is re- quired. Table 1.1 summarizes the displacement fields for various laminate theories discussed herein. The Generalized Zigzag Theory is chosen to be used in this thesis because it results correct kinky in-plane displacements and continuous transverse stresses. More importantly, the Generalized Zigzag Theory has a constant number of degrees of freedom. 1.3 Nonlinear Analysis All laminate theories discussed in the previous section can be used for both linear and non- linear finite element analyses: A linear analysis is used when deformation is small, material response is linear, and boundary conditions remain unchanged during the course Of deforma- tion. In general, nonlinear analysis would be otherwise considered. There are two categories .....m _ of nonlinearity: material nonlinearity and geometric nonlinearity [8]. Material nonlinearity is associated with nonlinear elastic or plastic deformations. Geometric nonlinearity occurs when a structure is subjected to large strains and/or large rotations. In this thesis, the mainrfocus is on geometrically nonlinear problems under static loading conditions. The essential feature of geometric nonlinearity is that equilibrium equations must be written with respect to an instantaneous state [14]. A large deformation problem can be analyzed using either Lagrangian (material) description or Eulerian (spatial) description. The Lagrangian description is also called total Lagrangian. When this approach is used, movements of material particles are described with respect to the original or undeformed configuration. In other words, regardless how large the strain and rotation are, all displace- ment differentiat ions and integrations are performed with respect to the original frame. As deformation becomes larger and larger, more and more terms (usually nonlinear) must be added to the strain-displacement relations in order to account for the nonlinearity. When the Eulerian description is used, movements of material particles are described with respect to the current or deformed configuration. In actual implementation, the Eu- lerian approach takes a form that is usually called updated Lagrangian. In this approach, differentiations and integrations are performed with respect to the deformed configuration. The current deformed configuration is also used as the reference state prior to the next increment of the solution. After the incremental solution is obtained, the reference state is updated and then the solution proceeds to the next increment. It is noted that although different formulations may exist when using different ap- proaches (one may be more complicated than another), final solutions to a problem should be identical. In the total Lagrangian approach, the kinematic relations are always non- linear because the deformation is usually given by a displacement field. In the updated Lagrangian approach, deformation can be described either by a displacement field or by a velocity field (see Section 2.2 for details) [48]. When a velocity field (such as the rate-of- deformation tensor) is utilized to describe the deformation, the kinematic relations become linear. When dealing with laminated composites, displacement fields are in general very complex. Therefore, it is preferred to use linear kinematic relations. 1.4 Formulation for Large Deformation Lea Although the advancement in the computational techniques for structural analysis is very significant in the last two decades, the development of finite element formulations for lam— inated composites subjected to large deformation is very limited. Many commercial pro- grams such as ABAQUS, LS-DYNA3D, PAM-CRASH, and RADIOSS CRASH, and publi- cations [2, 11, 24, 68, 78] use the Reissner-Mindlin Theory with various updated Lagrangian approaches. Most of them imply that their large deformation finite element formulations are valid not only for isotropic materials but also for laminated composites. However, although its prediction on overall behavior of structures, such as transverse deflections may be ac- ceptable, the Reissner-Mindlin Theory gives incorrect in-plane displacement and transverse shear stresses for laminated composites. As mentioned before, numerous studies have been presented using different lam- inate theories to improve the accuracy of simulating laminated composites in linear analysis. However, when their techniques were extended to large deformation anal- ysis of symmetric or unsymmetric laminated plates and beams subjected to bending, most of them used a total Lagrangian approach and the von Kdrmén nonlinear strains [7, 12, 26, 34, 35, 63, 66, 69, 73, 74]. The van Kérmén nonlinear strains are a simplifica- tion of the Green (Lagrangian) strain tensors with some nonlinear terms eliminated. Most researchers used them for small rotation and small strain nonlinear problems, although the theory can be used for moderately large rotations. Liao and Reddy [40] used a three- dimensional degenerated shell element along with the Green strain tensor and the second Piola-Kirchhofl' stress tensor, which is a total Lagrangian approach, to study post-buckling behaviors of stiffened composite shells. Kweon, et a1. [30, 31] also used the Green strain tensor and the second Piola-Kirchhofl' stress tensor along with the First-order Shear Defor- mation Theory to study the postbuckling compressive strength of graphite/epoxy laminated cylindrical panels. In their book, Palazotto and Dennis [55] used a forth-order shear defor- mation theory, in addition to the Green strain tensor and the second Piola-Kirchhofl' stress tensor, to formulate a shell element. 1.5 Problem-Solving Methodology As seen in the previous section, there has been lack of a shell element with computational efficiency that would accurately describe behaviors of laminated composites under large deformation (large rotation and large strain). Therefore, a continuum-based shell element basedpn the Generalized Zigzag Theory is proposed in this thesis. To avoid the complexity of involving nonlinear strain-displacement relations, the formulation for large deformation adopts an updated Lagrangian approach based on the rate-of-deformation tensor and the 'I‘ruesdell“ rate of Cauchy stress. The rate-of-deformation tensor is a linear formulation. 10 When it is used with the Truesdell rate of Cauchy stress, a symmetric stiffness matrix can be achieved. More importantly, the rate-of-deformation tensor is not a simplification of any strain-displacement relations. Hence, it can be used for investigations involving large strain and/or large rotation. The finite element formulation in this thesis leads to a four-node shell element, which has seven degrees of freedom at each node. , i If" The finite element formulation was programmed as a subroutine called LACOS (LAmi- nated COmposite Shells) using FORTRAN and was then linked to ABAQUS/ Standard as a new addition to its element library. The user defined element is called U101 in the element library. It follows the naming convention required by ABAQUS/ Standard. The elements of each finite element model are then assembled in global coordinates and solved iteratively by the ABAQUS / Standard solver. ABAQUS / Standard is a general nonlinear finite element package that has been developed and maintained by Hibbitt, Karlsson and Sorensen(HKS), Inc. for more than two decades. The concept of employing user subroutines is to make the most use of the commercial package’s existing ftmctions, i.e. solver, post-processing, etc., while the user is still able to tailor the program for specific applications. The subroutine developed in this thesis can also be used in other similar finite element packages with slight modifications. 1.6 Organization of the Thesis A flow chart of the thesis organization is illustrated in Figure 1.1. We begin in Chapter 2 by establishing some basic knowledge about objective stress rates and the rate-of-deformation tensor before we start deriving governing equations. The rate form of static equilibrium equations is also discussed to show its essential differences from equilibrium equations. In Chapter 3, derivations of the governing equations for a structure undergoing large defor- 11 mation (both large strain and large rotation) are described. In Chapter 4, incremental displacement fields using the Generalized Zigzag Theory are presented. In Chapter 5, de- scriptions of finite element formulation for the U101 element are given. In Chapter 6, various numerical studies are presented to evaluate the performance of U101. Finally, conclusions and recommendations are given in Chapter 7. Table 1.1 Summary of various laminate theories Theory Displacement Fields D.O.F. Remarks 1. Classic Plate Theory u(1:,y, 2) = 110(35, y)- 211).: v(z,y,z) = v()($,y) — zwyy 3 Global what/,2) = “Mid/l 2. Reissner-Mindlin 11(21. y. 2) = U003, y) + 3111x0541) Theory U(a:,y, z)“ — ’U0(.’17 ,y) + ztby(z,y) 5 Global w($9y7 Z).- _ 1110(37 y) 3 E 3. Generalized High-order u($,y,2)— — Uo($ my, y)+ Z¢I($7y) + ZQCAIEJJ) + Z 431%!» 3“?“ Dem’ma‘m” The” um, z) = way) + z¢. + 22w, y) + may) 7 Global (third order) U)(CU, yaz Z): w(l($ay) 4. Generalized Higher- order Shear Deformation ”(3 ’zy’ :zzoz “1(33, y) Theory ( 4th order and 2 1 1 4 GI b l . - o a higher) (2) U<$ 31,2 _:z 711(33 y) I (”1+ )+ I i: 0 10(3) yiz)=w0$( ay> ' 5. Generalized LayenNise u"($, y,z 2:) 5(z ,y) + 2qu + zzu§(x, y) + z“u§(:r,, y) Theory (thud order) (1) vk(:c, y,z z): affix, y) + zv’f + 22UI§(3:, y) + z3v§(m, y) 4(k+1)+1 Local “4131/, Z): MU< 9y) 4 6. Generalized Zigzag uk($, y,z z): uI§(x ,y) + 2qu + z2u2(z, y) + 2311;5(1‘, y) Thea third order 1 . Global- “ I I I m, 11.2 z): (is. y>+ at + 22v2($i y) + am. y) 7 Low. w<$vya Z): “10(1),?” 7. Generalized Higher- M order Zigzag Theory (4th ”km" y’ 7’) : ”GIL 3’) + ZUI($* y) + 2211411.!) i=2 order and hi her) (1),(2) - g k M 1. [2(M+1)-4]+1 G'Obal' v (:r,y, z—) — 115(3, y) + zv1(:z:, y) + Z 2 v,(:r,y) L003 i=2 wlf-Tayi 2): 1110(1)?!» 8. Quasi-layewvise Theory (third order) (1),(3) kII ’Zy’ F2514; ] 3 . 7 Global- [ vk(x,y,z)— - funny) Local I ' =0 1 w(:r,y, 2'): “10(1: y: 9. Global-local n 3 ~ (third order) (1),(4) :0 Globa,_ n i 11 Way. 2): (£10m vm($ y)+ +(€k) ”UM-7&9) +22 ill-(Iris!) Loca' i=0 1003,9374) = wll(mt y) 10. Double Superposition m n 3 . Theory (third order) (1)|(5) uktra y, z) 2 (£16) 115422.11) + (6k) 15543373!) + ($01)“???in 'I' :0 2111,,(1: n. 3 . Global- "U"($.y. 2)= (Eklm vm(at y)+ (5.) vibe, y) + (éklpvflay) + ZZ‘vl-(ziy) 13 Local Mani/,2) = wo(rr,y) (1) k=number of layers (2)M = number of terms in the polynominal (3) k=0 in any two out of four terms (4) m.n=1,1,2, or 3; math (5)m,n,p=1,2, or 3; mensp l3 comaanmqams £35 mo 98:0 Boa 93 28330238 mo .33 fin 23mg cosmiccou EoEom. BEE \ . _ z / new 207.3050 r Lr iiiiiiiiiiiii lllllllllllllll‘ Eco—.... usuEN , 20.83am finenwtm> 22635 5:22 / accumucoEEoomm _ _ _ \ a " uo~=EocoO 9:. 9.6: 8.2“. . EoEoomamE 558905 accumtm> 35395 .o 98:33 3:39.000 ’______._____———-___-— F .2920 I 53.2.6 mm_Eo>m_-_md:O $85980 $55234 «529w 065280 morons 8:223. 35:0 3 lllllllll \ / v 5520 228 new 32.9528 597. 59958229: 20:25.23 82m 5296:”. \—————_____—_—_——-——_—_/ 8.2”. 238230 iiiiiiiiiiiiiiiiiiiii I Illlliiliiilicul Tull ‘ _ 39% >5an " 3 2a.. :5an _ I I u 323 3250 30:8 " .0 82 :03th cogs—5.21903: _ I 3 I _ 8:95. 89.2 " ESE _mcmE_< _ u 8ch u 59$ :35 n N 5&ch 83:55. / \ I lllllllllllllllllllllllllllllllllll \ Chapter 2 LARGE DEFORMATION ANALYSIS 2.1 Introduction Once the Generalized Zigzag Theory is chosen as the displacement field, there are many different ways of establishing kinematic and kinetic relations for large deformation analy- sis. Careful considerations must be given and choices must be made before we start the development of governing equations and the subsequent formulation of a shell element that deals with laminated composites. In this chapter, we introduce the principle of objectivity (or material frame indifference) for stress rate and the rate-of-deformation tensor. This will enable us to describe the motion of a particle in a continuum at any instant of time during the process of a large deformation. Numerous studies have been done on this topic, for example, Fung [19] and Malvern [48]. Here, we only outline some of the important concepts needed for future use. Since the rate form of static equilibrium equations is quite different from the equilibrium equations themselves, one section is devoted to present different ways 14 of obtaining the rate form of static equilibrium equations. 2.2 Descriptions of Kinematic Relations When dealing with geometrically nonlinear problems, numerous measures of strain are avail- able. However, most theoretical works and computer programs utilize the following three kinematic relations: 1. The Green strain tensor (also called the Lagrangian strain tensor) . hi)" I’ll/'4}. (11’ J45) 3“,} .. ... 1' all“ an] auk auk , Ar}? 1‘ \1 6'] - 2 (an + BXi + 0X.’ 6X1) : It fl" Z -. " ”2.1) 2.. The Almansi strain tensor (also called the Eulerian strain tensor) . . _ (I ~ rift} '5; = .1. _% + 22]. _ flg‘i 1- 3 i I (2.2) J 2 sz 6223' 633' 8221' 3. The rate-of-deformation tensor (also called the stretching tensor or the velocity strain “3 ‘. '2 .t- {j tensor) D” = . y _ _.L 2.3 , U 2 (6171' + 6.733) ( ) In Eq.(2.l), X.- are components of Lagrangian coordinates. They are also called material coordinates because they describe the material particles with respect to the original or undeformed configuration. An approach is called total Lagrangian when the Lagrangian coordinate system is used to describe the motion of a particle. In Eqs.(2.2) and (2.3), 2:, are components of Eulerian coordinates. They are also called spatial coordinates because they describe the material particles with respect to the current or deformed configuration. An approach is called updated Lagrangian when the Eulerian coordinate system is used to describe the motion of a particle. The two sets of coordinates are related by $3 = X3 + ug (2.4) 16 where u,- are components of a displacement vector. It is noted that both the Green and Almansi strain tensors have nonlinear, coupling terms. The rate-of-deformation tensor is linear, but it is not derived by removing nonlinear terms from the Green or Almansi strain tensors. It is also noted that the rate-of-deformation tensor is calculated from velocity fields while the Green and Almansi strain tensors are calculated from displacement fields. The rate-of-deformation tensor expresses a change of the displacement vector from the current state to the next state along the loading history. Remark 2.1 The von Karmdn nonlinear strains used in many large deformation finite element formula- tions is a simplified, special case derived by removing some nonlinear terms from the Green strain tensor. It is mostly used for small rotation and small strain nonlinear problems. 2.3 Rate of Deformation In this thesis, we choose the rate-of-deformation tensor in Eq.(2.3) as the kinematic relations for large deformation analysis. When dealing with complex displacement fields for laminated composites such as the Generalized Zigzag Theory, it is preferred to use linear kinematic relations. More importantly, the rate-of-deformation tensor is not a simplification of any strain-displacement relations and, hence, can be used for problems involving large strains and/or large rotations. The following is a brief introduction of the rated-deformation tensor. Since our focus in the current study is on the instantaneous motion of a continuum, we describe the motion of a typical particle P inside the continuum by choosing the spatial description and use velocity of the particle as functions of its instantaneous position ($1,132, 1:3) in space. As shown in Figure 2.1, we consider two infinitesimal neighborhood particles P1 and P2 17 with instantaneous coordinates x and x + dx respectively. The relative motion of the two particles from one time instant t to another time instant t+ At can be completely described by a tensor quantity called velocity gradient (LU) defined as follows. 511i L,- = — J 62:,- E U331 (2.5) Lij in Eq.(2.5) can be additively decomposed into two tensors L3} = 031' + WU (2.6) where 1 T 1 Do = 5(Lij + Ln) = 505,1 + v1.1) 5 Um) (2-7) 1 T 1 Wu = 5(L3'J' - Lij) = ‘2-(vm‘ - 12,-,3) E inJ] (2.8) Dij = Dj,‘ and Wij = -Wj3' (2.9) The symmetric tensor ng is called the rate-of-deformation tensor, and the skew-symmetric tensor W3,- is called the spin tensor. The rate-of-deformation tensor 0,5 is a well defined quantity; it vanishes when the continuum performs a rigid-body motion. Therefore, it is an objective quantity. On the other hand, when there is no rigid body motion, the spin tensor Wij becomes zero. 2.4 Descriptions of Kinetic Relations In describing the kinetic relations, a frame indifference (also called objectivity) condition must be satisfied. When using a stress or strain tensor to describe the response of a material, it must be frame indifferent; otherwise, no constitutive relation measured from physical material tests can be established to accomplish the calculation of material response. In addition, a condition called energy conjugate must also be fulfilled. A stress is ”energy 18 conjugate" to the strain if its scalar product with that strain gives equivalent work (energy) to that in a reference frame [15]. Although the Cauchy stress is the energy conjugate to the time integration of the rate- of-deformation tensor, neither its material nor time derivative is frame indifferent [19]. The Cauchy stress a,- j(;r1, .132, 2:3, t) is a time-dependent stress field in a continuum and is referred to a fixed reference coordinate system. It is a true measure of the stress state inside the deformed continuum. The material derivative of the Cauchy stress (1122212912,, sit at 62:." (2.10) indicates the time rate change of a typical stress component at a particle of the contin- uum. For a stressed continuum performing rigid body rotation, neither the “999933313“ (Bag/(9t) nor the material derivative (dagj/dt) of the Cauchy stress vanishes identically. This can be seen in an example illustrated in Figure 2.2 [19]. A bar is subjected to simple tension and rigid body rotation about the z axis. At one instant, when the bar is parallel to y-axis, a; = 0 and ory 91$ 0. At another instant, when the bar is parallel to z-axis, 0; ¢ 0 and 0,, = 0. Accordingly, a rigid body rotation changes the stress tensor while the stress state is unchanged inside the bar. Thus, neither Bag/3t nor dog-J- /dt can serve as an appropriate stress rate measure to be related simply to the rate-of-deformation Dij- In other words, it is impossible to establish a constitutive relation between Djj and dUij/dt (or 80,-]- /8t). This is why an objective stress rate must be introduced. 2.5 Measure of Objective Stress Rate As described by Prager [57], the objective stress rate must vanish when a stressed continuum performs a rigid body motion. Apparently, this restriction is not severe enough to lead to 19 a unique definition of objective stress rate. A variety of definitions of stress rate have been proposed in the literature on mechanics of continua. Among the many frame-indifferent stress rates of the Cauchy stress, the Jaumann (also called Zaremba—Jaumann-Noll) rate and Truesdell rate of Cauchy stress are commonly used. Although many lecturers have contgbutedfitheirweffofts in discussing the superiority of one stress rate to the others, [16, 18, 20, 23, 50, 57, 58], the results were inconclusive. Theoretically, the result of using any one of them should be identical from a continuum mechanics point of view. However, depending on formulations and applications, one stress rate may be more suitable than the other. In this thesis, the Truesdell rate of Cauchy stress is used. When it is used with the rate-ofigformation tensor in an updated Lagrangian approach, a symmetric stiffness matrix can result in the finite element formulation. In the following, the Truesdell rate and Jaumann rate of Cauchy stress are briefly in- troduced. (1) Truesdell Rate of Cauchy Stress The relationship between the symmetric second Piola-Kirchoff (2nd PK) stress tensor and the Cauchy stress tensor is given by [49] W S3)" = JF‘EIO’HFJI (2.11) where ng is the second Piola-Kirchoff stress defined with respect to the material coordinates (X 1,X2,X3), Irij is the deformation gradient, and J = det (Fl-J). Using the method of consistent linearization shown in Appendix A, we have Elsijl = Ellerilladefil)+J£1Fillladefill + J(F;;l)c[a,.,](1r]7‘) + J(F,;‘)a,.,c[FJ;1] (2.12) 20 Now choose the frame of reference X to be the instantaneous motion x, then J=1 and (17,;1)=5,,- (2.13) and from Eqs.(A.9) and (A.ll) cw] = Auk). and cm; = —Au,,,- (2.14) Hence .. '~ f [I a i ' C(Sij] éf Vafj = A0,] + Unguch - O'uAUjJ — Airman]- (2.15) or 11’ :15. . A031 = Vafij - UijAukyc + (7ngij + Auualj (2.16) where Vafj is the Truesdell rate of Cauchy stress and Am,- is the Cauchy stress rate. (2) Jaumann Rate of Cauchy Stress Let us consider a field of flow with velocity components v.- attached to a rectangular Cartesian frame of reference (1:1, 9:2, :33). For convenience of discussion, we follow Fung’s work [19] and take the origin of the coordinate system at a generic point P in the flow field. Let (3’1, 1’2, 23) be another rectangular Cartesian frame of reference that has the same origin at P and rotates with the continuum at an incremental rotation (spin tensor) W, where the components of ng are _1 912. %) _l(§a m) _l(% 211.2) W23-2(6$3 632 ,W31—2 6:121 61:3 ,W12—2 01:2 811 (2.17) Let x: coincide with x,- at any instant of time t. Then the stress tensor at P is 0:1-(t) = cry-(t). Referring to the rotating axes x: at a later instant of time t+dt, let the stress at the particle P be denoted by o{j(t + dt). Then the Jaumann rate of Cauchy stress is defined as l V = . —- ’.. —— I.. a dltlgh dt[0,](t + dt) 0,1(t)] (2.18) 21 Now, the coordinates 2:; and 1:, are related by 87- . 8 - .L‘: = r,- + ,—,'—.rkdt = (0,). + idt) 1:)c (2.19) - dzrk (hi); 1 1a.. - Decomposing the velocity gradient tensor 803/6415); in the equation above into the rate-of- deformation tensor(D,-J-) and the spin tensor (WU) using Eq.(2.6), we obtain 2:; = (6,). + W,,.dt)x,. (2.20) where 0,1 = 0 for rigid body rotation. The stress tensor at particle P at the instant t + dt, with reference to the fixed coordinates 2:1, is 0,,(1 + dt) = a.,-(t) + 33141 (2.21) Transforming the stress tensor of Eq.(2.21) through the coordinate transformation Eq.(2.20) into the 2;; axes, we obtain ago + dt) = (5,, + W,,dt)(6,, + qudt) [.:,..,(t) + 3351111] d d0}! 2 Accordingly, from Eq.(2.18), we obtain da- V __ U , . . . Remark 2.2 The use of the Jaumann rate is to view the stress-strain relation from the standpoint of an observer in a moving material frame relative to which the local rotation vanishes. In other words, the Jaumann stress rate provides an objective measure of the change in stress viewed from a frame rotating with the material. 22 Remark 2.3 The generalized Hooke’s law can be applied as the stress-strain relation between the Jaumann stress rate and the rate of deformation in the coordinates aligned with the material principal directions, assuming the material is linear elastic. The engineering constants are Young ’3 moduli, Poisson’s ratios, and shear moduli that are measured from simple tests such as uniascial tension or pure shear tests. 2.6 Rate Form of Static Equilibrium Equations 52,. . . [25. -. o )1) At any instant of time, without considering the body force and acceleration of a continuum, equilibrium requires that the total force acting on the body must vanish, that is /t3'dS = / 0.1-1;de = 0 (2.24) S S where S is surface of the body and the traction t,- = aijnj with n, being a unit normal to the surface. The divergence theorem allows Eq.(2.24) to be written as [951W=0 (“9.3 l (2.25) V 631' " This equation must be valid for an arbitrary volume V (every portion of the body is in equilibrium), thus, we obtain the familiar stress equations of equilibrium 803' Z’ 6121' = 0 (2.26) As given in Lee [36] and Osias and Swedlow [52], to maintain the deforming body in equi- librium during a loading history, it is required that the time rate of the net applied force be zero. Thus, we have the material derivative of Eq.(2.24) as d . - Eds 1,-dS) = /S(t,-dS+t,-dS) =0 (2.27) 23 Applying the divergence theorem to Eq.(2.27) again will lead to the stress rate equations of equilibrium. Alternatively, we can take the material derivative of Eq.(2.26) directly (also / Q fi (fr (\\, use Eq.(2.lO)), that is f}, Q”; ’3',” ’ d aaij) 820” 6203']- "’ —‘ = = 0 2.28 dt ( 6171' max,- + axkaxj vk ( ) Eq.(2.28) leads to the following equation provided that 01'ch ,- is continuous because org-JU- becomes zero when total stress equilibrium 0,5,]- = 0, i.e. d 6031' 00;”) _ __ = _.l :0 2.29 3’61 0, l: r.. J ’11 3:'\ 93‘; v Usin E . 2.10) ain, we have ...»ml :2. ...,- -~’—- .- ...u-vv» ‘ . :- ‘ " ”' ..~ g Q( as ..-... M (I, ,2,“ J . .. (ti/0,5) _ 603.5 _ 620,3 v _ 30’ng __ 0 (2 30) 6t ,3- - ax, axkoz, : 613,, 8.7:; _ ° (a. ) . .. {3‘ r} , .C - ‘ ‘3 w *" Finally, we have i 4) {$72} ,) ’ 5U": VF? 1°11 -2‘1‘1'2'2 =0 1.5,(‘35L : o (2.31) 6.7:,- Bxk Buy This equation is in a rate form, and, thus, it governs the stress field irrespective of defor- mation magnitude and material structure. Satisfaction of Eq.(2.31) not only implies total 39).} '- 9 . stress equilibrium but also assures that, given an equilibrated stress field, equilibrium is maintained in the presence of time varying loading [36]. 24 Ah 3" ‘Hk v+dv Figure 2.1 Illustration of relative motion between two particles in space A)’ 25 time = t —— rigid body rotation -—-> time = t+At 0,,Oy' ottoy‘ +— » y ’ x, x iOY'Oy v Or-Ox-o OY=03H-0 0v=ay.¢0 ox=oy.¢0 Figure 2.2 Illustration of time and material derivatives of Cauchy stress under a rigid body rotation Chapter 3 GOVERNING EQUATIONS 3.1 Introduction Starting from equations of motion, a variational approach is performed to convert the gov- erning equations into variational equations. Since the variational equations are nonlinear, linearization based on Tayler’s expansion is required to simplify the nonlinear equations, ‘1 -( resulting in linear approximation of the variational equations. Subsequently, the 'Iisuesdell rate of Cauchy stress and the rate-of-deforiiiation tensor are introduced into the formulation to form the final linearized variational equations. Once the linearization formulation is com- pleted, a finite element method can be used to convert the linearized variational equations into finite element formulation, which will lead to a symmetric stiffness matrix. The stiffness matrix can be used in a Newton-Raphson iteration scheme. Therefore, nonlinear solutions can be obtained systematically by using many small linear increments. This procedure is consistent with an updated Lagrangian scheme [8]. It is noted that the procedure presented herein is general and not restricted to any type of deformation (finite or infinitesimal) or structure (bulky or thin-walled). 26 27 3.2 Differential Equations Consider components of a displacement vector u,(.r1,.c2,.rg, t) that satisfies the equations of motion (90" .. 6—2 + bi = P’Ui 3] (3.1) throughout the interior of the body which has a current volume V over the time interval t E (0, T), and is subjected to the following conditions: displacement (essential) boundary condition “i = 9411.32.33. t) 00 Pg“ traction (natural) boundary condition h,(:rl,;rg,:r3, t) = Ugjflj on F)“, and initial conditions ui(X1,X2.X3,0) = "9(X1,X2,X3) a,(X1,X2,X3,0) = v?(X1.X2.X3) where F95 UI‘h‘. = 5 P95 ml“);i = 0 and S is the surface of the continuum. Some variables appeared in Eqs.(3.1) through (3.5) are defined as follows: X ,- are components of the material coordinates, (3.3) 28 2:,- are components of the spatial coordinates and .r, = X, + u,(;r1,;r2, 33, t), org-(x1, 2:2, 223, t) are components of the Cauchy (true) stress tensor, b,(2:1,x2, 1:3, t) are components of the body-force vector per unit volume, p(;rl,.r2, x3, t) is the mass density, and n, are components of the unit normal vector relative to the boundary I‘m. Since the dot over a variable represents time differentiation, 11,- indicates velocity of a ma- terial point and ii,- acceleration of the same point. Furthermore, the indices i and j denote Vrufl 4" Cartesian coordinates relative to a fixed reference frame and they range from 1 to 3. Re- peated indices imply a summation over the range. 3.3 Variational Formulation The variational form of the equations of motion, Eqs.(3.1), can be written as follows / 5111' (% + bi - pu,)dV = 0 (3.6) V 6221' where 6a,- is an arbitrary weight function and must satisfy the homogeneous form of the essential boundary conditions, i.e., C" i live 4 J ,_J _\ 3,4... 3'2 7': {4] J 7 ‘3' i it; 6n,- = 0 on F9.. (3.7) U Integrating by parts of Eqs.(3.6) leads to {V .,.« 662:.- ' ' 6mde+/6u,hdF f6 —0‘,-jdV =/ du,pu,dV (3.8) v It is noted that Eq.(3.8) is a result of transferring differentiation of u,“ to 6m, thus equalizing the continuity requirement on u.- and 6a,- and weakening the requirement on u,- [5]. Now, instead of obtaining an analytical solution to Eq.(3.1), we try to find a numerically approximated solution by using Eq.(3.8) and its associated boundary conditions. LA 1. V \..>‘:‘," cl: ¢- ‘0 29 3.4 Variational Equations In this stage, it is to find u,(.rl,.r2,.r3, t) that satisfies the variational equation F“‘(u.) — F‘"‘(u.-) = Mm.) (3.9) where, by referring to Eq.(3.8) in the previous section, .‘x’ .‘ .t. ’,. .blf‘” " /' [Fext(u,') = [6a,h,dV+/6u,h,df‘ (3.10) V I‘ - adu- ant(u,-) = fv—aj’udV (3.11) M(i1,-) = / amt-W (3.12) v i {_"I Some remarks regarding the above equations are listed below. 1. Components of the variational displacement vector 611,- should satisfy appropriate continuity conditions and 6n; = 0 on P9.- (3.13) 2. Eq.(3.9) is subjected to the following initial conditions u.(X1,X2. X3,0) = u9(xi, X2, X3) a.(X1,X2,X3.0) = v?(X1. X2.X3). (3.14) 3. The traction boundary conditions have been absorbed in the variational process as the surface contribution (fr) in Eq.(3.10). 3.5 Linearization of Variational Equations The governing equation, Eq.(3.9), of a continuum undergoing large deformation (both large strain and large rotation) is generally nonlinear. The consistent linearization procedure described in Appendix A can be used to find the linear approximation of the equation. 30 Taking the function .7: in Appendix A and recalling Eq.(3.9) yields J:(Iv+l)— FexttL' v+l)_ Fant(I:)+l) _ N[(le+l) = 0 (315) The linearized variational equations can be written as arm] + £1F""‘] + £[M] = Fezt(x:’) - p"‘(:c§’) — M(:'i':’) (3.16) Now, if we restrict the derivation to a static case and ignore the contribution of the body- force, we let arm] '= 0, i. "z o, " £[M] = 0 and Mot?) = 0 V (3.17) Eq.(3.16) is then reduced to 5mm] = We?) - F‘“‘(z:-’) (3.18) Therefore, only £[F‘n‘] in Eq.(3.18) needs to be further discussed. By using Eq.(3.11), Eq.(3.18) becomes 66a.- var,- 8—613, 6__Xkl‘ Vo a—Xk a—_j _3_6u,- __1 -1 -1 [V0 {11ijan 1» J" + Pk] cps-)1” + pk] a ”j£[J]}dVo (3.19) arm] = c —o.-,-dv] =£ In Eq.(3.19), we introduce dV = Jdi/b, where V0 is the original (initial) volume of the continuum. By defining L[a,~,-] “.5! A0,,- ‘ (3.20) 31 and using Eqs.(A.9) and (A.11 ), we obtain [Wow ‘ . A as - 6A V cm“) =/V a;;{ 5;; 6;——3F,,'.la,, 0” +17, [A0,] +Fk, a ,-Auu}dV (3.21) Using 6611i F_1 _ 66713" 6X), km '_ 62m BAum _1 _ BAum 6X1 F” — 32:]- (3.22) we obtain '1‘ ~. C , .. 'nt Vaa—iut 0v 1:.[17'1 ] = —j{AO,j + 0 vJ-Aukk- UikAuj,k}dV (3.23) As explained in Chapter 2, the Cauchy stress rate A0,, is not an objective measure of stress; it cannot be involved in any constitutive equation directly. Therefore, Eq. (3. 23) is ...-a...— not useful unless the Cauchy stress rate is further specified. Later in this section, we will u-” revisit the definition of the Truesdell rate of Cauchy stress defined in Eq.(2.16). Recall Eq.(2.5) L" - —— = 12,-J- (3.24) In a static loading condition, the terms ”time increment” and ”velocity” really indicate, respectively, an increment along the loading path and the corresponding increment of dis- placement (Am). Therefore, we rewrite the velocity gradient, rate-of-deformation tensor and spin tensor in Eqs.(2.5), (2.7) and (2.8), respectively, as follows 4 t L,, ; 952-1- : 11,-, (3.25) 0,,- = ~21-(L ,-,- +L ,,) = é-(Auw- + Au“) 2 A116,, (3.26) W,, -;-(L,-,- — L3,) = -21-(Au,-,, — 2111,.) 2 Au“ (3.27) Let‘s define the Truesdell kinetic relations as V0‘ 11 = CitjklAUUcJ) (3-28) 32 where Cf,“ is determined experimentally or transformed accordingly. Here we transform the constitutive equation Cf,“ from the generalized Hookie’s law Czjkl by using the following formula 1 Cf,“ = Cijkl + 0:),de — ~2-(Uyk5jz + 0351* + 03955“ + ayldik) (3.29) The detailed derivation of the above formula is given in Appendix. C. From Eqs.(2.16) and (3.28), we have ,' rt 5 ,. u y A0,, = :,,,,Au(,,,, — [GAILNC + afiAuJ-J + Auuofj (3.30) \~../ 1"“ «1. «ms-"‘5" ,/ Now let's recall the final form of £[Fm‘] from Eq.(3.23) and denote A0,?)- = A0,,- + 0:11AM; — okauJ-J, (3.31) Then Eq.(3.23) becomes 5"“— l - 6611,- "t —— _ !. £[F‘ ]— V 6:5,- Aa,,dV (3.32) By substituting Eq.(3.30) into Eq.(3.31) A0,} ijklAu(k,,) + Anglo]; + aflAuJ-j — okauj,k = ijHAu(k,1) + aaAuu (3.33) . l [I .1 ‘ -. Therefore, ’ ' , .. \ £[F”“‘] = / (SuidCEJ-HAuuAdV + / 6Uf,j0’;"lAUj,ldV (3.34) v v Later in the finite element formulation, the first term of Eq.(3.34) will lead to the material stiffness matrix K "m“ . For most of engineering materials we are going to discuss later, Cf,“ possesses both minor and major symmetries, i.e., 1 _ c _ 1 _ t Cum - jikl - 01ch - 01m: and Citjkl = Clem (3.35) 33 Therefore, only symmetric part of (Suw- is required, i.e., . 1 - - duh,” = §(O’Ug,j + OUjJ) (3.36) Hence, the computation of K m“ is quite similar to those of small deformation analysis. The second term of Eq.(3.34) will lead to the geometrical stiffness matrix K 980'". Since the components of the of, possess major symmetry, a symmetric K 9”" will be obtained from finite element formulation. From Eqs.(3.10), (3.11), (3.18) and (3.34), the final linearized q: variational equation becomes , fin , ,7 1 w b.“.( [V 6u(,,,c:,k,Au(,,,,dv+ [V 6m,,a;,Au,,,dV = 1 0611, v . [F 61“,:th — [V (9.2:; aijdv (3'37) \ fl" \\ 3.6 Cauchy Stress Update \ (\ ., 1 J 1 In the previous sections, both displacement and stress components are expressed as incre- J mental forms when a consistent linearization approach is taken to approximate the varia- A) tional equations. As discussed in Chapter 2, the Cauchy stress rate. cannot be calculated directly because it is not objective. On the contrary, the Truesdell rate Of Cauchy stress can be evaluated via a suitable constitutive equation. Certainly, this cannot be accomplished J until the displacement increment is identified. Once the Truesdell rate is obtained, Eq.(2. 16) can be used to calculate the Cauchy stress rate. It is noted that both Cauchy stress and Cauchy stress rate are referred to a fixed reference frame during deformation. Therefore, once the Cauchy stress rate is obtained, it is to be added to the Cauchy stress that has been accumulated over previous increments to become the current Cauchy stress. 34 3.7 Newton-Raphson Iteration The Newton-Raphson iteration scheme is used for incremental solutions. A step-by-step solution procedure using the equations presented in the previous sections is described below. 1. In the beginning of current increment, a set of trial displacement increments Au, is assumed. 2. The Truesdell constitutive equation is calculated using Eq.(3.29). 8' l i - 3. The rate-of-deformation and velocity gradient tensors are identified using Eqs.(3.26) and (3.25), respectively. 4. The Cauchy stress rate is evaluated using Eq.(3.30). 5. Every term in Eq.(3.37) is then calculated. 6. If the residual (right-hand side of Eq.(3.37)) satisfies the convergency criterion (see Appendix D), update the Cauchy stress by using of,“ = of, + Ao,,-, and start the next increment. 7. If the convergency criterion is not satisfied, solve the system of simultaneous equations, Eq.(3.37), for a displacement correction vector Au? 8. Let An? = Am + Auf. 9. Repeat the steps starting from Step 3 by using Au? as the new displacement incre- ment, instead of Au... Chapter 4 INCREMENTAL DISPLACEMENT FIELD 4.1 Introduction An incremental displacement field based on the Generalized Zigzag Theory is presented in this chapter. It needs to be established before the rate-of-deformation and velocity gradient tensors can be evaluated. The major advantage of the Generalized Zigzag theory as described in Chapter 1 is that its variables are independent of the total number of layers for a composite laminate. The theory was originally proposed by Li & Liu [38] in studying infinitesimal deformation of a composite laminate. In order~ tobe used for large deformation '“‘ fi‘h». analysis, the linear displacement field must be written in an incremental form. The major “...—0.,” difference between the infinitesimal deformation and the large deformation analyses is due to ) the fact that, in discussing infinitesimal deformation, the spatial coordinate system always coincides with the material coordinate system in a continuum. However, in dealing with ‘0“ large deformation, the two coordinate systems must be described separately. 35 36 4.2 Displacement Field Shown in Figure 4.1 for a continuum, when a particle at position P at time t is deformed to a new position P’ at time t + At, components of the displacement increment of the particle at time t can be expressed as follows. Auka, y, z) = AU5(I.y) + Ania-.3112 + Amt-.1112? + mm. 11123 (4.1) Avk1x,y.z) = Avie. y) + Avie. y). + Am, 11122 + August/12‘ (4.2) Aw"(ar,y,2) = Awo(:r.y) (43) where 2:, y, z are spatial Cartesian coordinates of the particle P at the instant of time t, k is a layer-number index, where the bottom layer corresponds to k = 1, Aug, A123 and Awo are translational displacement components in 2:, y, and 2 directions, respectively, Au’f and Av]c are first-order rotational displacement components about y and z axes, re- spectively, and Au,- and Av,- (i = 2, 3) are higher-order rotational displacement components about y and z axes, respectively. 4.3 Displacement Continuity Conditions Imposing displacement continuity condition at each laminate interface, we have, for let” ml layer where It = 2, ..... , n with n being the total number of layers, Au"—1 Au’c z=zk azzqE Auk"1 = Avk z=zk z=zk Awk“ = Aw" = Awo (4.4) 2:74,, 37 By using Eqs.(4.1),(4.2) and (4.3), Eqs.(4.4) can be written as Au],E — Au];—1 2 (Au'f'1 — Au’f)zk A116c — A116"l = (Av’f—l — Avflzk If we let Au], = Auo and A116 = Avo Eqs.(4.5) become, for k=2,3, ..... , n, k . . A116c = Auo + X:(Au{-1 — Au{)z,- i=2 k . Av],c = Avo + X:(Av{"l — Av{)z,- i=2 Now let’s introduce the kinetic relations for the k“ layer of the laminate, r. " r L); c /-.. .i ._ Ur, . _ i {011:1} = [0"]{0} where T k k k I: I: k k {(4 >} = {01.1.01 1,.1 instead.) . [0"] is the constitutive equation, i.e. ( r . Ch 012 01:. 0 0 Cfs sz C22 053 0 0 056 [C k] _ Cfs 053 Cis 0 0 Cit 0 0 0 Cf, 0 0 oooocgso CfG 62:6 03:6 0 0 C56 l f f I I, .’ Au!) fan“; 6‘ " (4.7) (4.3) (4.10) ( 38 and {D} it the rate-of-deformation tensor, - 8Auk/83: BAvk/ay BAwk/Bz {D} = l 1 (4.11) 6A1)" / 32: + aAwk/By BAu" /6z + BAw" / 8:1: 1 BAuk/By-l-dAvk/Bz ) Remark 4.1 Later in the development of a finite element, the (2:,y,z) coordinate system will be so chosen that it coincides with the coordinate system of each element. The coordinates (1:,y,z) in general do not coincide with the material principal directions for each lamina. Remark 4.2 In Eq.(4.10), 025 = 0 implies that the out-o ~plane shear moduli, G13 and 053, are equal. As can be seen in the following derivation, when interlaminar shear stress continuity conditions are imposed, it will cause strong coupling between in plane displacement components if 0:5 7‘ 0. 4.4 Interlaminar Shear Stress Continuity Conditions Impose the following interlaminar shear stress continuity conditions: k-l k T]: 1 = T521 2:33 2:21: (4.12) 39 By using Eqs.(4.8), (4.9), (4.10), and (4.11), Eqs.(4.12) become 1' off 0 Bank-1 /32 + aAwk-l/ay 0 egg-1 aAuk-l/az+aawk~l/ax .-. - -k p CL 0 6A1)“ / 62 + BAwk/By 0 0:5 6Auk/8z + BAw" / 6:1: h z=Zk By expanding Eq.(4.13), we obtain, for k = 2, 3, ..., n, Cf4Av’f - CfflAvl '1 = 29kzkAv2 + 39kngv3 + OkAwW CfisAu’lc — Cg‘s'lAul'f‘1 = ZkakAuz + 3kagAu3 + QkAwm where Q), = 0:51—01; and 9k = Cir-1‘ 0:4 Giving the following new definition, Au] = Aul and Av,l = A111 Eq.(4.14) can be rewritten as .1, gr - (68 [(5 :1? , Au’f = FfAm + F§Au2 + F§Au3 + FiAwo.x Av'f = L’fAm + L§Av2 + L§Av3 + Limo... where, for k = 1, 1",1 =1, 1“?1 =0, F31=O, Ir",l =0 11:1, Lg=o, 131:0, L1=o (4.13) (4.14) (4.15) (4.16) (4.17) (4.18) 40 and, for k = 2,3, ..... , n, F," =-.a,,1n’,-l I." =0,,L’f“ F; = 111172 ‘1 + 2(Qk/C§5)zk L" = bkLg“ + 2(e,./C.1=,)z,c F; = 11,13?“ + 3(nk/c§5)zg Lf‘ = 310;“ + ask/clap; (4-19) Ff = akzrf-1 + (121/Cg.) L: = b.1211“ + (eh/Ci.) a), = ngl/Cg‘s, bk = Cir/054 4.5 Free Surface Shear Stress Conditions In most laminated composite studies, surface shear stresses are set free. By imposing the free shear stress conditions on the top and bottom surfaces of the composite laminate, we have =0 z=3n+l 1.35;) = 71(2) 2:21 = 0 (4.20) z=zn+l I T.](:Z) = Tag?) 2:21 By using Eqs.(4.8) and (4.16), Eqs.(4.20) can be expressed as ' 1 Du D12 [3112 [ E11 E12 A01 ' = (4.21) D21 D22 A113 L 5'21 5'22 A100,: . [ ) F11 F12 A02 H11 H12 A01 , = (4.22) F21 F22 A113 1. H21 H22 . Away where Du = 23n+1+ Ft? 012 1‘ 32,2”,1‘1- F? 021 = 2251 022 = 32% E11=—F1n E12=—(Fr+l) E21=—l. E22=—I (4.23) F11 = 22n+1+ L3 F12 = 3Z§+1+ L3 F21: 221 F22 = 32? H11=-L? H12=-(2+I) H21=—I H22=-I 41 Solving Eq.(4.21) and (4.22) for Aug, A113, Avg, and Avg, we obtain A‘UQ = AlAul + AgAwoJ, A113 == 81AU1+ BgAwOJ A1)? = C(Av1 + CQAw0,y, A173 = DlA‘Ul + DgAwoy (”J-24) where A1 = (471m(EllD22 - E21012) A2 = (35571512022 - E22012) B1: 73:7);(321011- E11021) B2 = (3.2117(322011‘ E12021) (det)1 = 011022 - 021012 7‘ 0 CI = 7%”:(3115'22 - ”211712) 02 = W1H12F22 -’ H22F12) D1 = W(H21F11- H11F21) D2 = min—,(HnFu - H12F21) (4802 = 17111722 - F21F12 9‘ 0 (4.25) Substituting Eqs.(4.24) into Eqs.(4.17), we have Au’f = RfAu1+R§Awmz Av’f = 013111 + ogawo, (4.26) Substituting Eqs.(4.26) into Eqs.(4.7), we have Aug = Aug + SfAul + S§Awo,z Avg .—. Avo+PfAvl +P2"Awo,y (4.27) where, for k=1, 121:1, 115:0, 0[=0, 021:0 (4.28) 511:0, 521:0, 103:0, P21=0 and, for k=2,3 ...... ,n, M=N+mm+a@,%=w+aw+efi 011: = L’f + 6'ng + DlLk, 0" = L: + 6'ng + Dng 42 k k 51:2(4-1—14). s:=}:( #13,), l=2 (:2 k k Pf = 2(01-1 — o‘,)-.~,, P; = 2(05-l — 0.3).;7 (4.29) l=2 l=2 Finally, substituting Eqs.(4.24), (4.26) and (4.27) into Eqs.(4.1), (4.2) and (4.3), we have ' . the final incremental displacement field , . ,, fix". Aukt’r. y. z) = Auo(r. y) + (51‘ + Riz + AIZ2 + 31231A"1($)y) + (537 + R52 + A222 + 8223) [W (4.30) Av"(a:, y, z) = Avo(:r, y) + (Pfc + O’fz + 0122 + Dlzs)Av1(:r, y) + (P; + Oé‘z + 02.22 + 0223) [W] (4.31) Awk(:r, y, z) = Awo(:z:, y) (4.32) It is noted from Eqs.(4.1), (4.2), (4.3), (4.30), (4.31) and (4.32) that the total number of variables in the incremental displacement field have been reduced from layer-number dependent to layer-number independent (seven variables in the final form). [Remark 4.3" When we use the free surface shear stress conditions, some errors will be introduced. The error may come from the following two aspects: 1. The applied load may not necessarily perpendicular to the laminate ’s top and bottom surfaces in the beginning of the loading history. 2. During deformation, the applied load may not be always perpendicular to the lami- nate’s top and bottom surfaces. However, if free surface shear stress conditions are not imposed, the total number of inde- pendent variables in the incremental displacement field will increase from seven to nine. 43 5x Figure 4. 1 Coordinate systems for a laminate shell Chapter 5 FINITE ELEMENT FORMULATION 5.1 Introduction The purpose of this chapter is to establish a solution procedure for general three-dimensional shell problems by spatially discretizing the governing equations derived in Chapter 3 via a finite element method. It is to transform the linearized variational equations into a system of linear algebraic equations by using the assembly of construction-identical element . - _ contributions. Once a linearization formulation is completed, the introduction of finite elements into the governing equations is merely a matter of selecting appropriate interpolation functions to approximate the unknown variables element by element. While the choice of appropriate interpolation functions and associated integration schemes, eg. full, reduced, or selectively reduced integrations, are still topics of many ongoing research endeavors, we do not attempt to settle the matter here. Instead, we start with a four-node, quadrilateral shell element 44 45 using a four-point integration scheme, and each node has seven degrees of freedom. In fact, the linearized governing equations we have formulated can be used with any new shell elements once they become available. Due to the dissimilarity between laminated composites and isotropic materials, a method called Zigzag Jacobian is proposed. The purpose of this technique is to emphasize the fact that the in-plane displacements of laminated shells are in a zigzag fashion through the laminate thickness. It is different from the Reissner-Mindlin theory, which assumes linear displacements through laminate thickness. Through the calculation of the Zigzag Jacobian matrix, kinematics of laminated shells can be expressed more precisely. Although the Zigzag Jacobian matrix requires nodal displacements at each element on the interfaces, it can be done in a postprocessing procedure. Hence, the total number of degrees of freedom for each element will not change. 5.2 Incremental Strains Recall Eqs.(4.30), (4.31), and (4.32) in the previous chapter, the incremental displacement equations can be rearranged as Auk = (Aug + S, Au1 + SkaAwo) + (RfAu1+aRgaA-—-w—o)z + (AlAu1+ .42 Mfg) +(B Au1+ 32 63:") 23 (5.1) Av" = (Avg + P,” Av) +P 66;”°)+ (01131;1 + 01332:“) 2 + (ClAvl + C26W)+(D(Av1 + 0266;00 23 (5.2) Awk = Awo (5-3) 46 Define a new vector {Ack } by rearranging the entries of the rate-of-deformation tensor in Eqs.(4.11) as BAuk/B'J: 6A1)" /6y v A C)! .Lb v {A‘Ek} = BAuk/dy + BAN/6:1: 3A1)" /32 + 613w" /8y t BAuk/az +0Awk/Bx J Combining Eqs.(5.1), (5.2), (5.3), and (5.4), it yields {Ask} = {ACO}+{A70}+{AK1}Z + {Arc2}z2 + {Arc3}z3 (5.5) where 0A:g + 55613:, + S; 628ng 2 %m+stfi%r+ss%%#+%m+Pr%m+Pé%-atm > (56) A {A60} 0 0 {A70} = 4 0 > (5.7) {M} {AKQ} {A53} A; \ 47 31:63:. + RI; ('92ng 0f 465‘; + 05 fiflygw an an 62; Rffit+offiu+(ag+0§)fi > 2C1Av1 + 2025—3551 2A1Au1+2Agi‘3—xm AlaA:l +1426;ng Clam), +C26263wg 211% +01% + (A2+02)Q3-$§n l 301Av1+3023~é§n 331m” + 332%? 310A!“ 4.320262% D1 6A0] +0232?“ 3175““ + 01 The: + (32 + Dz) #19sz 0 v 0 5.3 Incremental Displacement Gradients Incremental displacement gradients can be expressed as follows: 7 V 6Au£c /6:r v v (...) away 1 BAuf /62 J 1': 1,2,3 (5.9) (5.10) (5.11) where Let where V (Auk) {Aw0} {Awl} {sz} {A1273} 48 1 Aug Aw" aAu" /8:r aAu" /6y BAuk/az {Awo} + {Awl}z + {A5172}:a2 + {At.:73}z3 r = 4 f 121%!) + 35%?“ R’f 0m + 3562A? 1 2A1Au1+ 2.422%:3m J l A1135“ ”2%!“ 141%“ + fig-5W 1 381Au1 + 382%? 318A:l +826;ng 310A!“ +8262Awn 0 ‘ V 8Aun +Sk6Aul +Sk02‘9Aw 137"“ + Sfifi‘i + 35779;“sz R’fAul + #513? v v (5.12) (5.13) (5.14) (5.15) (5.16) 49 Similarly, BAvk/Bz V(A"k) = BAvk/By 6Avk/Bz {A1100} + {Awl }z + {A]z[MD,]dV+ / [MB]T[C‘<*)]z[MD,]dV V V + /V[MDI]T[C‘("’122[MDl]dV+ /V[MDg]T[C‘(k)]z2[MA]dV 58 + /V[MDQ]T[C‘("’]22[MB]dV+ /V[MA]T[C‘(“]22[MDgldV + L[MB]T[C‘(k’]z2[MDg]dV + [V[MD(]T[C‘(")]z3[MDg]dV + fv[MDg]T[C‘<"’]z3[MDl]dV+ /V[MD;;]T[C‘("’]23[MA]dV + /V[M03]T[C‘(")]z3[MB]dV+ /V[MA]T[C"(")]23[MD;;]dV + /V[MB]T[C‘("’]23[MD3]dV+ /V[MDI]T[C‘(")]24[MD;;]dV + [V [M03]T[C‘(’°)]z4[MDl]dV + [V[MDg]T[C‘(")]z‘[MDgldV + [v [MDg]T[C‘(")]25[MDg]dV + L[MD3]T[C‘(k)]z5[MDg]W + [V [M03]T[C"(")]25[M03]dV (5.63) (3) Geometric Stiffness Matrix Substituting Eqs.(5.48), (5.49), and (5.50) into the last three terms on the left-hand side of Eq. (5.59), components of the geometric stiffness matrix can be expressed as [ng‘m‘] = [V ([GU0]T + [00,]Tz + (002?; + [GU3]Tz3) [a] ([0110] + [GU1]z + [002122 + [GU3]23) av + [V ([GVolT + [GI/1]“: + [014)T22 + [GV;]T23) [a] ([GVO] + [014); + [oi/2122 + (016,123) W + /V[GWO]T[a]dV = [V [GU0]T[a][GUo]dV+ [v [GU0]T[a]z[GU1]dV + / [GU1]T[a]z[GUo]dV+ / [GU1]T[0]22[GU1]dV V V + / [GU0]T[a]z2[GU2]dV+ / [GU2]T[a]22[GUo]dV V V + j;[GUo]T[a]z3[GU3]dV + /V[GU3]T[a]z3[GUo]dV + / [GU1]T[a]z3[GU2]dV+ / [GU2]T[0]23[GU1]dV V V + / [GU2]T[a]z4[GU2]dV+ / [GU,]T[a]z4 [0031.117 V V 59 + flamflpwamwiw /V[GU2]T[a]zS[GU3]dV + /, ,[GU31T1rrlzslGUgldv + /V[GU3]T[0]26[GU3]dV + /V[GV0]T[0][GVo]dV+ /V[GV0]T[a]z[GV(]dV + / [CV1]T[a]z[GVo]dV+/ [GVO]T[a]22[GV2]dV V V + / [0V2]T[a]22[GI/b]dV+ / [GVl]T[a]22[GV(]dV V V + /V[GVo]T[a]z3[GV3]dV+ /V[GV3]T[a]z3[GI/b]dv + / [GVI]T[a]z3[GV2]dV+ / [GVg]T[a]z3[GV1]dV V V + / [GVI]T[a]z4[GV3]dV+ / [Gv3]T[a]z4[Gv,)dv V V + / [GI/2]T[a]z‘[GV2]dV+ / [GVg]T[a]z5[GV3,]dV V V + / [GVg]T[a]z5[GV2]dV+ / [0V3]T[o]z6[GV3]dV V V + [v [GWO]T[0][GWo]dV (5.64) (4) External Force Vector The first term on the right-hand side of Eq.(5.59) can be written in the following form: [r {6115} {h}dr‘ = [S (aukh, + 50th, + 6wkhz) 115‘ (5.65) where h;, h,,, and h, are tractions in the element 1:, y, and 2 directions, respectively, and S‘ is current top and bottom surface areas. By using Eqs.(5.29), (5.30), and (5.31), components of the external force vector become {Fm} = / (({xvo}T + {XU,}Tz + {2(a)}T 22 +{X113}T 23) h, S. + ({XVo}T + {XV1}T Z + {XV2}T 22 + {XV3}T 23) h” + {XW0}T 5.):15' (5.66) >5 i T Bilinear Shape Functions 1v. -,i-(l+§a§)(l+nan) Hermite Cubic Shape Functions 21.11 -.',-<1+§.:><1+n.2n)<2+4.§+n.n-§2-n21 H02 = %§a(l+§.§)(t+nan)(§a§-l) ’ 2 95.3 ( 7‘n5(1+§.§)(1+nan)(rtn-l ) Figure 5.1 Shape functions ‘Z 61 .4 """" ” k=2 6: .y: .5 1|“ (element coordinates) x Figure 5.2 Illustration of element and natural coordinate systems for a typical element Chapter 6 NUMERICAL STUDIES 6.1 Introduction Following the formulation of using the Generalized Zigzag Theory for large deformation analysis in the previous chapter, a finite element program named LACOS (LAminated COmposite Shells) was programmed using FORTRAN. It was then linked to the commercial software ABAQUS/Standard as a new addition to its element library. This particularly user-defined element was called ”U101”, which followed the naming convention required by ABAQUS/ Standard. Several case studies are performed in this chapter. The case studies are conducted not only to validate the U101 element but also to evaluate the new element in various applications. The entire chapter is divided into two major parts: linear solutions and nonlinear solutions. Linear solutions are acceptable in many structural mechanics problems. Linear problems that have analytical solutions play a primary role in the validation of a finite element for- mulation based on a new theory. These linear problems not only help identify programming errors but also help filter numerical problems such as shear locking. 62 63 Nonlinear solutions in laminated composite analysis occur when structures are either subjected to large deformation or made of special, laminate stacking-sequences such as unsymmetric and angle—ply laminates. In the second part of this chapter, both experimental and analytical cases resulting in nonlinear solutions are examined. All the calculations of study cases were performed on a HP/ C200 workstation with double-precision arithmetic. 6.2 Cross-Ply Laminate Under Cylindrical Bending As discussed in Chapter 4, the present work is based on the Generalized Zigzag Theory (ZIGZAG). Since there exists an exact solution for a simply supported, cross-ply laminate subjected to sinusoidal transverse-pressure, it is the first objective of this study to compare the theoretical solutions from ZIGZAG with those from other laminate theories. The linear solutions were presented by Pagano [53] by solving the problem based on a plane strain assumption. The results are expressed in terms of displacements and stresses through the laminate thickness. Among the various laminate theories available, the Interlaminar Shear Stress Conti- nuity Theory (ISSCT) [32] represents a class of layerwise theories with the total number of unknown variables dependent on the number of layers. It therefore requires very high computational time in solving the variables by finite element or other numerical methods. A third-order shear deformation theory presented by Lo, Christensen and Wu [44, 45] represents a family of ”global” (instead of layerwise) theories that use high-order polynomi- als to describe the displacement fields through the laminate thickness, namely High-order Shear Deformation Theories (HSDT). Like ZIGZAG and HSDT, theFirskorderSIiéTar De- formation Theory (FSDT)_also has five unknown variables and has been now commonly _. A c‘.u.—-v .1.‘ I used in commercial software such as ABAQUS, LS—DYNA3D, PAM-CRASH, and RADIOSS 64 CRASH. The theoretical solutions to the Pagano's problem based on the ZIGZAG, ISSCT, HSDT, and FSDT can be obtained by followingthe Navier solution procedure using double- . adm- ‘ " ~n”-or—a~~"‘"" Fourier series. u “ '5. “ta—... w.“ »-° A simply supported, [0/90/0] laminate under cylindrical bending is taken as a case study. The distributions of in-plane displacement is shown in Figure 6.1 for various laminate theories. The composite laminate of the length-to-thickness ratio equal to 4 is used in this study. As shown in the figure, the in-plane displacement distributions are measured at the end of the laminate. The values are normalized by a factor involving elastic constants, thickness of laminate, and the intensity of applied pressure. The kinky patterns from ZIGZAG and ISSCT are clearly illustrated and match with the elasticity solution closely. The solution from HSDT is a smooth cubic polynomial curve deviating from the elastic solution. Also shown in the figure is the result from FSDT, which is only a straight line of a constant slope. The distribution of the in-plane normal stress through the laminate thickness is depicted in Figure 6.2. Again, solutions from ZIGZAG and ISSCT follow the elasticity solution closely, while HSDT gives much smaller values near the interfaces and FSDT detours to the opposite side of the elasticity solution. It is noted that due to dissimilar material properties (in the sense of different fiber orientations) between adjacent layers, the in-plane normal stress distribution is not continuous through the thickness. The distributions of the transverse shear stress at the laminate end is shown in Fig- ure 6.3. In this case, ISSCT still matches the elasticity solution very well, while ZIGZAG gives a continuous distribution with certain inaccuracy. However, the distributions of the transverse shear stresses from both HSDT and FSDT are discontinuous through the lami- nate thickness. or“, 65 It is clear from the above examples that ISSCT and ZIGZAG present reasonably accurate results when compared to the linear elastic solutions. However, the results from HSDT and FSDT are inaccurate, especially in the transverse shear stresses. Although ZIGZAG is less accurate than ISSCT, it becomes clear that ZIGZAG is a better choice when considering both numerical accuracy and computational efficiencies. Remark 6.1 As the ratio of length-to-thickness of a laminate increases, e. g. a thin laminate, the difl'er- ence of transverse shear stress between ZIGZAG and elasticity becomes insignificant. As an example for the above statement, the maximum error for ZIGZAG is 5.4% (shown in Figure 6.4) for the ratio of length-to-thickness equal to 10, while 30% for the ratio of length-to-thickness equal to 4. The above example of cylindrical bending is investigated again by using U101 elements with a 10x1 mesh. The results in displacements and stresses are compared with the theoret- ical solutions from ZIGZAG. The purpose of this study is to evaluate the U101 element and to estimate its feasibility in solving more complicated structural problems where theoretical solutions can not be obtained. Shown in Figure 6.5 are the in-plane displacements. Excellent agreements exist between the finite element results at the integration points and the corresponding theoretical solu- tions of ZIGZAG. The finite element solution and the theoretical solutions are expressed by filled circles and a solid line, respectively. Excellent agreement can also be found in Figures 6.6, 6.7 and 6.8 for comparisons of in-plane normal stress, transverse shear stress and transverse deflection, respectively. 66 [0/90/O] cross-ply laminate S=L/h =4 L= 254.0mm=10.0inch E11=6.96Pa=1x106p31. n 1522:1725:GPa-25x10° psi q=qosm(T) 612:3.45:GPa -0.5x10° psi 11 Gza=1.=386Pa 0.211106 psi [— V12 = V23=0.25 x - "-555 ------------ :—» 2:2/11 U ___ E22U(O,Z) L b hqo Figure 6. l In-plane displacement distribution of a S=4, [0/90/0] laminate under cylindrical bending 67 Elasticity b [0/90/0] cross-ply laminate S=L/h=4 L : 254.0 mm : 10.0 inch 511: 6.9 GPa :1x10° psi 1522 : 172.5 GPa : 25x10° psi G12 : 3.45 GPa : 0.5x10° psi 6123 : 1.38 GPa : 0.2x10° psi V12 = v23 = 0.25 10. 15. 20. q = qosin(-’°1_‘-) "-fifi ............. ..J L > Figure 6.2 In-plane normal stress distribution of a S=4, [0/90/0] laminate under cylindrical bending / - —— Elasticity 4' .L /’ I ’4 .4 .8 1.2 1.6 2.0 2.4 [0/90/0] cross-ply laminate S = L/ h = 4 +2 L : 254.0 mm : 10.0 inch E11: 6.9 GPa :1x10° psi 1522 : 172.5 GPa : 25x10° psi 0 = qosin(%) 6,2 : 3.45 GPa : 0.5x10° psi Gm : 1.38 GPa : 0.2x10° psi .1." V12 = V23 = 0.25 X _ “—5545 ———————————— .—: Z=Zlh _= txz(o'z) L 5 112 QO «I Figure 6.3 Transverse shear stress distribution of a S=4, [0/90/0] laminate under cylindrical bending 69 [0/90/0] cross-ply laminate S = L/ h = 10 ‘2 L = 254.0 mm = 10.0 inch E11 = 6.9 GPa =1x10‘ psi 1522 = 172.5 GPa = 25x10° psi 0 = qosin(%) 612 = 3.45 GPa = 0.5x10° psi 623 = 1.38 GPa = 0.2x10° psi in V12 = V23=0.25 x "W ------------ ——» L N isz/n __ 1550.2) Clo d X2 Figure 6.4 Transverse shear stress distribution of a 8:10, [W90/0] laminate under cylindrical bending 70 o o U101 ZIGZAG -.5 .0 5 1.0 3 [0/90/01 cross-ply laminate Z S=L/h=4 ‘ i. = 254.0 mm = 10.00iench En = 6.9 GPa =1x1 psi . 101 522 = 172.5 GPa = 25x10° psi q ' q°s'"(‘|'-') 6.2 = 3.45 GPa = 0.5x10" psi 623 = 1.33 GPa = 0.2x10° psi v12 = v23=0.25 E=zln U = E220(0.Z) th Figure 6.5 Comparison of U101 and the Generalized Zigzag theoretical solutions of in-plane displacement distribution 71 .4 _ o o U101 ZIGZAG [0/90/01 cross-ply laminate S = L/ h = 4 A2 L = 254.0 mm = 10.0 inch En = 6.9 GPa =1x10° psi 1522 = 172.5 GPa = 25x10'3 psi q = qosin(l°‘._—) (3.2 = 3.45 GPa = 0.5x105 psi 62;; = 1.38 GPa = 0.2x10“ psi in v12 = v23 = 0.25 fifi x _ --— ------------ q——> Z Z/l'i 1 -:O(L/212) L h Figure 6.6 Comparison of U101 and the Generalized Zigzag theoretical solutions of in-plane normal stress distribution .4 \\ .. \ \ -.2 / -.4 7/ / 72 .0 .5 [0/90/0] cross-ply laminate S = L/ h = 4 L = 254.0 mm = 10.0 inch E11 = 6.9 GPa =1x10‘3 psi 1522 = 172.5 GPa = 25x10“ psi (3.. = 3.45 GPa .-. 0.5x1da psi 623 = 1.38 GPa = 0.2x10‘3 psi V12 = V23 = 0.25 i=z/h 2.0 t!!;(0,Z) Qo XI 0! Figure 6.7 Comparison of U101 and the Generalized Zigzag theoretical solutions of transverse shear stress distribution 2.5 73 1 4 " 1| 1| .2 , | 1| 0 o O U101 ' ZIGZAG 1| 1| . 2 " 1| 1| -.4 1| 1 -3 5 -3 0 -2 5 -2.0 -1 5 -1 0 - 5 [0/90/0] cross-ply laminate s = L/ h = 4 A2 L = 254.0 mm = 10.0 inch E11: 6.9 GPa =1x10‘3 psi 522 = 172.5 GPa = 25x10a psi q = q same) 812 = 3.45 GPa = 0.5x10° psi ° L 623 = 1.38 GPa = 0.2x10° psi V12 = V23=025 E=z/h _ 100x E22xh3xW(L/2.Z) L ”H W: 'l QOXL4 Figure 6.8 Comparison of U101 and the Generalized Zigzag theoretical solutions of transverse deflection 74 6.3 Bending of Rectangular Laminate The finite element program is next applied to a square plate made of isotropic material and subjected to uniformly distributed transverse-pressure. Due to the biaxial symmetry of the structure, only one quarter of the plate was modeled and symmetric boundary conditions were imposed. Three different mesh densities (2x2, 4x4, and 8x8) were investigated for convergency. Also, three length-to—thickness ratios (S = a/h = 5, 10 and 100) were exam- ined for she-ailoclging. The normalized, central deflection is summarized in Table 6.1. The material properties and the method of normalization are also given in the table. ABAQUS/S4R [2, 3] is a four-node, doubly curved, reduced integration shell element with'hourglass control and is based on the First-order Shear Deformation Theory (FSDT). The element also considers finite membrane strain and thickness change. As shown in the previous study, elements based on FSDT do not have correct displacement and stress cal- culations through the laminate thickness, especially in—plane displacements and transverse shear stresses. When reasonable transverse stresses are desired, the ABAQUS/S4R element uses in-plane stresses and equilibrium equations to recover them. The four-node, 28 degree- of-freedom element formulated by Palazotto and Dennis [55] is a fully integrated element based on a fourth-order shear deformation theory. Also shown in the table are the closed- form, Navier solutions generated by Reddy [62]. Reddy’s results were based on a third-order shear deformation theory. Again, it should be pointed out that recovering techniques for transverse shear stresses through equilibrium equations have been used by Palazotto and Dennis and Reddy. It can be seen from Table 6.1 that the difference of transverse deflection between U101 and Reddy’s closed-form solution is within 1% when an 8x8 mesh is used. Based on a study of various length-to-thickness ratios, U101 does not present any shear locking when the thickness of laminate decreases. 75 As a second example, the same structure is made of a zero-degree, orthotropic laminate. By using an 8x8 mesh, the normalized, central deflection and stresses of three length-to- thickness ratios are shown in Table 6.2. The material properties and normalization factor are also given in the table. The three-dimensional elasticity solutions were given by Srinivas and Rao [67]. Reddy's closed-form solutions [62] were also given for comparison. In general, the solutions from U101 agree very well with the three-dimensional elasticity results. The difference of the transverse deflection between U101 and three-dimensional elasticity is within 1%. The results from the Classical Laminate Theory (CLT) are also given in the table. It is shown that error increases as the length-to-thickness ratio increases because of the negligence of the shear deformation effect. In other words, CLT is not suitable for thick composite laminates. The in—plane stress 0,, is taken at an integration point near the center of the laminate, while the transverse shear stress 0w is chosen at an integration point near the laminate boundary. The stresses from U101 do not perform as well as the transverse deflection when compared with the three-dimensional elasticity solutions. It is believed that mesh refinement near the stress sampling locations is required because stress gradients around these regions are very large, especially when the laminate is thin [55]. The difference is also believed to be due to the fact that the stresses from U101 are recorded at integration points, not exactly at the boundary. The next case deals with a simply-supported [0/90/0], rectangular laminate under si- nusoidal transverse-pressure. The normalized transverse-deflections at the center of the laminate are summarized in Table 6.3 for four different length-to-thickness ratios. Also given in the table are material properties, loading conditions, and the normalization factor. The threedimensional elasticity solutions were given by Pagano [54]. An 8x8 mesh was used in the finite element solutions based on U101, ABAQUS/S4R and Palazotto and Den- 76 nis [55]. It is observed from Table 6.3 that U101 outperforms ABAQUS/S4R and Palazotto and Dennis, and even Reddy’s closed-form solutions at small length-to-thickness ratios. The superiority diminishes when the laminate becomes thin. Table 6.1 Normalized, central transversedeflection of a simply supported, isotropic, 77 square laminate under uniform transverse-pressure S=a/h Mesh U101 ABAQUS/S4Fi Palazottoai Reddy' Dennis 2x2 0.0314 0.0450 0.0240 100 4114 0.0430 0.0448 0.0425 8x8 0.0441 0.0444 0.0443 0.0444 I 2x2 | 0.0428 | 0.0478 0.0481 I 10 | 4x41 0.0458 | 0.0488 0.0488 | ] 8x8] 0.0484 | 0.0488 ] 0.0487 | 0.0487 | 2x2 ] 0.0500 ] 0.0549 | 0.0538 | 5 | 4x4 | 0.0525 ] 00539 | 0.0538 ] | 8x8 | 0.0533] 0.0538 J 0.0538 ] 0.0535 "’ theoretical solutions A y Simply supported 1‘ ''''''' t' ''''' ‘ isotropic material : ] E = 68.9 GPa u. l '3 8 . °i ' 1: =10.0x10 psi 5, la v -03 8' ' a ‘ ' 3' symmetry ] g alb=1.0 .o "l . J > , _>.~ I 2‘ x uniform pressure=q g g; l g normalized deflection = ] ’17) g : '6 W,,l13E/qa4 -‘ g ] l : simply sharia-n80- "- a 78 Table 6.2 Normalized deflection and stresses of a simply supported, orthotropic, square laminate under uniform transverse-pressure S = a / h 20 10 7.143 U101 10394 686.7 190.6 ABAQUS/S4Fl 10486 691 .4 191 .9 W Palazotto & Dennis 10448 889.8 191.7 Reddy‘ 10450 689.5 191.6 SD Elasticity' 10443 688.6 191.1 CLT' 10250 640.7 166.8 U101 138.9 34.34 17.47 ABAQUS/S4Fl 142.8 35.34 17.81 c Palazotto & Dennis 144.6 36.12 18.41 Reddy‘ 144.3 36.01 18.34 3D Elasticity' 144.3 36.02 18.35 U101 10.23 5.11 3.63 ABAQUS/S4R 10.27 5.12 3.64 1: Palazatto 8. Dennis 8.66 5.07 3.70 Reddy' 10.85 5.38 3.81 30 Elasticity' 10.87 5.34 3.73 " theoretical solutions orthotropic material [0] 5,, = 143.82 GPa = 20.83 x 106 psi 52 = 75.43 GPa = 10.94 x106 psi V12 = 0.44 8,2 = 42.08 GPa = 8.1 x 10° psi 1313 = 25.58 GPa = 3.71 x 10‘ psi G23 = 25.58 GPa = 3.71 x 10° psi 8 / b = 1.0 n 8 x 8 mesh uniform pressure = q normalized results: W = 23.2x106W0/qh o = oy (0,0.h/2) / q 1: = tyz (0,-b/2, 0) / q simply supported a Table 6.3 Normalized, central transverse-deflection of a simply supported, rectangular [0/90/0] laminate under sinusoidal transverse-pressure S = a/ h U101 ABAQUS/S4Fl Palazotto 8 Reddy‘ 30 Elasticity" Dennis 4 2.742 3.164 2.647 2.641 2.82 10 0.912 0.936 0.864 0.862 0.919 20 0.604 0.613 0.595 0.594 0.610 100 0.503 0.509 0.507 0.507 0.508 " theoretical solutions [0/90/0] laminate E1, = 258.8 GPa = 37.5 x 10" psi 1 y 522: 10.3 GPa = 1.5x 10'3 psi v12 = 0.25 simply supported 612 = 5.2 GPa = 0.75 x 106 psi {' ““““ 613 = 2.1 GPa = 0.3 x 10° psi I E‘ 8 . '0 l o ‘3 623:2.1 GPa=0.3x 10 ps1 g. E 1: a = 203.2 mm = 8.0 inch 8:] g. 5%. b = 609.6 mm = 24.0 inch 3 I a > 3 X 3 mesh fil symmetry ] 2'5 x sinusoidal pressure = ,5] . g q,sin[u(x+a/2)/a] sin[n(y+b/2)/b] : : normalized deflection = l ] 100w,h352/ qoa‘ l _________________ i simply supported a 80 6.4 Shear Locking Shear locking is a condition in which the element stiffness is immensely overestimated for any reasonable mesh density, thus yielding near-zero solutions [68]. It is observed that some fully integrated, Co shell elements that are based on the First-order Shear Deformation Theory are vulnerable to the shear-locking phenomenon [80]. Using a straight beam as an example, the deflection of a thin beam should only be governed by the bending stifl'ness matrix because the transverse shear deformation is neg- ligible. However, in the case of shear locking as the beam becomes slender, the traverse shear stiffness matrix grows in relation to the bending stiffness matrix as a result of enforc- ing the constraint of zero transverse shear strain. Accordingly, the traverse shear stiffness matrix acts as a penalty function that yields near-zero transverse-deflection. Even if more elements are used in this case, every additional element usually brings equal amounts of degrees of freedom and constraints when the full integration scheme is used. Therefore, using more elements will not help alleviate the problem. The proceeding arguments can also be expanded to plate problems. The U101 element does not have a shear-locking problem. The key reason is that the Hermite cubic shape functions are used in the U101 element to describe the transverse deflection 10. Using the straight beam case as an example again, the U101 element has Aw, Aw ,3, and A111 as degrees of freedom at each node; hence, the constraint of transverse shear strain is imposed on A111,; and Aul but not directly on the transverse deflection Aw. The aforementioned case presented in Table 6.1 was used again to verify, to some extent, this claim. The new length-to—thickness ratio was chosen as 1000 ( =a/ h) for the U101 element via an 8x8 mesh, while the other conditions were kept the same. The result from the U101 element is 98.8% of the' theoretical solution by Timonshenko and Woinowsky-Krieger [71]. 81 6.5 Patch Tests The patch test was first given by Irons [9, 25]. Many authors have emphasized the impor- tance of this test [8, 14, 21, 22, 79]. The consensus is as follows. 1. Passing the patch test is a sufficient condition to convergency but not a necessary condition. Under a convergent situation, the behavior of the real structure can be reproduced as the size of element decreases. 2. It is important to use irregular element shapes in constructing a patch because one 6- nite element formulation may pass the patch test in certain special mesh configuration but not in others. In the patch test, the elements are assembled in such a way that at least one node is completely surrounded by elements. In performing the test, a displacement field that provides an arbitrary state of constant strain, or a consistent nodal loading, is applied to the boundary nodes. Nodes not on the boundary are neither loaded nor restrained. Solutions for degrees of freedom of all nodes that are not prescribed are sought, and the strains (or stresses) in all elements are computed. The element passes the patch test if the resulting displacements at the internal nodal points correspond with the applied displacement field and the computed strains and stresses at every point in every element agree with analytical values to the limit of computer accuracy [14, 21]. As given in Table 6.4, the U101 element passed the membrane patch test. However, the U101 element did not pass the bending patch test when using the same mesh configuration. The failure in passing the bending patch test is believed to be caused by the mixed approach used in the U101 element [55]. That is, the U101 element uses the nonconforming Hermite cubic shape functions for the variables of Aw, AwJ, and Awm while bilinear shape functions 82 for the variables of Au, Av, Aul, and A01. The reason for using the mixed approach is because Aw needs C l continuity, while Au, Av, Am, and A01 need only C'0 continuity. In the bending patch, a pure bending mode is created, either by prescribing rotational degrees of freedom or by applying bending moments. It is impossible to have zero transverse shear strain when for example Aw; and Aul are interpolated differently. One way the problem may be conquered is to reduce the continuity requirement of Aw from C1 to C0 by artificially imposing a certain type of constraint on transverse shear strains, e.g. ”discrete Kirchhoff” [21]. The consequence of passing the membrane and failing the bending patch tests can be also viewed from the convergency point. The convergency is guaranteed for the U101 element under a membrane loading; however, it may or may not converge under bending. As can seen in the previous cases and later examples, the U101 element has not shown any unconvergency problem under bending conditions. Table 6.4 Results of membrane patch test coorcflnates Woretical soiutfins U101 solutions node (x,y) (u,V)x10'° (U,V)x10'3 . 1 (0.04, 0.02) (005? 0.04) (0.05, 0.04) 2 (0.18, 0.03) (0.195, 0.120) (0.195, 0.120) 3 (0.16, 0.08) (0200,0180) (0.200, 0.180) 4 (0.08, 0.08) (0.120, 0.120) (0.120, 0.120) linear, elastic material 5,, = 1.0 x 10° V12 = 0.25 h = 0.001 at all exterior nodes: 0 = 10°(x + y/2) v = 10'°(y + x/2) at all nodes: W = 0 AV 0.12 ] 1 T . 0.24 84 6.6 Large Deflection Effects in Unsymmetric Cross-ply Com- posite Laminates Due to the existence of bending-stretching coupling, nonlinear deflection occurred even when an unsymmetric composite laminate is under a small loading. Sun and Chin [70] presented some theoretical solutions for an unsymmetric cross-ply laminate subjected to uniform transverse-pressure. Their analysis was based on the Kirchoff-Love hypothesis, in which the effects of transverse shear deformation were not taken into account. They also used a large deformation theory in you Kc’irmén sense. Accordingly, a special, linear ordinary differential equation with constant coefficients was solved using a plane strain assumption for a. laminate under a pin-pin boundary condition. The relation between the uniform pressure and the normalized, maximal transverse- deflection is shown in Figure 6.9 for a length-to—thickness ratio 5 of 225. The uniform pressure ranges from 0 to 68.95 kPa (10 psi). The theoretical solution of Sun and Chin is illustrated as a solid line. The solutions from U101, ISSCT (Interlaminar Shear Stress Continuity Theory [32]), and ABAQUS/S4R are all based on the finite element analysis of a, 10301“ “mesh It is observed in the figure that the results are consistent with one another because the effect of transverse shear deformation can be neglected in the thin laminate. It is also noticed that the solutions from a positive loading and a negative loading are different. In addition, the in-plane total reaction force at the pinned end is calculated and shown in Figure 6.10. Again, because a large 3 of 225 is used, good agreements are found among all solutions. The same structure was subjected to higher pressure of up to 6.89 MPa (1000 psi). The relations between the uniform pressure and the normalized transverse-deflection are shown 85 in Figure 6.11. It is noted that the results split into two groups when the pressure is higher than 2.76 MPa(400 psi). The results from U101 and ABAQUS/S4R are in one group, while the rest are in the other. The largest difl'erence in transverse deflection between the two groups is about 2.6% when the pressure is 6.89 MPa (1000 psi). This result may be due to the use of different large deformation theories. It is noted that U101 and ABAQUS/S4R use an updated Lagrangian approach, while both ISSCT and Sun and Chin use a total Lagrangian approach with von Kérmdn nonlinear strains, which is a simplied version of the full-term Green strains. The difference among the curves from the negative loading is not noticeable. The next study investigates a similar structure with a smaller length-to-thickness ratio S, i.e. 11.25. The uniform transverse-pressure ranges also from 0 to 6.89 MPa (1000 psi). This loading is in fact small because the laminate is ”thicker” now. The purpose of this study is to examine the efl'ect of transverse shear deformation. It can be seen from Figure 6.12 that the transverse deflection from Sun and Chin is smaller than the rest of the theories that consider transverse shear deformation. The difference can be as much as 27.5% at a positive loading of 6.89 MPa (1000 psi). Some interesting observations can be found from the same case for the in-plane, total reaction force shown in Figure 6.13. When the loading is positive, the reaction force remains negative. Until the the loading reaches 3.45 MPa (500 psi), the reaction force starts to increase positively. It eventually becomes positive. A negative reaction force in this study actually indicates that the laminate pushes the end supports instead of pulling them. The last set of curves are given in Figure 6.14 in which a positive transverse-pressure of 68.95 kPa (10 psi) is applied to the same structure using various length-to-thickness ratios (S=L / h). The largest transverse-deflections from all finite element formulations are 86 normalized by the corresponding theoretical results from Sun and Chin. By doing so, the results of Sun and Chin are actually a horizontal line of 1.0 in Figure 6.14. Again, due to the effect of transverse shear deformation, larger transverse-deflections were found in U101, ISSCT, and ABAQUS/S4R for thicker laminates. However, these results approach Sun and Chin’s solution as the laminate becomes thinner. 87 91981) -10 -8 -6 -4 -2. 0 2 4 6 8. 10 5 J 4 —— Sun&Chin ' G- - 0 U101 A----A ISSCT 3- 13- - -a ABAQUS/S4R 2 Z 1 W 0 -1. -2. -3. -4. -5. -68.9 -55.2 -41.4 -27.6 -13.8 0.0 13.8 27.6 41.4 55.2 68.9 9009) [0/90] cross-ply laminate 42 E11 = 137.9 GPa -.- 20.0x10°psi 222 = 9.7 GPa: 1.4x 10° psi Giz=st=Gis= 4.8 GPa = 0.7 x 10°psi V12=0.3 . q h t:i?.%°m':‘."‘:.i.i."fi.. Hi llHlll i S=L/h=225 _-_]« 4] .............. mesh=10x1 W=Wnulh L w Figure 6.9 Small load vs. transverse deflection curve of a S=225, [0/90] unsymmetric laminate 88 9198') -10. -8. -8. -4. -2. 0. 2. 4 8. 8 10. 2670. -6 2225. .5 1780. .4 3 1335. .3313 X x z 890. .2 —— Sun&Chin 445-09 e-- «e U101 .1 A----A lSSCT B---El A8AQUS/54Fi 0. i I L 0 -68.9 ~55.2 -41.4 -27.6 -13.8 [0/90] cross-ply laminate E11 =137.9 GPa = 20.0 x 10° psi E2 = 9.7 GPa =14 x 10° psi Giz=cza=eie= 4.8 GPa = 0.7 x 10°psi V12 = 0.3 L = 228.6 mm = 9 inch h = 1.02 mm = 0.04 inch S = L I h = 225 mesh = 10 x 1 Nx = total in-plane force at supported and 0.0 13.8 27.6 41.4 55.2 68.9 9 (RP!) ¢\¢—- 3133111. Figure 6.10 Small load vs. in-plane force curve of a S=225, [0/90] unsymmetric laminate 89 q 11103035” -1.0 -.8 -.6 -.4 -.2 .0 .2 .4 .6 .8 1.0 25. f l l i :. 20.4 — SUD & Chin @- - e U101 1 A----A ISSCT 5" e---e ABAQUS/S4R 10. / -25. -6.89 -5.52 -4.14 -2.76 -1.38 0.0 1.38 2.76 4.14 5.52 6.89 9 (MP!) [0/90] cross-ply laminate 2 E11 =137.9 GPa = 20.0 x 10°psi 4 E22 =9.7GPa=1.4x10° psi Giz=st=G13= 4.8 GPa = 0.7 x 10°psi V12 = 0.3 0...”... L11 1.18.1.1}. l- 1" S=L/h=225 mesh=10x1 W=Wmlh L D Figure 6.1 1 Large load vs. transverse deflection curve of a S=225, [0/90] unsymmetric laminate 90 1: x103(psI) -1.0 -.8 -.6 -.4 -.2 .0 .2 .4 .6 .8 1.0 .6 l l l A _ — Sun&Chin ' , '5 e-e U101 X / A----A ISSCT . .4 ‘ l3---El ABAQUS/S4Fi [ -.2 -.3 -6.89 -5.52 -4.14 2.76 -1.38 0.0 1.38 2.76 4.14 5.52 6.89 9(MPa) ‘2 [0/90] cross-ply laminate E11=137.9 GPa = 20.0 x 10°psi E22=9.7GPa=1.4x10° psi G12=623=G13= 4.8 GPa = 0.7 x 10°psi 11:212.“... iii [111110 i“ h =20.32 mm :08 inch _ IX S=L/h=11.25 Kg] mesh=10x1 _ L i w=w /h max Figure 6. 12 Large load vs. transverse deflection curve of a S=11.25, [0/90] unsymmetric laminate 91 1: x103(psi) -1.0 -.8 -.6 -.4 -.2 .0 .2 .4 .6 .8 1.0 89.0 20. 71.2 18. 53.4 v 12. ’5‘ 35.8 8. i s V a '2, 17.8 4.2 3" x x 2‘ 0.0 0.2 47.8 — Sun & Chin _4. e- — e U101 15- - - -A ISSCT -355 (3- -—£i ABAQUS/S4Fi , -8. -53.4 * 12- -6.89 -5.52 4.14 -2.76 -1.38 0.0 1.38 2.78 4.14 5.52 6.89 9 (MPa) [0/90] cross-ply laminate 2 E11 = 137.9 GPa=20.0x 10° psi A E22 =97 GPa: 1.41110a psi 6 Giz=st=Gia= 4.8 GPa = 0.7 x 10 psi V12 = 0.3 L = 228.6 mm = 9 inch h = 20.32 mm = 0.8 inch 8 = L/ h = 11.25 mesh = 10 X 1 Nx = total in-plane force at supported and .lilifii 1" Figure 6.13 Large load vs. in-plane force curve of a S=11.25, [0/90] unsymmetric laminate 92 —--- Sun 8. Chin o- - o U101 A- - - -A ISSCT (3- - —-El ABAQUS/S4Fl Llh [0/90] cross-ply laminate E11=137.9 GPa = 20.0 x 10°psi E22: 9.7 GPa =14 x 10° psi G12=Gza=G13= 4.8 GPa = 0.7 x 10°psi V12 = 0.3 L = 228.6 mm = 9 inch S = L/ h W = max. deflection normalized by Sun 81 Chin’s solution Figure 6.14 Length-to-thickness ratio vs. transverse deflection curve of a small loading, [0/90] unsymmetric laminate 93 6.7 Comparison of Laminates Under Large Deformation with Experimental Data There are very few articles in the literature that investigate laminated composites under large bending-deflection using experimental approaches. The tests done by Zaghloul [74] and Kennedy [76, 77] in the early 703 presented some expermental data for symmetric and unsymmetric laminates subjected to large transverse-deflection, The symmetric and unsymmetric laminates were made of an epoxy matrix and unidirectional glass fabrics by a hand-lay-up process. Both clamped and simply-supported boundary conditions were performed in the tests. A uniform, transverse air-pressure was also applied to both square and rectangular laminates placed in an air chamber. The maximal transverse-deflection at the center of the laminate was measured with a mechanical dial-gauge. In addition to the experimental work, Zaghloul and Kennedy also presented a simple phenomenological method [74, 75] to determine the laminate constants from the properties of fiber and epoxy. The calculated constants were then employed in a set of governing nonlinear equations, which were then solved by a finite difference technique. The governing equations were based on a plain stress condition and also accounted for the Kirchhoff-Love hypothesis and von Karmdn nonlinear strains. Three studies were performed using U101 and ABAQUS/S4R with a 12x12 biaxially symmetric mesh. The first study was for a square [60]4 orthotropic laminate with clamped boundaries. The maximal transverse-deflections (normalized by the laminate thickness) and the setup of the study are shown in Figure 6.15. It is found that the transverse deflection from U101 is about 6% less than the experimental data at 13.8 kPa (2 psi), while ABAQUS/S4R is about 9% less. 94 The second study is for a square [0/90]2 laminate with simply supported boundaries. The results of the normalized, central deflection are depicted in Figure 6.16. In addition to U101 and ABAQUS/S4R, the numerical results obtained by Zaghloul using the aforementioned finite difference method is also presented. Again, all the numerical results tend to have less deflection than the experimental data. At 7.6 kPa (1.1 psi), the central deflection of ABAQUS/S4R is 6.5% less than that of the experimental data, followed by U101 8.0% less, and Zaghloul’s 13.6% less. As mentioned above, Zaghloul’s formulation was based on the Kirchhoff-Love hypothesis, which did not consider transverse shear deformation. The third study is for a square {—604/604] laminate with clamped boundaries. By following the same setup as used in Figures 6.15 and 6.16, the normalized, central deflections are shown in Figure 6.17, in which ABAQUS/S4R and U101 have 7.2% and 9.5% less deflections than the experimental data, respectively. In establishing the comparison between the experimental data and the numerical results, the testing setup and specimen fabrication should be carefully examined and their differences should be identified. Although Zaghloul [74] gave quite detailed descriptions regarding the testing setup and specimen fabrication, there were some unclear aspects. For example, there were not enough details to support that the boundary conditions used in the numerical analysis could truly represent those used in the experiments. Besides, the accuracy of air-pressure measurements in the experiment was unknown. With these issues in mind, the numerical results within 90% accuracy of the experimental data should be considered satisfactory. El 95 (1 (psi) .0 .5 1.0 1.5 2.0 4.5 1 — 0101 4.0 ----- ABAQUS/S4R o 0 Experiment 3.5 l. 3.0 y ' _ . o ........... 2.5 1* ,/" O / 2.0 1 / 2 y/ 1.5 7 1.0 / .Z ' 0.0 3.4 6.9 10.3 13.8 q (”’60 [60/60/60/60] laminate 6 y E11 = 20.68 GPa = 3.00 X 10 DSl clamped E22: 8.83c3F>a=1.28x106 psi --------- ...————————-—l 612:2.55 GPa=0.37x10° psi 1 l 1 G23=1.386Pa=0.20x10° psi : g : ea: 1.38 GPa=0.20x10° psi : 1 : v12: 0.32 E 1 mmet : 1 8 a=b=304.8mm=12inch ' 3y '3’ ' --------- J b E h=1.7mm=0.066inch g g. : g s = a / h = 182 0 o 1 ° uniform pressure = q E 1 mesh = 12 x 12 a; : l W = Wmax / h F a ’| clamped Figure 6.15 Normalized, central transverse-deflection of a [60l4 laminate under uniform pressure 96 Cl (psi) .0 .2 .4 .6 8 1.0 5.0 r l ----- ABAQUS/S4R 4 - - - Zaghloul '0 o 0 Experiment g ' .1 3 5 - ‘ ‘ x/ . //’ ’ i 3.0 */" o, , ‘ 7 /' ’ I W 2.5 ' ' . ’ I 2.0 , , ’ 1.5 ,. , r . / 1.0 '7 I a I // .5 ff f" 0/ .0 0.0 1.4 2.8 4.1 5.5 6.9 1: (kPa) [0/90/0/90] laminate y E11: 20.68 GPa = 3.00 x 106 psi simply Supported E22: 8.83 GPa=1.28x10° psi --------- T----------, (312 = 2.55 GPa = 0.37 x 106 psi : 1 l 823:1.388Pa=0.20x10° psi 8 . : l g G,3=1.3BGPa=0.20x10°psi t : : : t v12: 0.32 § l 5 met I . l é a=b=304.8mm=12inch (3' ym fl --------- -: be?) h = 1.9 mm = 0.078 inch 3 g: 1 E S = a/ h = 158 g g j g uniform pressure = q "’ E ] mesh = 12 x 12 «’1‘ : w = Wmax / h a +] simply Supported Figure 6.16 Normalized, central transverse—deflection of a [0/9012 laminate under uniform pressure 97 90280 .0 5 1.0 15 2.0 3.0 -——-- U101 2.5 ----- ABAQUS/S4Fl o 0 Experiment 2.0 W 1.5 1.0 . / ./ .0 0.0 3.4 6.9 10.3 [604/604] laminate y E11: 20.68 GPa = 3.00 x 10° psi T- damped E22: 8.83 GPa = 1.28 x 10° psi G112 = 2.55 GPa = 0.37 x 10° psi 883:1.38 GPa = 0.20 x 10° psi 6,3 =1.38 GPa = 0.20 x 10" psi v12: 0.32 a = b = 304.8mm = 12 inch 8 symmetry _________ .. h = 2.4 mm = 0.095 inch E a s = a / h = 128 ° ’5 uniform pressure = q E mesh = 12 x 12 5‘ W = Wmax / l‘l ‘ a > clamped Figure 6. 17 Normalized, central transverse-deflection of 8 [604/604] laminate under uniform pressure 13.8 clamped Chapter 7 CONCLUSIONS AND RECOMMENDATIONS 7 .1 Conclusions This thesis shows the problem-solving methodology improving the accuracy and efficiency of the finite element simulation of composite laminates under large deformation. An updated Lagrangian approach is employed to analyze structures subjected to large deformation using the rate-of-deformation tensor and the Truesdell rate of Cauchy stress. The Generalized Zigzag Theory (ZIGZAG) presented by Li and Liu [38] is used to account for the laminated composites. Although it can provide excellent linear solutions, ZIGZAG has never been used for finite element formulation. Furthermore, there has never been an attempt to combine the updated Lagrangian approach and any Zigzag Theories for large deformation of laminated composites. The conclusions of this thesis are listed below: 1. A general procedure of deriving the governing equations for structures undergoing 98 99 large deformations (large rotations and large strains) was established. The procedure was an updated Lagrangian approach, which utilized the Truesdell rate of Cauchy stress tensor and the rate-of-deformation tensor. This procedure was general enough that it was not limited to any specific type of displacement fields. . A set of incremental displacement components was successfully derived from the Gen- eralized Zigzag Theory (ZIGZAG). Although the theory described the displacements of laminated composites layer by layer, the total number of unknown variables re mained constant, i.e. independent of layer number. . Based on the ZIGZAG, a four-node shell element (designated as U101) was formulated. It had seven degrees of freedom at each node. The kinky in-plane displacements through the laminate thickness could be naturally displayed. All stress calculations were based on constitutive equations. The transverse shear stresses were continuous through the laminate thickness; no shear correction factor was necessary, and no recovering process through equilibrium equations was required. . An innovative method of constructing the Jacobian matrix was demonstrated. The Zigzag J acobian provided a consistent way of describing both kinematics and geometry within each element. . Extensive benchmark problems were presented by using the U101 element. The pur- pose of the studies was to investigate the performance and to gauge the accuracy of the element. The element gave excellent convergency and accuracy in problems that dealt with isotropic and laminated composites under linear and nonlinear deforma- tions. In addition, the element offered reasonable agreement with some experimental data without any shear locking problem. Although it did not pass a bending patch 100 test, there was no sign of convergency problems. 7.2 Recommendations To bring the U101 element presented herein to more complete fruition, the following topics are recommended for future studies: P1 . Lower continuity requirement for transverse displacement Aw: Instead of using the C 1 Hermite cubic shape function, one may want to use only the C0 continuity shape function for Aw. The use of shape functions of the same order for all degrees of freedom seems to be more consistent in logic. Damage models for delamination analysis: Delamination is a primary damage mode in laminated composites. Implementation of a damage model capable of identifying delamination is critically important to laminated composite analysis. Buckling and post-buckling analyses: Improving the buckling resistance is a challenge task in designing shell-type composite structures. A standard setof benchmark prob- lems should be given to evaluate the capability of the U101 element in buckling and post-buckling analysis. The inclusion of a damage model into the U101 element is a prerequisite in this study. A selectively reduced integration scheme: The U101 element uses a four-point Gauss- Legendre integration scheme. If the number of integration points can be reduced while the degree cf accuracy is maintained, the computation time can be greatly shortened. Triangular element: For a composite structure with irregular boundaries, triangular elements are better than rectangular elements in describing the boundaries. 101 6. Dynamic analysis: Impact response is an important concern in laminated composite analysis. Dynamic analysis is necessary when inertia efl'ect has to be considered. 7. Independent finite element program: The U101 element is currently implemented as a subroutine in ABAQUS/Standard. Every time the program is to be executed, the FORTRAN subroutine needs to be compiled and linked to the main code. This process is very time consuming. Besides, the utilization of data storage space is not optimized. If the element is to become more versatile, a stand-alone program would be desired. APPENDICES APPENDIX A CON SISTENT LINEARIZATION The meaning of consistent linearization here is to use a systematic process to obtain linear approximation of a set of nonlinear equations within a small neighborhood of some known“— approximate solution to the nonlinear problem. Let us consider the motion of a material(in this case, it is also equal to the mesh motion) is given by x]? = X,- +uf, (A.1) and the new motion is 1:?“ = z¥+Aui 1, 11:”? = X,“ +111: ‘1' A211,] l H ‘ (A2) in which the superscript v is interpreted as the increment ”counter”. Now consider an abstract nonlinear function, .7, of components of a vector field If“, that satisfies the following nonlinear equation 102 103 0+1 1' Assuming that f(.r ) is sufficiently smooth (i.e. differentiable) in the neighborhood of a given (known) ” point”, 932’, it maybe expanded in a Taylor series about If to obtain \/ 1'1 -, ..g 1 (,1. \ J:(.'r'-’+1) = fa?) + 527(3? + cAu,) + RU?) (AA) ' s=0 in which the sum of first two terms of the right hand side of Eq.(A.4) is the linearized part of .7: and it is a measure of the rate of change of .77 in the direction of Au, at the point :63. And R is the remaining higher order term, which by Taylor‘s theorem vanishes faster than the linearized part as Au,- —> 0. c is a scale parameter. Denoting the second right-hand-side term by as d 1' ‘ ()1 cm: 45514248214) 0 A (4.5) e: By neglecting the R term and employing Eqs.(A.3), (AA) and (A.5), we obtain the following linear approximation (locally) to the original nonlinear problem - Eq.(A.3) -£[fl.:=f(x¥) ' " . (A.6) where the only unknown is Au,-, all other quantities being evaluated at the known reference point 23?. Iterative use of Eqs.(A.6) and (A.2) may be identified with the Newton-Raphson solution algorithm, where Eq.(A.6) is used to computer the iterative change (Au,) and Eq.(A.2) updates the solution between increments. To simplify subsequent writing, the subscript 2:}? is dropped with an understanding that Elf] = Llflxf (A-7) The following few examples of [[7] will be helpful for the later derivation. 1. the deformation gradient, F},- = 0x,/6Xj, _ BAu, d 6 2:? + (Am) £[Ej] = { ( } (“-0 — an 2% 6X,- 104 2. the Jacobian, J = det (F11), _ d 6(1):) + cAu,))} £[Jl _ d. {det( an (:0 = PAM-J- pass,- (11.9) 61" 3. inverse of the deformation gradient, F51, FF-l = 1 £[F]F‘1 + F£[F'1] = 0 £[F‘1] = 4141117111-1 (A.10) - Therefore, the components of £[F‘1] are: aAu -1 _ -1 k —1 APPENDIX B SOME SPECIAL MATRICES FOR FINITE ELEMENT FORMULATION [X U0] is a 1 by 28 vector and its non-zero components are listed as follows, for 1' = 1, 2, 3, 4: XU0(1,7(1' — 1) +1) = N,- XU0(1, 7(1' — 1) + 3) = 53(61111/611) XU0(1,7(1' — 1) +4) = 5M, XU0(1,7(:' — 1) + 6) = S; (ails/as) XU0(1,7(1' — 1) + 7) = s; (airs/ax) [X U 1] is a 1 by 28 vector and its non-zero components are listed as follows, for i = 1, 2, 3, 4: XU,(1,7(1' — 1) +3) 35(37'11'1/53) XU,(1,7(1' — 1) +4) = R’fN, 105 106 XU1(1, 7(1‘ - 1) + 6) = R}; (min/62:) XU,(1,7(1' — 1) + 7) R}; (sing/as) [X U2] is a 1 by 28 vector and its non-zero components are listed as follows, for i = 1, 2, 3, 4: XU2(1,7(1' — 1) +3) = A2(6’H1‘1/0:r) XU2(1,7(1' — 1) + 4) = Ami,- XU2(1, 7(1' - 1) + 6) A2 (min/8.1:) XU2(117(i-1)+ 7) = 42(5711'3/393) [X U3] is a 1 by 28 vector and its non-zero components are listed as follows, for 1' = 1,2,3, 4: XU3(1,7(1° — 1) +3) = 32(3711'1/02') XU3(1,7(3-1)+4) = BlN,’ XU3(1, 7(1' — 1) + 6) = 82(31'112/03) XU3(1,7(1° - 1) +7) 32 (3716/ 53 l [X V0] is a 1 by 28 vector and its non-zero components are listed as follows, for i = 1,2, 3,4: XV0(1,7(1' - 1) +2) N,- X1/b(1,7(1’—1)+ 3) = P2" (min/09) XV0(1,7(1' — 1) + 5) PfN, XVo(1, 7(1' — 1) + 6) = P2’°(611.-2/69) xvo(1, 7(1' — 1) + 7) = P41681069) [X V1] is a 1 by 28 vector and its non-zero components are listed as follows, for i = 1, 2, 3,4: XV,(1, 7(1' — 1) + 3) = 05(61111/011) 107 XV,(1,7(i—1)+5) = 01w, XV,(1,7(1' — 1) + 6) = 05(6H12/6y) ”111.711 - 1) + 7) = 051618.061) [X V2] is a 1 by 28 vector and its non-zero components are listed as follows, for i = 1, 2, 3, 4: XV2(117(i-1)+3) = Cg(6’H,~1/6y) XV2(1, 7(1' - 1) + 5) = GIN,- XV2(117(i-1)+6) C'2 (min/31.1) XV2(117(1'- 1)+ 7) 02(671113/39) [X V3] is a 1 by 28 vector and its non-zero components are listed as follows, for i = 1, 2, 3, 4: XV:1(117(i-1)+ 3) = 02(37‘11/39) XV3(1, 7(11— 1) + 5) = 0,111,- XV3(117(i-1)+6) = 02(37'11'2/39) XV3(1, 7(1' - 1) + 7) D2 (Wis/69) [X W0] is a 1 by 28 vector and its non-zero components are listed as follows, for 1' = 1, 2, 3, 4: XWo(1,7(i — 1) + 3) = 7111 XW0(1,7(i — 1) + 6) = 7112 XWo(l,7(i - 1) + 7) = 7113 [MA] is a 5 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: MA(1,7(1' — 1) + 1) = 6N,/62: MA(1, 7(1' — 1) + 3) S; (02H1'1/0x2) 108 MA(1,7(1‘ - 1) +4) -_— Sf (BM/(9.7:) MA(1,7(1' — 1) + 6) = 55(627112/611?) MA(1,7(1‘ — 1) + 7) = 5; (ems/as?) MA(2,7(1'—1)+2) = (am/81) MA(2,7(1‘—1)+3) = P2"(627-l,-1/6y2) MA(2,7(i—1)+5) = P,*(6N,-/ay) MA(2,7(1'-1)+6) = P2" (627112/6‘92) MA(2,7(i—l)+7) = P§(a27i,-3/6y2) MA(3,7(1‘—1)+1) = away MA(3,7(i-1)+2) = 6N,/8:c MA(3,7(1‘ —— 1) +3) = 55(621lu/asay) +P2"(621£,-1/6:l:0y) MA(3,7(1‘—-1)+4) = sf(aN,-/ay) MA(3,7(i—1)+5) = Pf(aN,-/as) MA(3,7(1‘ - 1)+6) = s§(a’7r,-2/axay) +12; (ems/asap) MA(3,7(1‘— 1)+7) = s; (6211,3/62611) +19," (ems/asap) [MB] is a 5 by 28 matrix and its non-zero components are listed as follows, for 1' = 1, 2, 3, 4; MB(4,7(1' — 1) + 3) (02+1)(6Hil/ay) MB(4,7(1' — 1) + 5) = ofN, MB(4,7(1' - 1) + 6) = (0§+1)(8’H12/6y) MB(4, 7(1' — 1) + 7) (05 + 1) (ans/69) MB(5, 7(1' - 1) + 3) (17!; + 1) (dun/ax) MB(5, 7(1' — 1) + 4) = R’fN, 109 MB(5,7(z'—1)+6) = (R’;+1)(afl.—2/as) MB(5,7(z'—1)+7) = (ng+1)(aH,-s/as) [M D1] is a 5 by 28 matrix and its non-zero components are listed as follows, for i = 1,2,3, 4: 11104176“ - 1)+3) R; (027-13'1/632) M041, 7(1’ — 1) + 4) [If (BM/63:) M041, 7(5 — 1) + 6) 12’; (8271,4622) M041, 7(i — 1) + 7) R; (@2711'3/332) M042, 7(i - 1) + 3) 0!,c (6274,4613) M042, 7(i — 1) + 5) 0’: (BM/By) M042, 7(i - 1) + 6) 0!; (ems/6y”) M042, 7(i — 1) + 7) 0!; (ems/6y?) M043, 7(i - 1) + 3) (R2 + 02) (awn/62:63,) M043, 7(1' — 1) + 4) R1 (aNi/By) M043, 7(i — 1) + 5) 01 (am/as) M043, 7(i — 1) + 6) (R2 + 02) (62Hi2/616y) M043, 7(1‘ — 1) + 7) (122 + 02) (0271;3/62631) MD) (4, 7(1' — 1) + 3) 202(6714/617) M044, 7(i — 1) + 5) 2C1N, M044, 7(i — 1) + 6) 202 («ms/61y) M044, 7(2’ — 1) + 7) 202 (mas/By) M045, 7(i — 1) + 3) 242 (Wu/31:) M045, 7(1' — 1) + 4) 244v, _=1fn-h . In. 110 MD1(5,7(3' — 1) + 6) = 242 (malts/as) MD45, 7(t' - 1) + 7) 242(aw,3/as) [M02] is a 5 by 28 matrix and its non-zero components are listed as follows, for 1' = 1, 2, 3, 4; M041,7(i - 1) + 3) = A2(82H,1/6x2) M041,7(i — 1) + 4) = ,41 (BM/62:) M041,7(i—1)+6) = 42(azuss/as2) MD2(1,7(i-1)+7) = A2(62H.-3/6a:2) M042,7(i-1)+3) = 02(627ii1/6172) M042,7(i—1)+5) = C1 (BM/6y) M042,7(i—1)+6) = 02(6221,2/6~y2) M042,7(i — 1)+7) = Cg (wuss/6y?) M043, 7(1— 1)+3) = (A2+C2)(62’H.-1/8$6y) M043,7(t'-1)+4) = A1(6N.-/6y) MD2(3,7(i-1)+5) = 01(aN,/as) M043, 7(i— 1) +6) = (A2+Cg)(627-(.-2/82:6y) M043, 7(i— 1)+7) = (A2+Cs)(62u,s/asay) M044,7(i—1)+3) = 3D2(6’H.-1/6y) M044,7(i—1)+5) = 304v,- MD2(4,7(i—1)+6) = 302(37112/314) MD2(4,7(i-1)+7) = muons/617) M045,7(i—1)+3) = 332(6114/as) MDQ(5, 7(1. - 1) +4) = 381Ni 111 M045, 7(t' - 1) + 6) 30.2 (Wig/Bx) M045, 7(i — 1) + 7) 332 (miss/as) [M D3] is a 5 by 28 matrix and its non-zero components are listed as follows, for i = 1,2, 3,4: MDg(1,7(z'—l)+3) = 02(62714/as2) MD3(1,7(i — 1) + 4) = 81(6Ns/62) M041,7(i—1)+6) = 32(62uss/as2) MD3(1,7(i—1)+7) = 32(827-143/822) M03(2,7(t'—1)+3) = 02(621tsi/ay2) MD3(2,7(i—1)+5) = D1 (BM/0y) M042,7(t'—1)+6) = D2(027{.-2/6y2) MD3(2,7(z‘—1)+7) = D2(62H.-3/6y2) MD3(3,7(i—1)+3) = (024-02) (ems/assay) M043,7(t'—1)+4) = 31(0Ni/6y) MD3(3,7(z'—1)+5) = 046117462) MD3(3,7(i—1)+6) = (32+02)(32u,~2/as6y) M043, 7(2’ — 1) + 7) = (132 + 02) (yum/3363’) [0%] is a 3 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: GUO(1,7(i — 1) +1) = aNi/az' GU41, 7(i — 1) + 3) = s; ((927-14/832) GU41, 7(i — 1) + 4) sf (am/as) GU0(1,7(i-1)+6) = 55(62Hi2/622) GU41,7(i - 1) + 7) GU42,7(i - 1) + 1) GU42, 7(i — 1) + 3) GU0(2,7(i - 1) + 4) GU42, 7(i — 1) + 6) GU42, 7(1' — 1) + 7) GU0(3,7(i — 1) + 3) GU43,7(i — 1) + 4) GU43,7(i — 1) + 6) GU43, 7(i — 1) + 7) GU41,7(t' — 1) + 3) GU41,7(i - 1) +4) GU41,7(i — 1) + 6) GU41,7(i — 1) +7) GU42, 7(i — 1) + 3) GU1(2,7(1' — 1) +4) GU1(2,7(3' — 1) + 6) GU1(2,7(3' -‘1)+ 7) GU1 (3, 70‘ — 1) + 3) GU1(3,7(i - 1) +4) GU1(3,7(3' — 1) +6) 112 5; (6271mm?) away 55 (621-14/62-611) Si (0Ni/0y) s; (ems/assay) 5; (0214.3 /6$6‘y) R: (ans/ax) R’fN, Rt (miss/ax) R; (37113/ 327) [GU 1] is a 3 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: R”; (ems/as?) R’i (aNs/az) RI; (3271,4622) 0'; (6211.3 ms?) R; (627-14 /6x8y) 6’: (aNs/ay) 3'; (627-10 /6$6y) 0’; (6211.3 /82:8y) 2,42 (631,-, /as) 2A1N, 2A2 (3740/51?) GU1(3,7(i - 1) + 7) 113 2A2 (37113/ 61-) [GUg] is a 3 by 28 matrix and its non-zero components are listed as follows, for z' = 1, 2, 3, 4: GU41,7(i — 1) + 3) GU2(1,7(z’ — 1) +4) GU41,7(i — 1) +6) GU41,7(i — 1) + 7) GU42,7(i — 1) + 3) GU42,7(i — 1) + 4) GU42,7(i — 1) + 6) GU42,7(i — 1) + 7) GU43,7(i - 1) + 3) GU43, 7(1‘ — 1) + 4) GU43, 7(1‘ - 1) + 6) GU43, 7(1‘ — 1) + 7) A2 (62%,) we?) A1 (am/as) 42 (8211.2 /31:2) 42 (6226-468) A2 (ems/assay) At (BM/5‘11) 42 (62%2/31'817) A2 (62%3/82631) 332 (87-1,) /0:c) 381M- 332 (8%2/02) 3B2 (3710/ 31?) [GU3] is a 3 by 28 matrix and its non-zero components are listed as follows, for 3' = 1, 2, 3, 4: GU3(1,7(1' — 1) +3) GU41,7(i — 1) +4) GU41,7(i — 1) +6) GU3(1,7(i - 1) +7) GU42,7(i - 1) + 3) GU3(2,7(3' — 1) + 4) 32 (6271,1/6222) 131 (6N,- /6:1:) 32 (ems/ax?) 32 (ems/ax?) 32 (emu/assay) Bl (aNi/By) GU3(2,7(i — l) + 6) GU3(2,7(3' — 1) + 7) G'Vo(1,7(i — 1) + 2) GVo(1,7(i — 1) + 3) G‘Vo(1,7(i — 1) + 5) GVo(1,7(i — 1) + 6) GV41,7(i — 1) + 7) GV42,7(i — 1) + 2) GVo(2,7(i — 1) + 3) GV42,7(i — 1) + 5) G'Vo(2,7(i — 1) + 6) GV42, 7(i — 1) + 7) GV43,7(i - 1) + 3) GVo(3,7(z' — 1) + 5) GVo(3,7(i — 1) + 6) GV43, 7(t’ — 1) + 7) GV41,7(i — 1) +3) GV41,7(i — 1) +5) 114 B2 (32%2/3’3511) 32 (32H13/326y) [GVO] is a 3 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: aim/as- P2" (6211,1/61631) Pf (6N,- /0:7:) P; (6221.6 /3$(9y) P2" (6211.3 /8:7:6y) 3Ni/5‘y P2" (emu/6y?) P1" (3Ni/0y) P2" (6271.2 My?) P2" (ems/as?) 0% (3741/ 6y) O’fN, 02‘ (37112/ By) 02 (37(13/ (931) [0V1] is a 3 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: 05 (027ii1/6zr2) 0’: (aNs/ax) GV,(1,7(1’ — 1) +6) GV,(1,7(2’ — 1) + 7) GV,(2,7(i — 1) + 3) GV,(2,7(3’ — 1) + 5) GV,(2,7(i — 1) + 6) GV,(2,7(1' — 1) + 7) GV,(3,7(2‘ - 1) + 3) GV43, 7(1' — 1) + 5) GV,(3,7(i — 1) + 6) GV43, 7(i — 1) + 7) GVg(1,7(i — 1) + 3) GVgU,7(i - 1) + 5) GV41,7(i — 1) +6) GV41,7(i — 1) + 7) GV42,7(t‘ - 1) + 3) GV42,7(i — 1) + 5) GV42,7(i — 1) + 6) GV42,7(i - 1) + 7) GV2(3,7(1’ — 1) + 3) GV43,7(i — 1) + 5) GV2(3, 7(3' -— 1) + 6) 115 202 (Wu/3y) 2011)],- 202 (3710/53!) 202 (37113/ all) [GV2] is a 3 by 28 matrix and its non-zero components are listed as follows, for 3' = 1, 2, 3, 4: 02 (827-!4/83637) G, (BM/82:) G2 (ems/assay) G, (6221,, /asay) 02 (ems/as?) C'i (BM/011) G2 (62%,, /6y2) G2 (62%,, /8y2) 302 (6H,, /0y) 30,1v, 302 (min/3y) 116 GV2(3.7(i - 1) + 7) = 302(37'113/311) [0V3] is a 3 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: GV3(1,7(i—1)+3) = 02(627-{4/62631) GV41, 7(i — 1) + 5) = 0, (BM/3x) GV3(1,7(i—1)+6) = 02 (6211,2/az6y) GV3(1,7(i—1)+7) = 02(62H,3/az6y) GV42,7(i-1)+3) = 02(0’713/6172) GV42,7(i—1)+5) = 01(6Ni/6y) GV3(2,7(i—1)+6) = 02(02Hi2/6y2) GV42, 7(t' — 1) + 7) = D2 (Ema/6y?) [GWo] is a 3 by 28 matrix and its non-zero components are listed as follows, for i = 1, 2, 3, 4: GWo(1,7(i—1)+3) = (ms/as GWo(1,7(z’-1)+6) = ems/6s GW41,7(i—1)+7) = ans/6s GW42,7(i—1)+3) = ans/6y GW42,7(i—1)+6) = Mia/6y GW42, 7(3' — 1) + 7) Wis/3y APPENDIX C TRUESDELL IN CREMENTAL CONSTITUTIVE EQUATION Combining Eqs.(2.16) and (2.23), we obtain Vat V ij = aij + OEAUk’k - UfiAUJ',‘ + Auuafj + ngazj - Okakj (C.1) Using Eqs.(3.28) and (4.8), Eq.(C.1) can be expressed as ijtlAulkJ) = Common + 053“” “ “3M“ — aglAmJ + Au[t',k]0:j + akauW] = CijltlAu(k,l) + Oi’jAuch - 011416.: ‘ ”ilAuiJ l v ,, 1 s v + 5 (aikAuUM + 0,),5lkAuw1) + 5 (ajkAuliJcl + ajk‘sk‘AuliJl) = CijklAu(lt,l) + 031130”: - ”HAUL! ‘ 071mm 1 + 2 (03416.11 + UilAulul + o?,,Au,,-,s + “ikAulukll = CijklAuUcJ) + Gig-Au“, - aflAuJ-J — aflAuu 117 118 1 + Z(G,‘-"Auj,z -— afiAulJ + aylAuu —— aflAuM "Av—”A-VA--'-’A + 03k “M at]: “k.3+0)k “th 01k Wu) 1 = CijklAu(k,1) + 03AM”: + §(-03A“j,t " UilAujJ l - aykAuch — ijAu,,k) + 2(03A‘UjJ — GfiAuu +o”A-—”A -+’-’A«—”A- jl “3.1 Ojl "1,: 03’: "1.1: aik We.) + aykAu,,k—G}’kAuk,,-) 1 = CijlclAu(k,l) + 50351:! (Aunt 'l' Aulsk) 1 — :1- (Gz’kAuJ-J, + 03AM; + akaui-s + UflAws') 1 _ Z (callus,- + 034319.12“ 03km“ + ”ilAuiv‘) = CijklAu(k,l) + Oi’j‘lflAWkJ) l — Z (03“st + 0361']; + 03"];63'1 + 03316“) (Auksl + Aulsk) = CijklAu(k,l) + Ui’j‘WAWkJ) 1 _ 5 (one), + 036,-), + ogkos + 7,396.4) A110,,” ((3.2) Therefore, we obtain Cijkl = Cijkl + 035“ l — 2 (034%; + 036,-), + 0&6.) + 03163:) (0'3) APPENDIX D ABAQUS/Standard. CONVERGENCY CRITERIA D.1 The Solution of Nonlinear Problems While the commercial program ABAQUS/Standard is still evolving, the objective of this chapter is to document the convergency criteria used in the current thesis. The description of convergency criteria is reproduced from ABAQUS/ Standard User’s Manual [2] Sec.8.2.1. ABAQUS/Standard (ABAQUS in subsequent writing) uses a direct, Gauss elimination method to solve a system of simultaneous linear algebraic equations. Usually, many sets of simultaneous linear algebraic equations must be solved to obtain nonlinear solutions. The nonlinear load-displacement curve for a structure is shown in Figure D.1. The objective of the analysis is to determine this response. In a nonlinear analysis the solution cannot be calculated by solving a single system of nonlinear equations, as would be done in a linear problem. Instead, the solution is found by specifying the loading as a function of time and incrementing time to obtain the nonlinear response. Therefore, ABAQUS breaks 119 120 the simulation into a number of load increments and finds the approximate equilibrium configuration at the end of each load increment. Using the Newton-Raphson method, it often takes ABAQUS several iterations to determine an acceptable solution to each load increment. D.2 Steps, Increments, and Iterations 1. The time history for a simulation consists of one or more steps. The user defines the steps, which generally consist of an analysis procedure option,_loading options, and output request options. Different loads, boundary conditions, analysis procedure options, and output requests can be used in each step. 2. An increment is part of a step. In nonlinear analyses each step is broken into in- crements so that the nonlinear solution path can be followed. The user suggests the size of the first increment, and ABAQUS automatically choose the size of the sub- sequent increments. At the end of each increment the structure is in (approximate) equilibrium and results are available for writing to the restart, data, or results files. 3. An iteration is an attempt at finding an equilibrium solution in an increment. If the model is not in equilibrium at the end of the iteration, ABAQUS tries another iteration. With every iteration the solution that ABAQUS obtains should be closer to equilibrium; however, sometimes the iteration process may diverge - subsequent iterations may move away form the equilibrium state. In that case ABAQUS may terminate the iteration process and attempt to find a solution with a smaller increment size. m1 D.3 Convergency Consider the internal (nodal) forces, I, and the external force, P, acting on a body. The internal load acting on a node are caused by the stresses in the elements that are attached to that node. For the body in equilibrium, the net force acting at every node must be zero. Therefore, the basic statement of equilibrium is that the internal force I, and the external force, P, must balance each other: P—I=0 (on The nonlinear response of a structure to a small load increment, AP, is shown in Figure D.1. ABAQUS uses the structure’s tangent stiffness, K0, which is based on its configuration at no, and AP to calculate a displacement correction, ufi, for the structure. Using ufi, the structure’s configuration is updated to no. ABAQUS then calculates the structure’s internal force, Ia, in this updated configuration. The difference between the total applied load, P, and In can now be calculated as \ m=P—n (am where R, is the force residual for the iteration. If R, is zero at every degree of freedom in the model, point a in Figure D.1 would lie on the load deflection curve and the structure would be in equilibrium. In a nonlinear problem Ra will never be exactly zero, so ABAQUS compares it to a tolerance value. If R0 is less than this force residual tolerance at all nodes, ABAQUS accepts the solution as being in equilibrium. By default, this tolerance value is set to 0.5% of an average force in the structure, averaged over time. ABAQUS automatically calculates this spatially and time-averaged force throughout the simulation. If R0 is less than the current tolerance value, P and 1,, are considered to be in equilibrium 122 and 11,, is a valid equilibrium configuration for the structure under applied load. However, before ABAQUS accepts the solution, it also checks that the largest displacement correction, 112, is small relative to the total incremental displacement, Ana = 11,, — no. If u: is greater , than a fraction (1% by default) of the incremental displacement, ABAQUS performances another iteration. Both convergence checks must be satisfied before a solution is said to have converged for that load increment. If the solution from an iteration is not converged, ABAQUS performs another iteration to try to bring the internal and external forces into balance. First, ABAQUS forms the new stiffness, K0, for the structure based on the updated configuration, ua. This stiffness, together with the residual Ra, determines another displacement correction, ufi, that bring the system closer to equilibrium (point b in Figure D.1). ABAQUS calculates a new force residual, Rb, using internal forces from the structure’s new configuration, 13),. Again, the largest force residual at any degree of freedom, R4,, is compared against the force residual tolerance, and the displacement correction for the second iteration, 113, is compared to the increment of displacement, Aub. If necessary, ABAQUS performs further iteration. D.4 Automatic Incrementation Control. By default, ABAQUS automatically adjusts the size of the time increments to solve non- linear problems efficiently. The user needs to suggest only the size of the first increment in each step of the simulation, after which ABAQUS automatically adjusts the size of the increments. If the user does not provide a suggested initial increment size, ABAQUS will attempt to apply all of the load defined in the step in a single increment. For highly non- linear problems, ABAQUS will have to reduce the increment size repeatedly to obtain a solution. 123 The number of iterations needed to find a converged solution for a time increment will vary depending on the degree of nonlinearity in the system. With the default incrementation control, the procedure works as follows. If the solution has not converged within 16 iterations or it the solution appears to diverge, ABAQUS abandons the increment and starts again with the increment size set to 25% of its previous value. It then attempts to find a converged solution with this smaller load increment. If the solution still fails to converge, ABAQUS reduces the increment size again. This process is continued until a solution is found. If the time increment becomes smaller than the minimum defined by the user or more than five attempts are needed, ABAQUS stops the analysis. If the increment converges in fewer that than iterations, this indicates that the solution is being found fairly easily. Therefore, ABAQUS automatically increases the increment size by 50% if two consecutive increments require fewer that 5 iterations to obtain a converged solution. 124 Load D Displacement AUb .1. Figure D.1 First and second iteration of an increment BIBLIOGRAPHY BIBLIOGRAPHY [1] ABAQUS/Standard Example Problems Manual (1996), Volume I & [1, version 5.6, Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, RI. [2] ABAQUS/Standard User’s Manual (1996), Volume I, II & 11, version 5.6, Hibbitt, Karlson & Sorensen, Inc., Pawtucket, RI. [3] ABAQUS Theory Manual (1995), version 5.5, Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, RI. [4] ABAQUS/Standard Verification Manual (1995), version 5.5, Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, RI. [5] Altiero, N. (1987), Finite Element Method, Course Notes, Michigan State University, East Lansing, Michigan. [6] Barbero, E.J. and Reddy, J .N. (1990), ”An Accurate Determination of Stresses in Thick Laminates Using a Generalized Plate Theory,” International Journal for Nu- merical Methods in Engineering, 29, 1-14. [7] Barbero, E.J. and Reddy, J.N. (1990), ”Nonlinear Analysis of Composite Laminates Using a Generalized Laminated Plate Theory,” AIAA Journal, 28, 1987-1994. \[8] Bathe, K. J.( (1996), Finite Element Procedures, Prentice-Hall, Englewood Cliffs, New Jersey. [9] Bazeley, G.P., Cheung, Y.K., Irons, B.M., and Zienkiewicz, O.C. (1966), ”Triangular Elements in Bending - Conforming and Non-conforming Solutions,” Proc. Ist Conf. Matrix in Structural Mechanics, AFFDLTR-CC-SO, Air Force Inst. Tech, Wright- Patterson AF Base, Ohio, 547-576. [10] Belytschko, T. and Hughes, T.J.R. ed. (1983), Computational Methods in Mechanics - Volume 1: Computational Methods in Transient Analysis, Chapter 1, North-Holland Publisher, New York, New York. [1 1] Belytschko, T., Lin, J .I. and Tsay, C.-S. (1984), ”Explicit Algorithms for the Nonlinear Dynamics of Shells,” Computer Methods in Applied Mechanics and Engineering, 42, 225-251. 125 126 [12] Chia, CY. (1972), ”Large Deflection of Rectangular Orthotropic Plates,” Journal of Engineering Mechanics Division, Proceeding, ASCE, 98, 1285-1298. [13] Chia. CY. (1980), Nonlinear Analysis of Plates, McGraw-Hill, Inc., New York, New York. \/ [14] Cook, RD. (1989), Concepts and Application of Finite Element Analysis, 3rd ed., Wiley, New York. [15] Crisfield MA. (1991), Non-linear Finite Element Analysis of Solids and Structures - Volume 1: Essentials, John Wiley 85 Sons, New York. [16] Dienes, J.K. (1979), ”On the Analysis of Rotation and Stress Rate in Deforming Bodies,” Acta Mechanica, 32, 217-232. [17] Di Sciuva, M. (1992), ”Multilayered Anisotropic Plate Models with Continuous Inter- laminar Stress,” Composite Structures, 22, 149-167. [18] Flanagan, DP. and Taylor, L.M. (1987), ”An Accurate Numerical Algorithm for 'Stress Integration with Finite Rotations,” Computer Methods in Applied Mechanics and Engineering, 62, 305-320. N/ [19] thg, Y.C. (1965), Foundations of Solid Mechanics, Prentice-Hall, Englewood Cliffs, New Jersey. V! [20] Halleux, J .P. and Donea, J. (1985), ”A Discussion of Cauchy Stress Formulations for Large Strain Analysis,” Finite Element Methods for Nonlinear Problems: Europe-US Symposium, Bergan, Bathe and Wunderlich ed., Trondheim, Norway, 61-74. [21] Hinton, E. and Owen, D.R.J. (1979), An Introduction to Finite Element Computa- tions, Pineridge Press Limited, Swansea, UK. \ [22] Hughes, T.J.R. (1987), The Finite Element Method - Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cliffs, New Jersey. [23] Hughes, T.J.R. and Belytschko, T. (1995), Nonlinear Finite Element Analysis, notes of a short course, Palo Alto, California. [24] Hughes, T.J.R. and Liu, W.K. (1981), ”Nonlinear Finite Element Analysis of Shells: Part I. Three-dimensional Shells,” Computer Methods in Applied Mechanics and En- gineering, 26, 331-362. [25] Irons, BM. (1974), ”The Patch Test for Engineers,” Proc. Finite Symp. Atlas Com- puter Lab., Chilton, Didcot, England, 171-192. 127 V/ [26] Kapania, R.K. and Raciti, s. (1989). ”Nonlinear Vibration of Unsymmetrically Lam- inated Beams,” /em AIAA Journal, 27, 201-210. [27] Kapania, R.K. and Raciti, S. (1989), ”Recent Advances in Analysis of Laminated Beams and Plates, Part I: Shear Effects and Buckling,” AIAA Journal. v27n7, 923- 934. [28] Kardestuncer, H. (1989), Editor in Chief, Finite Element Handbook, Part 2, Chapter 6, McGraw-Hill, New York. [29] Knight, N.F. Jr. and Starnes, J.H. Jr. (1985), ”Postbuckling Behavior of Axially Compressed Graphite-Epoxy Cylindrical Panels with Circular Holes,” Transactions of the ASME, 107, 394-402. [30] Kweon, J .H. and Hong, CS. (1993), ”Postbuckling Analysis of Composite Laminated Cylindrical Panels Under Axial Compression,” AIAA Journal, v31n8, 1535-1537. [31] Kweon, J.H., Hong, GS. and Lee, LC. (1995), ”Postbuckling Compressive Strength of Graphite/Epoxy Laminated Cylindrical Panels Loaded in Compression,” AIAA Journal, v33n2, 217-222. [32] Lee, C.Y. (1991), A Study of the Interlaminar Stress Continuity Theories for Compos- ite Analysis, Ph.D. Dissertation, Dept. of Materials Science and Mechanics, Michigan State University, East Lansing, Michigan. [33] Lee, C.Y. and Liu, D. ( 1992), ”An Interlaminar Stress Continuity Theory for Lami- nated Composite Analysis,” Computers and Structures, v42n1, 69-78. [34] Lee, C.Y. and Liu, D. (1993), ”Nonlinear Analysis of Composite Plates Using Inter- laminar Shear Stress Continuity Theory,” Composite Engineering, v3n2, 151-168. [35] Lee, L.J., Huang, K.Y., Hwang, DC. and Pai, CK. (1993), ”Geometrical Nonlinear Analysis of Composite Cylindrical Shell Panels Subjected to Impact,” Finite Elements in Analysis and Design, 15, 135-149. ._ [36] Lee, Y.K. (1992), ”Shear Banding and Material Instability in Finite Deformed Com- pressible Materials: An Example in Thin Sheets,” International Journal of Plasticity, v8, 543-565. [37] Li, X. (1994), Investigation on Laminated Plate Theories, Ph.D. Dissertation, Dept. of Materials Science and Mechanics, Michigan State University, East Lansing, Michigan. [38] Li, X. and Liu, D. (1995), ”Zigzag Theory for Composite Laminates.” AIAA Journal, v33n6, 1163-1165. [39] Li, X. and Liu, D. (1997), ”Generalized Laminated Theories Based on Double Super- position Hypothesis,” International Journal for Numerical Methods in Engineering, v40n8, 1197-1212. 128 [40] Liao, CL and Reddy, J.N. ( 1989), ”Continuum-Based Stifl'ened Composite Shell Element for Geometrically Nonlinear Analysis,” AIAA Journal, v27n1, 95-101. [41] Liu, D. (1988), ”Impact-induced Delamination - A View of Bending Stiffness Mis- matching,” Journal of Composite Materials, 22, 674-691. [42] Liu, D. and Li, X. (1996), ”An Overall View of Laminate Theories Based on Displace- ment Hypothesis,” Journal of Composite Materials, 30, 1539—1561. [43] Liu, WK. (1988), Advanced Finite Element Methods, Course Notes, Northwestern University, Evanston, Illinois. [44] Lo, K.H., Christensen, R.M., and Wu, EM. (1977), ”A High-order Theory of Plate Deformation - Part 1: Homogeneous Plates,” Journal of Applied Mechanics, v44, 663-668. - [45] Lo, K.H., Christensen, B.M., and Wu, EM. (1978), ”Stress Solution Determination for High Order Plate Theory,” International Journal of Solids and Structures, v14. 655-662. ' [46] La, X. and Liu, D. (1992), ”An Interlaminar Shear Stress Continuity Theory for Both Thin and Thick Composite Laminates,” Journal of Applied Mechanics. 59, 502-509. [47] MacNeal, R.H. and Harder, BL. (1985), ”A proposed Standard Set of Problems to Test Finite Element Accuracy,” Finite Elements in Analysis and Design, 1, 3-20. [48] Malvern, LE. (1969), Introduction to the Mechanics of a Continuous Medium, Pren- \/ tice Hall, Englewood Cliffs, New Jersey. [49] Mase, GE. and Mase, GT. (1992), Continuum Mechanics for Engineering, CRC Press, Boca Raton, Florida. [50] Mattiasson, K., Bengtsson, A. and Samuelsson, A. (1985), ”On the Accuracy and Efficiency of Numerical Algorithms for Geometrically Nonlinear Structural Analysis,” Finite Element Methods for Nonlinear Problems: Europe-US Symposium, Bergan, Bathe and Wunderlich ed., '11'0ndheim, Norway, 3-23. [51] Mindlin, RD. (1951), ”Influence of Rotatory Inertia and Shear on Flexural Motions of Isotropic Elastic Plates,” Journal of Applied Mechanics, v18, 31-38. [52] Osias, J.H. and Swedlow, J.L. (1973), ”Finite Elasto-Plastic Deformation-I: Theory and Numerical Examples,” International Journal of Solids Structures, v10, 321-339. [53] Pagano, N .J . (1969), ”Exact Solutions for Composite Laminates in Cylindrical Bend- ing,” Journal of Composite Materials, v3, 398-411. [54] Pagano, NJ. (1970), ”Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates,” Journal of Composite Materials, v3, 20-35. 129 [55] Palazotto, AN. and Dennis, S.T. ( 1992), Nonlinear Analysis of Shell Structures, AIAA Education Series. [56] Pipes, R.B. and Pagano, N .J . (1974), ”Interlaminar Stress in Composite Laminates - An Approximate Elasticity Solution,” Journal of Applied Mechanics, 41. 668-672. [57] Prager, W. (1961), ”An Elementary Discussion of Definitions of Stress Rate,” Quart. A ppl. Math., v18n4, 403-407. [58] Prager, W. (1962), ”On Higher Rates of Stress and Deformation,” Mech. Phys. Solids, 10, 133-138. [59] Reddy, J.N. (1982), ”Analysis of Layered Composite Plates Accounting for Large Deflections and Transverse Shear Strains,” Recent Advances in Non-linear Computa- tional Mechanics, Chapter 6, 155-203, Pineridge Press Limited, Swansea, UK. [60] Reddy, J.N. (1984), ”A Refined Nonlinear Theory of Plates with 'Ii'ansverse Shear Deformation,” International Journal of Solids and Structures, 20, 881-896. [61] Reddy, J .N. (1984), ”A Simple Higher-order Theory for Laminated Composite Plates,” Journal of Applied Mechanics, 51, 745-752. [62] Reddy, J.N. (1984), Energy and Variational Methods in Applied Mechanics, Wiley, New York. [63] Reddy, J .N. (1990), ”A General Non-linear Third-order Theory of Plates with Mod- erate Thickness,” International Journal of Non-linear Mechanics, 25, 677-686. [64] Reddy, J .N . and Chandrashekhara, K. (1985), ”Nonlinear Analysis of Laminated Shells Including 'Il'ansverse Shear Strains,” AIAA Journal, v23n3, 440-441. [65] Reissner, E. (1945), ”The effect of 'Itansverse Shear Defamation on the Bending of Elastic Plates,” Journal of Applied Mechanics, v12, A69-A77. [66] Singh, G. and Sadasiva Rao, Y.V.K. (1987), ”Nonlinear Vibration of thick Composite Plates,” Journal of Sound and Vibration, 115, 367-371. [67] Srinivas, S. and Rao, A.K. (1970), ”Bending, Vibration, and Buckling of Simply Sup- ported Thick Orthotropic Rectangular Plates and Laminates,” International Journal of Solids and Structures, 6, 1463-1481. [68] Stanley, G. M. (1985), Continuum-Based Shell Elements, Ph.D. Dissertation, Dept. of Applied Mechanics, Stanford University, Palo Alto, California. [69] Sun, QT. and Chin, H. ( 1988), ”Analysis of Asymmetric Composite Laminates,” AIAA Journal, 26, 714-718. 130 [70] Sun, GT. and Chin, H. (1988), ”On Large Deflection Effects in Unsymmetric Cross- Ply Composite Laminates,” Journal of Composite Materials, 22. 1045-1059. [71] Timoshenko, S. and Woinowsky-krieger S. (1959), Theory of Plates and Shells, 2nd ed., McGraw-Hill, New York. [72] Toledano, A. and Murakamim, H. ( 1987), ”A High-order Laminated Plate The- ory with Improved In-plane Responses,” International Journal of Solids Structures, v23n1, 111-131. [73] Yip, Y.C. ( 1996), A 3-D Laminated Plate Finite Element with Zig-zag Sublaminate Approximations for Composite and Sandwich Panels, Ph.D. Dissertation, Dept. of Materials Science and Mechanics, Michigan State University, East Lansing, Michigan. [74] Zaghloul, M.S.A. (1972), Analytical and Experimental Investigations of Filamentary- Composite Laminated Plates, Ph.D. Dissertation, Dept. of Civil Engineering, Univer- sity of Windsor, Windsor, Ontario, Canada. [75] Zaghloul, M.S.A. and Kennedy, J.B. (1974), ”Material Constants of Filamentary- Composite Laminates,” Aeronautical Journal, v78n766, 464-467. [76] Zaghloul, M.S.A. and Kennedy, J.B. (1975), ”Nonlinear Behavior of Symmetrically Laminated Plates,” Journal of Applied Mechanics, 234-236. [77] Zaghloul, M.S.A. and Kennedy, J .B. (1975), ”Nonlinear Behavior of Unsymmetrically Laminated Plates,” Journal of the Engineering Mechanics Division, 169-186. [78] Zhu, Y. and Zacharia, T. (1996), ”A New One-point Quadrature, Quadrilateral Shell Element with Drilling Degrees of Freedom,” Computer Methods in Applied Mechanics and Engineering, 136, 165-203. [79] Zienkiewicz, 0.0. and Taylor, R.L. (1989), The Finite Element Methods - Volume 1: Basic Formulation and Linear Problems, 4th edition, McGraw-Hill International (UK). [80] Zienkiewicz, 0.0. and Taylor, R.L. (1991), The Finite Element Methods - Volume 2: Solid and Fluid Mechanics, Dynamics and Non-linearity, 4th edition, McGraw-Hill International (UK). [81] Zienkiewicz, O.C., Too, J. and Taylor, R.L. (1971), ”Reduced Integration Technique in General Analysis of Plates and Shells,” International Journal for Numerical Methods in Engineering, 3, 275-290. [82] Chiang, H.Y., private communication. [83] Hill. L., private communication. 131 [84] Lee, Y.K., private communication. [85] Nagtegaal, J ., private communication. "illiiliii]17111117111“