aka. 9...... . . .4. . .twnuura 7 1...... cm: 9.; It“. "2 ~ l!Till"!Ill/ll/l‘i/l/Vl/lillill/flit , 3 1293 01567 1356 This is to certify that the dissertation entitled A 3-D Laminated Plate Finite Element with Zig-Zag Sublaminate Approximations for Composite and Sandwich Panels presented by Yuen Cheong Yip has been accepted towards fulfillment of the requirements for Ph.D. Mechanics degree in mew Major professor Date ‘/’ 4/96 MSU is an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY ' Mlchlgan State Unlverslty PLACE I! RETURN BOX to remove this checkout from your mood. To AVOID FINES Mum on or boton dot. duo. DATE DUE DATE DUE DATE DUE l l l l l 1. l ... #- i________i.__lE:l \l I . \ MSU It An Affirmative ActloNEquol Opportunity IMRNOH Wows-m Jl TL l A 3-D LAMINATED PLATE FINITE ELEMENT WITH ZIG-ZAG SUBLAMINATE APPROXIMATIONS FOR COMPOSITE AND SANDWICH PANELS By Yuen Cheong Yip 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 1996 ABSTRACT A 3-D LAMINATED PLATE FINITE ELEMENT WITH ZIG-ZAG SUBLAMINATE APPROXIMATIONS FOR COMPOSITE AND SANDWICH PANELS By . Yuen Cheong Yip This dissertation describes the development of a new zig-zag theory and an eight-node brick finite element using zig-zag sub-laminate approximations which are suitable for the analysis of both thick and thin laminated composite plates and sandwich panels. The zig- zag theory employs the sublaminate concept, in which each computational layer (or sublaminate) contains several, even many, physical layers. The new element derived from this new theory has the topology of a three-dimensional (3-D) continuum-type element (an eight-node brick) although it is based on plate kinematics. This eight-node brick topology permits the stacking of elements through-the-thickness of a laminate, which is not convenient in conventional plate/shell elements whose kinematics are usually with respect to the reference mid-surface. This allows sub-layering for through-the-thickness refinement if higher accuracy is desired. This zig-zag sublaminate approach has the following desirable properties: (i) satisfies interlaminar transverse shear stress continuity; (ii) satisfies transverse shear tractions at the top and bottom surfaces of the plate exactly; (iii) tends to the correct prescribed load value for the transverse normal tractions at the top and bottom with refinement in the thickness direction; (iv) has a small and fixed number of degrees-of-freedorn per sublaminate; (v) uses traditional engineering degrees-of-freedom (displacements and rotations); and (vi) permits the use of adaptive techniques for through-the-thickness discretization when used with the eight-node brick topology. In thin plate applications, the parasitic shear locking phenomenon has plagued shear deformable elements Since their inception. This new eight-node brick finite element using zig-zag sub-laminate approximations utilizes the assumed strain field approach to extend the element to thin plate applications. The element shear strain fields meet both field- consistency and edge-consistency requirements. Numerical results demonstrate that this new element is robust, accurate and computationally efficient. The element also passes the membrane and Kirchhoff patch tests for plates. The current models Show excellent promise for efficient and accurate analysis of thick laminated composites and sandwich panels. Furthermore, the eight-node brick topology will permit the use of adaptive techniques for through-the-thickness discretization and will allow the coupling of the new elements to conventional 3D continuum-based elements. to my beloved wife, Regina Wai Cheng Hooi and daughter, Michelle Shi Yun Yip. iv ACKNOWLEDGENIENT S I would like to thank NASA Langley Research Center for their extensive support for this research under Grant NAG-l-159l. Dr. Jerry Housner and Dr. Alex Tessler were the grant monitors. A special thanks to Dr. Alex Tessler for his assistance in some of the work. I would like to acknowledge my doctoral committee members, Professors Thomas J. Pence, Dahsin Liu, Nicholas J. Alterio and Zhengfang Zhou for their guidance and suggestions during the course of the research. I must also thank the Computational Mechanics Research Group (CMRG) for their helpful discussions and contributions to my work, with special thanks to my colleagues Venkateshwar Rao Aitharaju, Yong-Bae Cho and Lorenzo Marco Smith for their friendship and useful discussions this past year. I am grateful to my wife, Regina for her understanding and support throughout the years that I have known her and to Professor Anthony Waas for if not for him, all this would not be possible. Finally, I must acknowledge my advisor, Professor Ronald C. Averill. In the preparation of this dissertation and throughout the whole period of my studies at Michigan State, I have had the good fortune to benefit from his guidance, advice and encouragement. I will always be grateful and indebted to him for giving me the opportunity as well as the financial support to pursue this study. TABLE OF CONTENTS LIST OF TABLES ....................................................................................................... viii LIST OF FIGURES ..................................................................................................... ix CHAPTER 1 INTRODUCTION AND BACKGROUND ................................................................. 1 1.1 Introduction ................................................................................................... 1 1.2 Background: Laminated Plate/Shell Theories .............................................. 4 1.3 Background: Laminated Plate/Shell Finite Element Models ........................ 9 1.4 The Present Study ......................................................................................... 11 1.5 Interpolation Theory: The C-Concepts ......................................................... 12 1.6 C0 versus C Finite Elements ..................................................................... 14 1.7 Interdependent Interpolation ......................................................................... 14 1.8 Organization of the Dissertation ................................................................... 15 CHAPTER 2 LAYERWISE BEAM ELEMENT ............................................................................... 17 2.1 Introduction ................................................................................................... 17 2.2 Theory Formulation ...................................................................................... 18 2.3 Finite Element Model ................................................................................... 24 2.4 Numerical Studies ......................................................................................... 27 2.5 Summary ....................................................................................................... 31 CHAPTER 3 MIN DLIN PLATE ELEMENT .................................................................................... 3S 3. 1 Introduction ................................................................................................... 35 3.2 Interdependent Interpolation ......................................................................... 38 3.3 Finite Element Model ................................................................................... 43 3.4 Consistent Shear Strain Fields ...................................................................... 45 3.5 Edge Consistency .......................................................................................... 47 3.6 Numerical Studies ......................................................................................... 50 3.7 Summary ....................................................................................................... 63 CHAPTER 4 LAYERWISE PLATE ELEMENT - THICK PANELS ............................................... 65 4.1 Introduction ................................................................................................... 65 vi 4.2 Theory Formulation ...................................................................................... 66 4.3 Theory and Finite Element Model ................................................................ 72 4.4 Numerical Studies ......................................................................................... 80 4.5 Summary ....................................................................................................... 87 CHAPTER 5 LAYERWISE PLATE ELEMENT - THIN PANELS .................................................. 100 5. 1 Introduction ................................................................................................... 100 5.2 Formulation ................................................................................................... 102 5.3 Interdependent Interpolation ......................................................................... 102 5.4 Field Consistent Shear Strain Fields ............................................................. 103 5.5 Consistent Penalty Constraints ..................................................................... 109 5.6 Consistent Transverse Normal Strain Field .................................................. 111 5.7 Edge Consistency .......................................................................................... 113 5.8 Numerical Studies ......................................................................................... 117 5.9 Summary ....................................................................................................... 124 CHAPTER 6 NON-LINEAR MODEL .............................................................................................. 141 6. 1 Introduction ................................................................................................... 141 6.2 Theory Formulation ...................................................................................... 142 6.3 Newton-Raphson Iteration Method ............................................................... 147 6.4 Numerical Results ......................................................................................... 150 6.5 Summary ....................................................................................................... 152 CHAPTER 7 CONCLUSIONS .......................................................................................................... 157 7.1 Review .......................................................................................................... 157 7.2 Area of Future Research ............................................................................... 160 BIBLIOGRAPHY .................................................................................................. 162 APPENDICES APPENDIX A - Inplane Displacement Field Constants (Layerwise Beam) ......... 171 APPENDIX B - Interdependent Interpolation ....................................................... 173 APPENDD( C - Shape Functions .......................................................................... 180 APPENDIX D - Transformed Shape Functions .................................................... 185 APPENDIX E - Interpretation of Stress and Strain Tensors in the Case of Small Strains/Large Rotations ...................................................... 186 APPENDIX F - Element Stiffness Matrix (von Karman strains) .......................... 190 APPENDIX G - Direct Stiffness Matrix (von Karman strains) ............................. 192 APPENDIX H - Tangent Stiffness Matrix (von Karman strains) .......................... 201 vii Table 1 Table 2 Table 3 Table 4 Table 5 Table 6 Table 7 Table 8 Table 9 Table 10 Table 11 Table 12 Table 13 Table 14 LIST OF TABLES Material Properties Used in the Numerical Analyses ............................ 28 Description of Laminate Sequence ........................................................ 29 Straight Cantilever Beam ....................................................................... 52 Curved Cantilever Beam ........................................................................ 5 3 Long Cantilever Beam Test .................................................................... 54 Results for Rectangular Plate, Simple Supports: Uniform Load ........... 55 Results for Rectangular Plate, Clamped Supports: Center Load ........... 56 Normalized Center Deflection: Uniform load ........................................ 63 Lamination Scheme for Random Five-Layer Panel ............................... 82 Material Properties for Random Five-Layer and Sandwich Panel ......... 82 Lamination Scheme for Sandwich Panel ............................................... 83 Material Properties for TACOM’S Composite Armored Vehicle (CAV) Panel ........................................................................................... 84 Results for Square Plate, Simple Supports: Sinusoidal Load ................ 86 Results for Square Plate, Simple Supports: Sinusoidal Load ................ 122 viii Figure 1 - Figure 2 - Figure 3 — Figure 4 - Figure 5 - Figure 6 - Figure 7 - Figure 8 - Figure 9 - Figure 10 - Figure 11- Figure 12 - Figure 13 - LIST OF FIGURES Schematic of Sublaminate and Layer Divisions in the Laminate Theory .................................................................................................... 18 Nodal Topology of Unconstrained Element (Top) and Constrained Element (Bottom) .................................................................................. 27 Schematic of Loading and Support Conditions ..................................... 27 Predicted Normalized Maximum Deflection and Normalized Maximum Inplane Normal Stress versus Number of Elements in the Mesh ....................................................................................................... 32 Normalized Center Deflection versus Span-to—Thickness Ratio of a Five-Layer Simply-Supported Beam Subjected to a Sinusoidal Load .. 32 Axial Displacement versus Normalized Thickness Coordinate at the End of a Simply-Supported Beam Subjected to a Sinusoidal Load ....... 33 Axial Stress versus Normalized Thickness Coordinate at the Midspan of a Simply-Supported Beam Subjected to a Sinusoidal Load .............. 33 Transverse Shear Stress versus Normalized Thickness Coordinate at x=L/10 of a Simply-Supported Beam Subjected to a Sinusoidal Load.. 34 Transverse Deflection versus Normalized Thickness Coordinate at Midspan of a Simply-Supported Beam Subjected to a Sinusoidal Load ....................................................................................................... 34 Sign Convention for Rotations and Stress Resultants ............................ 40 Topology of the Finite Element Model .................................................. 44 Quadrilateral Element Coordinate Description ...................................... 44 Kirchhoff Patch Test for Plates .............................................................. 51 ix Figure 14 - Figure 15 - Figure 16 - Figure 17 - Figure 18 - Figure 19 - Figure 20 - Figure 21 - Figure 22 - Figure 23 - Figure 24 - Figure 25 - Figure 26 - Figure 27 - Figure 28 - Figure 29 - Figure 30 - Straight Cantilever Beam ....................................................................... 51 Curved Cantilever Beam ........................................................................ 53 Long Cantilever Beam ........................................................................... 53 Geometry, Material Properties, Boundary Conditions and Theoretical Solutions for the Rectangular Plate Tests .............................................. 55 Normalized Center Deflection versus Span-to-Thickness Ratio of a Simply-Supported Square Plate Under Sinusoidal Loading (a/b = 1)... 57 Normalized Center Deflection versus Span-to-Thickness Ratio of a Simply-Supported Square Plate Under Sinusoidal Loading (a/b = 5)... 57 Normalized Maximum Normal Stress versus Span-to-Thickness Ratio of a Simply-Supported Square Plate Under Sinusoidal Loading (a/b = 1) ................................................................................... 58 Normalized Maximum Normal Stress versus Span-to—Thickness Ratio of a Simply-Supported Square Plate Under Sinusoidal Loading (a/b = 5) ................................................................................... 5 8 Clamped Circular Plate (Quadrant) ........................................................ 59 Normalized Deflection (W = w 161tD/ (Pa4)) Along Line of Symmetry (y=0) of a Thin (a/h=50) Clamped Circular Plate Under Center Load . 60 Normalized Deflection (W = w161tD/(Pa4)) Along Line of Symmetry (y=0) of a Thick (a/h=2.5) Clamped Circular Plate Under Center Load 60 Normalized Deflection (W = w64D/(qa2)) Along Line of Symmetry (y=0) of a Thin (a/h=50) Clamped Circular Plate Under Uniform Load 61 Normalized Deflection (W = w640/(qa2)) Along Line of Symmetry (y=0) of a Thick (a/h=2.5) Clamped Circular Plate Under Uniform Load ....................................................................................................... 61 Morley’s Acute Skew Plate .................................................................... 62 Schematic of Sublaminate and Layer Divisions .................................... 67 Topology of the LZZ3 Finite Element Model ........................................ 80 Schematic of Loading and Boundary Conditions for Example Problems ................................................................................................ 81 Figure 31 - Figure 32 — Figure 33 - Figure 34 - Figure 35 - Figure 36 - Figure 37 - Figure 38 - Figure 39 - Figure 40 - Figure 41 - Axial Displacement versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Random 5-Layer) Subjected to Sinusoidal Load ..................................................................................... 89 Transverse Deflection versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Random 5-Layer) Subjected to Sinusoidal Load ..................................................................................... 89 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Random 5-Layer) Subjected to Sinusoidal Load ..................................................................................... 90 Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Random 5-Layer) Subjected to Sinusoidal Load ..................................................................................... 90 Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Random 5-Layer) Subjected to Sinusoidal Load ..................................................................................... 91 Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Random 5-Layer) Subjected to Sinusoidal Load ..................................................................................... 91 Axial Displacement versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ....................................................................................................... 92 Transverse Deflection versus Normalized Thickness Coordinate of a Simply—Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ....................................................................................................... 92 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ....................................................................................................... 93 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load - Top Face Sheet ........................................................................... 93 Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 94 xi Figure 42 - Figure 43 - Figure 44 - Figure 45 - Figure 46 - Figure 47 - Figure 48 - Figure 49 - Figure 50 - Figure 51 - Figure 52 - Figure 53 - Figure 54 - Figure 55 - Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 94 Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load - Top Face Sheet .......................................................... 95 Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 95 Axial Displacement versus Normalized Thickness Coordinate of a Simply~Supported Square Plate (CAV) Subjected to Sinusoidal Load . 96 Transverse Deflection versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load . 96 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load . 97 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load - Outer Core .............................................................................................. 97 Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... 98 Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load . 98 Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... 99 Quadrilateral Element Coordinate Description ...................................... 103 Patch Test for Plates ............................................................................... 1 18 Clamped Circular Plate (Quadrant) ........................................................ 119 Normalized Deflection (W = wl61tD/ (Pa4)) Along Line of Symmetry (y=0) of a Clamped Circular Plate Under Center Load ......................... 120 xii Figure 56 - Figure 57 - Figure 58 - Figure 59 - Figure 60 - Figure 61 - Figure 62 - Figure 63 - Figure 64 - Figure 65 - Figure 66 - Figure 67 - Normalized Deflection (W = w640/(qa2)) Along Line of Symmetry (y=0) of a Clamped Circular Plate Under Uniform ............................... 120 Normalized Center Deflection versus Span-to—Thickness Ratio of a Simply-Supported Square Plate (Isotropic) Subjected to Sinusoidal Load ....................................................................................................... 125 Normalized Center Deflection versus Span-to—Thickness Ratio of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ....................................................................................................... 125 Normalized Center Deflection versus Span-to-Thickness Ratio of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load . 126 Axial Displacement versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 126 Transverse Deflection versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 127 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 127 Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load - Top Face Sheet .......................................................... 128 Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 128 Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 129 Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load - Top Face Sheet .......................................................... 129 Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... 130 xiii Figure 68 - Figure 69 - Figure 70 - Figure 71 - Figure 72 - Figure 73 - Figure 74 - Figure 75 - Figure 76 - Figure 77 - Figure 78 - Axial Displacement versus N orrnalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Transverse Deflection versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load - Outer Shell .................................................................................. Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ..................................................................................... Axial Displacement versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... Transverse Deflection versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load - Top Face Sheet .......................................................... xiv 130 131 131 132 132 133 133 134 134 135 135 Figure 79 - Figure 80 - Figure 81 - Figure 82 - Figure 83 - Figure 84 - Figure 85 - Figure 86 - Figure 87 - Figure 88 - Figure 89 - Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (Sandwich) Subjected to Sinusoidal Load ..................................................................................... Axial Displacement versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Transverse Deflection versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Inplane Normal Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load - Outer Shell .................................................................................. Transverse Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Inplane Shear Stress versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Transverse Normal Strain versus Normalized Thickness Coordinate of a Simply-Supported Square Plate (CAV) Subjected to Sinusoidal Load ....................................................................................................... Load-Deflection Curves for the Bending of a Simply-Supported, Square, Orthotropic Plate under Uniform Transverse Loading - Effect of Inplane Mesh Refinement .................................................................. XV 136 136 137 137 138 138 139 139 140 140 154 Figure 90 - Figure 91 - Figure 92 - Figure 93 - Figure 94 - Figure B.1- Load-Deflection Curves for the Bending of a Simply-Supported, Square, Orthotropic Plate under Uniform Transverse Loading - Effect of Through-Thickness Mesh Refinement ............................................... Load-Deflection Curves for the Bending of a Simply-Supported, Square, Orthotropic Plate under Uniform Transverse Loading - Comparison Among Different Theories and Experimental Results ...... Load-Deflection Curves for the Bending of a Clamped, Square, Symmetric Cross-Ply (0/ 90/ 90/ 0) Plate under Uniform Transverse Loading - Effect of Inplane Mesh Refinement .................... Load-Deflection Curves for the Bending of a Clamped, Square, Symmetric Cross-Ply (0/ 90/ 90/ 0) Plate under Uniform Transverse Loading - Effect of Through—Thickness Mesh Refinement. Load-Deflection Curves for the Bending of a Clamped, Square, Symmetric Cross-Ply ( 0/ 90/ 90/ 0) Plate under Uniform Transverse Loading - Comparison Among Different Theories and Experimental Results ............................................................................. Quadrilateral Element Coordinate Description ...................................... xvi 154 155 155 156 156 174 Chapter 1 INTRODUCTION AND BACKGROUND 1.1 Introduction AS analytical solutions are usually restricted to problems with simple geometry, loading and boundary conditions, the finite element method is therefore useful when dealing with complex problems. Using conventional three-dimensional (3-D) displacement-based finite elements is very computationally expensive. Yet in the analysis of composite laminates where local effects are important (e.g. embedded delaminations, ply—dropoffs, free edges, cracks or other regions of 3-D stress fields), one must use a theory based on 3-D kinematics. Most two-dimensional (2-D) plate/shell models are therefore inadequate. However, with the more recent advent of such 2-D (quasi-3-D) theories like the layerwise (e.g., Reddy (1987), Lu and Liu (1992), Liu (1995), Li and Liu (1996)) and zig-zag theories (e.g., DiSciuva (1985, 1987), Cho and Parmerter (1993)), these local effects are modeled satisfactorily. We can implement these more precise theories to develop a robust finite element model that is more computationally efficient than the conventional displacement-based continuum finite element model and yet describe the local kinematics adequately. The zig-zag theory has several important features that are ideally suited for our purpose. Firstly, it possesses only a fixed number of degrees of freedom, regardless of the pr: COI “It (‘03 mm mm number of layers in the laminate. Secondly, it is accurate and satisfies conditions such as (i) the continuity of transverse Shear stress through the thickness of the laminate and (ii) the homogeneous shear traction boundary conditions. Although these features imply an almost ideal combination of accuracy and efficiency, there is one significant drawback to this theory. The a prior assumptions of the displacement field require the transverse deflection degree-of-freedom to be C 1 continuous. This restricts the choice of element interpolation functions that can be used and is therefore very inconvenient for the development of general purpose finite element models based on these theories. The focus of this research will be on the development of a robust C0 finite element model based primarily on the high-order zig-zag theory. It incorporates ideas from the layerwise theories (Reddy, 1987; Barbero and Reddy, 1990; Li and Liu, 1996), the sub- laminate analysis concept of Flanagan (1994), in conjunction with an interdependent interpolation scheme of Tessler and Hughes (1983, 1985). The plate element will be developed in a consistent approach that features exact satisfaction of boundary and interface conditions for displacements and shear tractions. Averill and Yip (1996a) presented a generalized form of the high-order zig—zag theory for beams in which the continuity of the transverse deflection was relaxed by introducing a new variable along with an appropriate constraint condition. Based on this generalized form of the theory, a two-node beam element was developed that contained only C0 variables. The constraint condition was then enforced via a combination of analytical and numerical means by utilizing an interdependent interpolation scheme (Tessler and Dong, 1981) and the penalty method, resulting in a simple, efficient and robust finite element beam model. The new proposed plate finite element model will be based on modified forms of the high-order zig-zag theories. These elements contain only engineering-type degrees of freedom, displacements and rotations, and have the nodal topology of an eight-node brick element, though they are based on plate kinematics. This, therefore, allows our proposed finite element great versatility to model laminates of complex cross-sections. Depending on the technique used in modeling the laminate, this finite element formulation will enable us to achieve a solution similar to that of either the high-order zig-zag or layerwise theory. If all the layers in the laminate are modeled as an equivalent mathematical layer (or block), the current model becomes equivalent to the high-order zig-zag theories (DiScuiva, 1987; and Averill and Yip, 1996b). If the number of mathematical layers is equal to the number of physical layers in the laminate, it will yield the conventional discrete layerwise theory (Reddy, 1987 and Li and Liu, 1996). This type of plate element permits vertical stacking through the thickness of a laminate with all the proper continuity conditions except for transverse normal traction. In addition to being able to emulate the high-order zig-zag theory and layerwise theory, this unique form of the proposed plate element can also be utilized to emulate a concept very similar to the “local and global laminate” model (sub-laminate analysis) proposed by Pagano and Soni (1983) to model a laminate consisting of many layers. They came up with a global and local model whereby they divided the laminate into two parts; namely, (i) the local region - region that is of interest and (ii) the global region - the remaining region. The field equations in their model are based upon an assumed thickness distribution of stress components within each layer of the local region and an assumed thickness distribution of displacement components in the global region (effective properties are used). Although an apparent loss of accuracy occurs in the calculation of transverse Shear stress components Nil Na: 8581 [O b Star in the global model, increasing the number of sub—laminates will obviously enhance both the accuracy of transverse stress predictions as well as the local predictions. The type of sub-structuring technique described above is especially important when solid modeling of individual layers may not be feasible yet it is important to isolate a region of critical behavior. This efficient local-global modeling of composite structures has been employed widely (see for example, Lee and Liu, 1991; Kong and Cheung, 1995; Yu et al., 1995; Li and Liu, 1996a). 1.2 Background: Laminated Plate/Shell Theories Several recent survey papers summarize the approaches of plate theories (Bert, 1984; Noor and Burton, 1989, 1990; Kapania and Raciti, 1989; Reddy, 1990). They consider the problem of reducing the three-dimensional equations of the theory of elasticity to two- dimensional plate/shell equations. The two-dimensional theories are obtained from the three-dimensional elasticity theory by making a prior assumptions concerning the variation of displacements and/or stresses through the thickness coordinate of the laminate, yielding an equivalent single-layer plate theory. These theories based on the method of expansion with respect to the thickness coordinate can be found in the work of Mindlin (1961). The classical laminate plate theory (CLPT) is an extension of the classical plate theory (CPT) to laminated plates. In both theories, the inplane displacements are assumed to vary linearly through the thickness and the transverse displacement is assumed to be constant through the thickness (i.e., transverse normal strain is zero) - Reissner and Stavsky (1961), and Dong et al. (1962). With the anisotropic and non-homogeneous nature of composite laminates, important effects of transverse shear and normal deformation cannot be neglected, demanding improvements in the classical laminate plate at: 10;. her theory (CPLT), which is based on the well-known Kirchhoff assumptions. As a result, modeling approaches which take into account these effects have been the topic of research in the last two decades. A refinement to the classical laminate theory is provided by the first-order Shear deformation plate theory (FSDT). The first-order shear deformation theory is commonly known as the “Mindlin plate theory”, (Mindlin (1961) and Reissner (1945)). The first- order shear deformation theory yields a constant value of transverse shear strain through the thickness of the plate, and thus requires shear correction factors. Inclusion of first- order transverse shear effects is often adequate as shown by Yang et al. (1966), provided the laminate does not contain adjacent layers with drastically varying stiffness properties (see Averill, 1994; Sun and Whitney, 1973). Second and high-order shear deformation plate theories (HSDT) involve high-order expansions of the displacement field (Reddy, 1984; and Barbero and Reddy, 1990). Barbero and Reddy (1990) suggested a generalized strain-consistent third-order theory and demonstrated that the displacement fields of many other “equivalent single layer” third-order theories can be derived from it. The displacement field of HSDT accommodates the vanishing of transverse shear strains (and hence shear tractions) on the top and bottom of a laminate. This type of theory, however, generally results in incompatible transverse shearing stresses between adjacent layers. In summary, all the above theories (CLPT, FSDT, HSDT) differ only in their choices of assumed displacement fields. Termed as “equivalent single layer” approach, all those theories have a common drawback in that the thickness variation of displacements, and . . . . 1 . . thus strains, rs assumed to be continuous and smooth (1.6. C continuous). This characteristic will lead to a discontinuous interlaminar stress field because of different constitutive properties in adjacent plies. It therefore precludes the satisfaction of continuity of transverse stresses at interfaces between adjacent layers of different materials, and does not accurately reflect the kinematics in laminates. This may become Significant for laminates that contain drastically different constitutive properties in adjacent layers. To overcome the deficiencies of “equivalent single layer” theories, discrete-layer (or layerwise) theories (DLT) have been developed in which a unique displacement field is assumed within each layer (see for example, Reddy, 1987; Toledano and Murakami, 1987; Lu and Liu, 1992). These theories explicitly account for the layerwise construction of laminates and allow the possibility of achieving continuous transverse shear stresses at all layer interfaces. The layerwise theory by Reddy (1987) treats each layer within the multi-layered laminates as separate individual layers or combines several layers into sub-laminates. He assumes that the displacement field is expanded within each layer using Lagrangian shape functions used in conventional elasticity displacement finite elements. It is therefore a quasi-3-D model based on plate kinematics as opposed to truly 3-D kinematics. Each layer of the theory can be treated either as a mathematical or a physical layer. Several physical layers can be combined into a sub-laminate and treated as a mathematical layer wherever necessary (for computational savings). Because of the layerwise construction of this theory, it can therefore model local effects like delaminations, ply terminations (i.e. ply- dropoffs), ply splits, etc. The most significant aspect of this layerwise theory is that, compared to the conventional 3-D displacement finite element model, it has a data Vii DOI mu [0 5 kill: [nor COns C0m structure that saves computational time while giving exactly the same results for comparable meshes (see Robbins and Reddy, 1992). However, this theory contains a number of degrees-of-freedom that is proportional to the number of layers in the laminate. Consequently, for analysis of thick-section laminates involving hundreds of layers, the computational effort required will still be enormous. The number of degrees of freedom will grow rapidly with the number of layers and computer storage may force limitations on the total number of layers analyzed. The sub-laminate concept within the layerwise theory by Reddy (1987) that can be employed to model several layers as one equivalent single layer is an attractive feature to save computational costs and efforts. However, one drawback of this concept when several layers are combined is that the continuity of interlaminar shear stress conditions will be violated. Liu and his co-workers (1995, 1996) extended the work of Reddy (1987) by adding nodal rotations in addition to nodal displacements as independent variables. The nodal rotations, however, are not required to be continuous across the composite interface so as to satisfy interlaminar shear stress continuity, reflecting the correct nature of the laminate kinematics. The theory gives excellent results for both stresses and displacements and, more importantly, allows accurate predictions of transverse shear stresses from the constitutive equations instead of having them recovered from the equilibrium equations as commonly done. The discrete layerwise theory (DLT) developed using plate kinematics requires only 2- D finite elements. The advantage is that the element aspect ratio for these layerwise models is restricted only in the two dimensions (in-plane coordinates), as opposed to all three dimensions in the conventional elasticity finite elements. It is required for 3-D displacement finite elements to have an element aspect ratio of one or close to one for good converged results. To achieve this ratio, the thinness of each layer in the laminate will force the other two dimensions of an element to be small in the case of 3-D elements. Consequently, a very large number of 3-D elements will be needed to ensure a proper mesh There is another new class of laminate theory termed the zigzag theories which incorporates the layerwise construction (different constitutive properties in each layer) of laminates into the in-plane displacement kinematics. This was done by assuming piecewise continuous functions through-the-thickness for in-plane displacement components. These theories only involve a small number of degrees of freedom (depending on the order of approximation of the displacement fields) unlike the layerwise theories discussed above where the number of degrees of freedom is a function of the number of layers in the laminate. This reduction in the number of degrees of freedom circumvents the overwhelming difficulties and complexities associated with stress analysis of multi-layered composite laminates. The original idea can be found in the earlier works of DiSciuva (1985) where he assumed a global linear variation superposed on a piecewise linear variation through the thickness of the in-plane displacement. This theory is termed a first-order zig-zag theory (FZZT). Although the theory enforces interlaminar shear stress continuity at each individual interface analytically, the transverse shear stresses are constant through the thickness because of the low order assumed for the displacement field. In addition, it does not satisfy the traction conditions at the top and bottom surfaces. This drawback was overcome by extending the theory to include higher-order terms for the in-plane displacement field (see for example, Murakami, 1986; DiScuiva, 1987, 1992; Lu and Liu, 1992; Cho and Parmerter, 1993; Liu, 1995; Averill and Yip, 1996a). Refinement to the first-order zig-zag theory was done by superposing a global cubic function instead of a linear variation through the thickness onto the zig—zag function for the in-plane displacement component. This improved the prediction of structural response in thick, and especially, asymmetrical laminated structures (the form of DiSciuva’s theory (1993) is strictly valid only for symmetric laminates). These are termed high-order zig-zag theories (HZZT). Several other investigators have followed up with variations and slight improvements on the above zig-zag theories (Xavier, Lee and Chew, 1993; Ling-Hui, 1994; and Murakami, 1986). 1.3 Background: Laminated Plate/Shell Finite Element Models Modeling of thick-section laminates using conventional continuum-type finite elements normally involves extremely large numbers of these elements to achieve accurate predictions, and is therefore highly computationally expensive. Structural finite element models derived using first-order (FSDT) as well as high-order (HSDT) shear deformation theories (see Reddy, 1987) do not face this problem, but they are often inaccurate because of two factors. First, they often do not adequately represent the effects of transverse shear stresses and transverse normal stresses. Second, the commonly assumed linear variation of inplane displacement components through-the-thickness of the laminate is only valid in an average sense and is inadequate for accurately predicting the response of these structures, especially in critical damage analysis. Higher order effects are therefore needed for improved structural analysis. Bogdanovich (1991) presented a detailed discussion of the 10 various shortcomings of these computational models used in the analysis of thick laminated composite plates. A detailed discussion of various theories and computational models for composite laminates can also be found in Reddy and Robbins (1994). Discrete layerwise theory (Reddy, 1987; Barbero and Reddy, 1990; Robbins et al., 1991; Lu and Liu, 1992; Robbins and Reddy, 1993; Liu, 1995; Li and Liu, 1996) and layerwise constant shear theories (see Srinivas, 1973; Epstein and Glockner, 1977; Epstein and Huttelmaier, 1983; Barbero et al., 1990) were introduced to overcome the limitations of shear deformation theories (FSDT and HSDT). These theories are based on a distribution of displacements that is continuous through-the-thickness, while derivatives with respect to the thickness coordinate are not necessarily continuous. Discontinuous derivatives will result in discontinuous transverse Strains at layer interfaces. This in turn will allow the satisfaction of interlaminar transverse stress continuity. These theories are very accurate, but the “three—dimensional” elements derived from such theories generally have a large number of degrees-of-freedom associated with them. This is because the degrees-of-freedom are directly related to the number of physical layers in the laminate. In fact, the elements of the discrete layerwise theories of Reddy (1987, 1994) can be viewed as a ‘super-element’ of an assembly of continuum elements. These discrete layerwise theories can be used in conjunction with other less refined theories for more efficient local-global modeling of composite structures (see Lee and Liu, 1991; Kong and Cheung, 1995; Yu et al., 1995). Other notable recent “three-dimensional” elements for laminate analysis include the layerwise constant shear element of Barbero (1991) and Zinno and Barbero (1994) which have the form of an eight-node brick element with constant transverse shear stresses and deflection through-the-thickness. This ll topology allows the flexibility of mesh refinement through the thickness, has the versatility of modeling complex Shape structures, e.g. ply drop-off, lap—joints, tapered beams, etc., and also maintains compatibility with conventional three-dimensional continuum elements. Other models include those of Babuska et al. (1992) and Actis and Szabo (1993) which use the hierarchical approach and Bogdanovich (1991) which uses piecewise-polynomial deficient Spline functions to analyze composite laminates. 1.4 The Present Study This dissertation will introduce a new technical theory and associated finite element model for analysis of thick laminated composite and sandwich panels. This theory essentially combines the discrete layerwise and zig-zag theories mentioned above. The theory employs the sublaminate concept, in which each computational layer (or sublaminate) contains several, even many, physical layers. Within each sublaminate, an accurate approximation of the displacement field is employed that accounts for discrete- layer effects without increasing the number of degrees-of-freedom as the number of layers is increased. This is accomplished by satisfying analytically the continuity of transverse shear stresses at layer interfaces as well as Shear traction conditions at both top and bottom surfaces of each sublaminate. Because the resulting through-thickness variation of “in- plane” displacements takes the form of a_ piecewise nearly linear function, the theory is still essentially a zig-zag theory. The operative degrees-of-freedom in the theory are located at both the top and bottom surfaces of each sublaminate to facilitate the satisfaction of continuity conditions between sublaminates. Averill and Yip (1996b) have employed this scheme for a layerwise beam with success. This fundamental idea of recasting the proposed displacement field into 12 surface quantities rather than traditional mid-plane quantities without altering the mechanics of the plate problem can be traced to Flanagan (1994). Flanagan (1994) sub- divided a general laminate beam into sub-laminates to determine the strain energy release rates in composites. The accuracy and efficiency of the present theory is thus adaptable, depending on the number of sublaminates chosen to model a given laminate. For most global structural analyses, only one sublaminate is needed through the thickness of the entire laminate to attain the desired accuracy of overall structural response measures. More sublaminates may be used to increase the accuracy of overall structural response predictions or to capture local effects such as interlaminar stresses. Depending upon the ply stacking sequence and the type of global and/or local response measures that are sought, the optimal number of sublaminate approximations required for accurate analysis of thick multilayered composite laminates or sandwich panels will often be greater than one but far less than the number of the layers in the laminate. The adaptable nature of the current theory can thus be used to great computational advantage. 1.5 Interpolation Theory: The C-Concepts We will also introduce the basic concepts needed to understand why some finite elements fail while others succeed when modeling structural mechanics problems. The concepts of continuity and completeness as understood in finite element practice are first introduced. The continuity criterion ensures the compatibility of the displacement fields across element edges. This can be reasonably achieved by ensuring the continuous representation of the displacement fields across element edges if the element is CO. C1 elements like those derived from the Kirchhoff-Love theories of plates and Shells, where 511 hr the cor fon elem sing Shea state demc consi Pilma (01151} C 0112‘ in that [h CflSUre la}: Fa The elemfm. hat the] 13 strains are based on second derivatives of displacement fields, will require continuity of first derivatives of displacement fields in addition to the displacement fields themselves. The completeness criterion restricts the form of the assumed displacement field to ensure the element can produce the strain—free rigid body motion as well as constant strain state condition. These two criteria are seen to be insufficient to provide a complete and coherent formal basis for the displacement type finite element method. The poor behavior of conventionally formulated displacement-type C0 finite elements, generally called ‘locking’, is often attributed to the high rank and non- Singularity of the penalty-linked stiffness matrix. Shear locking is due to the inability of Shear deformable elements to accurately model the curvatures within an element under a state of zero transverse shearing strain (Averill and Reddy (1992)). Prathap (1994) demonstrated that the correct rank and non-singularity required emerges directly from the consistency of the discretized strain field approximations and are not necessarily the primary cause of the poor behavior of the elements. He introduced the concepts of consistency and correctness norms and collectively called them the C-Concepts, namely, continuity, completeness, consistency and correctness. The consistency criterion requires that the interpolation functions chosen to initiate the discretization process must also ensure that any special constraints that are anticipated must be allowed for in a consistent way. Failure to do so causes the model to be ‘locked’ out of the ‘correct’ solution. The C—Concepts of Prathap (1994) will be utilized in the development of our new finite elements based on the sub-laminate zig-zag theory. Numerical results will demonstrate that the new elements developed using this new strategy are very accurate and robust. 14 1.6 6’ versus C‘ Finite Elements The development of elements based on well-known classical theories like the Euler- Bemoulli beam theory, the Kirchhoff-Love plate theory and equivalent shell theories will result in C1 elements. These Cl theories will entail more restriction on the use of approximation functions in the trial space besides having limitations in the range of applicability. Cl continuity will require the continuity of slopes in addition to the continuity of fields across element boundaries. They therefore must rely on cubic or higher polynomials for the element shape or basis functions. AS such it is more attractive to have elements that require C0 continuity and therefore require only simple basis functions to satisfy the continuity of fields across element boundaries. It is with this in mind that our new elements are specially formulated using a penalty-type formulation to require simple 0 . . . C continuous basrs functions. 1.7 Interdependent Interpolation Most conventional displacement type C0 finite elements are prone to “lock” in the thin regime. This arises from the field inconsistency when equal order interpolations are used for field variables which appear in different orders of its derivatives in the strain field that has to be constrained. Thus, if one ensures that the strain field is consistently represented by a proper a priori choice of unequal order interpolations for the contributing field variables, there would not be ‘locking’. Tessler and Doug (1981) and Tessler and Hughes (1983, 1985) derived a family of Timoshenko beam and Mindlin-Reissner plate elements using this approach. However, unequal order of interpolations typically gives rise to elements with nodes having different degrees-of-freedom at each node. In the case of a Timoshenko beam 15 element with three nodes (two end nodes and a mid-node), it will have the deflection degree-of-freedom, w defined at all three nodes and the rotation degree-of-freedom, 0 defined only at the two end nodes. A simple quadrilateral Mindlin plate element will require eight nodes, four corners and four mid-side nodes. The deflection degree-of- freedom, w will again appear in all eight nodes and the rotation degrees-of-freedom, 0x and By only at the four corner nodes. Tessler and co-workers (1981, 1983, 1985) overcame this problem by condensing out the undesirable middle or mid-side nodes by explicitly enforcing the linear Kirchhoff mode in the case of the beam or through the use of four differential edge constraints, i.e. the vanishing of the linear transverse shear strain, in the case of a plate. This will result in an explicit deflection/rotation dependence. In the present formulation, interdependent interpolation will be used both to alleviate locking and to partially satisfy explicit constraints that are introduced to relax continuity of the variables. 1.8 Organization of the Dissertation The dissertation will be organized in such a way that each of the above concepts will be introduced gradually. In Chapter 2, the basic concepts of the new technical theory and the development of the finite element model will be introduced for the one-dimensional beam problem. The beam problem is the simplest case that could be used to demonstrate the underlying principles involved in the new technical theory. Chapter 3 will introduce the interdependent interpolation concepts to alleviate shear locking in Mindlin plate elements. In addition, the field and edge-consistency concepts will be introduced to identify the inconsistent Strain fields that is causing elements to be ‘locked’ out of the solutions for constrained problems e. g. the vanishing of the transverse 16 normal and shear strains in thin plates. Chapter 4 extends the one-dimensional laminate beam theory of Chapter 2 to a two- dimensional laminate theory for plates. In Chapter 5, the independent interpolation scheme and the field and edge-consistency concepts discussed in Chapter 3 will be adapted and used to develop a robust plate finite element model using this new technical theory. Von Karman non-linearities are added to the element model in Chapter 6. Actual experimental test results were used to compare the solutions predicted by the model. The conclusions are in Chapter 7. The chapter summarizes the major differences between the present theory and the existing discrete layerwise and zig-zag theories as well as their derived finite element models. It also includes a section that discusses possible area of research that can be pursued in the future. d: C01 all lift HE“ Chapter 2 LAYERWISE BEAM ELEMENT As a preliminary study, we will first consider finite element approximations of the Simplest structure, i.e. the beam element, using the new laminate theory with zig-zag approximations. This is important as it will clearly illustrate the principle that is used in deriving the new laminate theory. This chapter will also briefly cover the interdependent interpolation concept that is used in beams. 2.1 Introduction A new laminated beam theory is developed that can be considered a hybrid theory, combining the general layerwise theories with independent layerwise kinematic approximations and the efficient zig-zag theories in which the layerwise degrees of freedom are eliminated by enforcing transverse shear stress continuity conditions. In the new theory, the through-the-thickness approximation of the inplane displacement components takes the form of a layerwise theory in which each layer is really a sublaminate containing several, even many, physical layers. Within each sublaminate, a zig-zag through-the-thickness approximation of the inplane displacement components is taken in which the layerwise degrees of freedom are eliminated by enforcing continuity of transverse stresses. Shear traction conditions at the top and bottom of each sublaminate are also satisfied. The theory includes the effects of transverse normal strain, and is cast in a 17 18 form that is exceptionally well-suited for solution by the finite element method. An accurate and convenient four-node planar element with beam-type kinematics is developed and its utility is demonstrated. 2.2 Theory Formulation In the present theory, the layerwise construction of laminated composite beams is modeled as M sublaminates, with each sublaminate containing N m layers, where m is the sublaminate number (see Figure l). Z,z ,6 Pym/ """" 27)), 2 $333fi>§3§3§ 1"” z, ————— Subaunm 1 ----- 4 21g: Figure 1. Schematic of sublaminate and layer divisions in the laminate theory The total number of layers in the laminate is then: M . Ntotal = 2 Nm (1) m- The thickness and material stiffness properties of each layer are arbitrary, and it is assumed that adjacent layers are perfectly bonded together. In the following, 2 is a local (sublaminate) thickness coordinate with its origin at the bottom of the sublaminate, while Z is the global (laminate) thickness coordinate. The coordinate x is measured along the length of the beam. SL su' “he Sub l. h 19 If the width of the beam is small, then 0),, txy, tz U are negligible, and a state of plane stress can be assumed. The plane stress constitutive relations for the kth layer of the beam take the form: - - k k k - . k 0x() ”C(1)C(13)0 l8x() OZ = C(],;) Cg? O 8z (2) _sz_ o 0 C§§’_ .Yxa The infinitesimal strain—displacement relations are: (k) (’0 8m _ 3'11: 8m _ 9:9 x - x z — z (3) (k) (k) (k) a z aux In each sublaminate, an independent displacement field is assumed in which the through-the-thickness variation of inplane displacements is described by a cubic polynomial in the local thickness coordinate with a piecewise linear (or zig-zag) function superimposed upon it. The transverse deflection is assumed to vary linearly with the transverse coordinate. The displacement components of the nth layer within the mth sublaminate, where 1 S n S N m and 1 _<_ m S M , can be written as: “(m ")(x _k20 zku ("1)”)+ ”Z (Z- 2 KO") i=1 (4) ugm’")(x, z) = z)w(m)(x)(l Th Z)+w(m)(x)(h—z— ) where the subscripts b and t refer to the bottom and top surfaces, respectively, of the mth sublaminate, and hm is the total thickness of the mth sublaminate. It is possible to eliminate the degrees of freedom 6m) in Eq. (4) by enforcing the ”an: $11sz “here 20 condition of transverse shear stress continuity at each layer interface. The condition of shear stress continuity at the jth interface is: (mi) ("1.1+ 1) 1,, = I... (5) Making use of Eqs. (2-4) in Eq. (5), it is found that: dw . g, - a,(a, +311”), b,i22 (6) . dwt dwb ,. ,. +C‘(dx —E )+d,u3 where: (i) t—l C55 .. i=[(i+1)_111+ 20]] C55 j=l (I) i-l . C55 . bi =[ (i+l)_112Zi+ 2 b1] C55 j= (I) i—l (7) e.=[C55—1-zi+ a] ‘ (i+l) h 1 C55 ”' J: (i) r—l . C55 2 i = [ (1+1)—II3z,-+ 2d,] C55 j= Additional simplification of Eq. (4) can be achieved through satisfaction of the transverse shear traction boundary conditions at the top and bottom surfaces of each sublaminate. For sublaminate m, these conditions are: ("1.1) (m) T = T xz 2:0 b 8 mm.) _ (m) ( ) xz z= h - Tt "I where 12'"), rim) are the applied Shear tractions (or interlaminar Shear stresses, as the fu sul 21 case may be) at the bottom and top surfaces, respectively, of the mth sublaminate. The displacement field is cast in its final form by introducing the variables: (In) “hm = “5cm, l) = ’70 2:0 9 (m) (mst) ( ) u, = ux z=h,,, Eqs. (8, 9) can be solved either analytically or numerically so that the displacement field for the mth sublaminate can be expressed in terms of the operative degrees of freedom (all functions of It only): (m) (m) e(m) ("0 (m) Wm) 9(m) t ’ t (m) ub ,Wb , b ,Tb ,ut , T , , (10) where again the subscripts b and t refer to the bottom and top surfaces, respectively, of the sublaminate, and (M) (m) (M) __ dwb (In) _ dw, 9b -— ’t_Fx_ dx (11) The displacement field now takes the form: the the 9&0 COU pies Cap} 22 uim’") = u(m)[l +p(lm)z 2+pf,_m)z3 + ”2 (z— 2,1)a "0] i=1 n-l (m) (m) 2 (m) 3 ( ) + u, [—plz Z—pz Z -2 (Z‘Zi)aim ] i=1 +0(m)[—z+p(3m)z+ 2+ pftm)z3 + "2 (z— z )bim)] i=1 +9§m)[p<5m m) z2+ p(m)z3 + 2(z— z)c(-m )] (12) i=1 n—l 2 3 + tgm)[p£,m )z + pgm)z + pamz + 2 (z - Zi)dl(-m)] i=1 .‘A n-l 2 +r§m)[p 'g)z + p111); + 2(2— z)e(.m )] i=1 ("1. n) (In) 2 (m) 2 z hm t hm In Eq. (12), pk (k = l to 11) and al., bi, ci, di, ei are functions of the layer shear stiffnesses and thicknesses, and can thus be calculated a priori. Their functional forms are given in Appendix A. In the sublaminate displacement field described above, the functions of z that multiply the degrees of freedom of (10) can be viewed as shape functions that describe the through- the-thickness variation of these measures within the sublaminate. Thus, a layerwise laminate theory can be developed in the form introduced by Reddy (1987), where now each “layer” in the layerwise theory may contain several, even many, physical layers. Of course, this concept was always possible to envision with the layerwise theories, but the present formulation provides an extremely accurate and efficient approximation within each sublaminate. The sublaminate approximations are connected by imposing continuity dlSl 23 of displacements and transverse shear stresses at the interface between each sublaminate region. The theory thus has the following properties: (i) (ii) (iii) (M (V) continuity of transverse shear stresses through the entire thickness of the laminate is satisfied; shear traction boundary conditions on the top and bottom surfaces of the laminate are satisfied exactly; a piecewise (layerwise) continuous through-the-thickness variation of the inplane displacements is allowed, yet the number of degrees of freedom in each sublaminate is independent of the number of layers in that sublaminate; a linear through-the-thickness variation of the transverse deflection in each sublaminate explicitly accounts for transverse normal strains and stresses; only engineering-type degrees of freedom are used -- displacements and rotations, plus shear traction terms that are always known on the surfaces and unknown between sublaminates. The accuracy and efficiency of the present theory is thus adaptable, depending on the number of sublaminates chosen to model a given laminate. If only one sublaminate is used through the thickness of the entire laminate, then the theory falls into the class of zig-zag theories (see DiSciuva, 1985, 1986, 1987, 1993; Cho and Parmerter, 1993; Xavier, Lee and Chew, 1993; Ling-Hui, 1994; Murakami, 1986; Averill, 1994; Averill and Yip, 1996a). If the number of sublaminates is equal to the number of layers in the laminate, then the theory can be categorized with the more general layerwise theory of Reddy (1987). In this case, a cubic layerwise through-the-thickness variation of the inplane displacement components would yield extremely high accuracy (as in Lu and Liu, 1992). fre file. “her 24 The optimal number of sublaminate approximations required for accurate analysis of thick multilayered composite laminates or sandwich panels will often be greater than one but far less than the number of layers in the laminate. The adaptable nature of the current theory can thus be used to great computational advantage. 2.3 Finite Element Model Because first derivatives of wb and w, appear in the displacement field of Eq. (12), second derivatives of these variables are present in the strain energy functional, requiring them to be C 1 continuous. In the current finite element model, this continuity requirement is weakened by treating 0b and 0, as independent degrees of freedom and subsequently imposing the constraints dwb 8b = 91"}; = 0 (l3) dw, gt = et-E; = 0 via a combination of an interdependent interpolation scheme and the penalty method (Lynn and Arya, 1974; Zienkiewicz, Owen and Lee, 1974; Reddy, 1980). All degrees of freedom in the displacement field are then required to be only C0 continuous. The principle of minimum total potential energy is employed to develop the finite element model. For a constrained system, we have: anp = so + 5v + Sg’qgfijdx] nag 1.3)...) = 0 where U is the internal strain energy, V is the potential energy of external forces, and 71,, y, (14) 25 are penalty parameters that enforce constraints in Eq. (13). Substituting Eqs. (2, 3, 12, 13) into Eq. ( 14), the governing equations, boundary conditions, and displacement-based finite element model can be developed. While the finite element geometry could take the form of a two-node beam element with eight degrees of freedom per node, it is advantageous to use the topology of a four-node, planar element, as shown in Figure 2 (bottom). This topology allows the laminate thickness to be conveniently subdivided and modeled by multiple finite elements (representing sublaminates). It is thus possible to increase the accuracy of the finite element model, as needed, to capture through-the-thickness gradients and transverse (interlaminar) stresses. It is also possible to simulate delaminations using the redundant node concept. The element formulation takes advantage of an interdependent interpolation concept introduced by Tessler and Doug (1981) and recently used by other element developers (Averill, 1994; Averill and Yip, 1996a; Friedman and Kosmatka, 1993). Except for the element interpolation scheme discussed below, all aspects of the finite element formulation follow the well-known standard procedures (see, for example, Reddy, 1993), and the details are omitted. In the constrained element formulation, the transverse deflection degrees of freedom wb, w, are initially approximated using quadratic Lagrange interpolation functions, while all other degrees of freedom are expanded using linear Lagrange interpolation functions. Such an approximation scheme results in a six-node element in which each of the two midside nodes contain a single degree of freedom associated with the transverse deflection (see Figure 2, top). These midside nodes are eliminated by considering a modified form of the constraint conditions in Eq. (13): PM the the C0n. aim the 5 Wk 26 d dwb) _ 2319')? ‘ 0 d dw, _ 210*? l " 0 By substituting the element approximations into Eq. (15), the midside degrees of (15) freedom wb3, w,3 may be determined in terms of the other nodal degrees of freedom, resulting in a new interdependent element interpolation scheme: 2 2 wbs 2 wbij+ Z ijlT/j i=1 i=1 2 2 (16) (42,5 2 wtij+ 2 9,ij i=1 i=1 where P1, P2 are the linear Lagrange interpolation functions, 1711 = (lg-<1 -&2) N2 = Ego—£2) (17) are is the length of element e, and 7; is the natural axial coordinate in the element (-1 5 § 5 1) . If the constraints in Eq. (13) are viewed as being composed of a constant part and a linear (in x) part, then the above interdependent interpolation scheme satisfies the linear part of the constraints identically. This makes the penalty function technique in the present formulation more robust, because the penalty parameter enforces only the constant part of the constraints. This approach effectively increases the order of the element without introducing any additional nodes or degrees of freedom and eliminates the shear locking problem so an exact order of integration may be used. A consistent force vector is also obtained. 2.4 I) SUbje 27 “11’ WI] “:2, Wt2 W 0:1, 1:rl '3 922’ 7:2 O C “bl, wbl W123 “122’ sz 9“. Tb] 9w sz utl’ wtl “,2, "Q2 911’ Ill 9’2, 1’2 E::::i p- _ __..L.:] “bl, Wbi “b2, sz 91w Tb] 9w sz Figure 2. Nodal topology of unconstrained element (top) and constrained element (bottom) 2.4 Numerical Studies Numerical results are presented for bending of a simply-supported laminated beam subjected to a sinusoidally varying transverse load of magnitude one (see Figure 3). Comparisons are made between predictions of models based on first-order shear deformation theory (FSDT) (see Yang, Norris and Stavsky, 1966), the present layerwise zig-zag theory (LZZT), and an exact elasticity solution of Pagano (1969). Figure 3. Schematic of loading and support conditions The material properties used in the analyses are as listed in Table 1, where E0 is a 28 reference value that is chosen here to be 1.0x106. Laminates composed of combinations of these materials have been used previously (Toledano and Murakami, 1987), and provide a good test of the models because they are distinctly different in their stiffness characteristics. Material 1 is considered compliant in tension/compression and compliant in shear. Material 2 is stiff in tension/compression and stiff in shear. Material 3 is stiff in tension/compression and compliant in shear. Table 1. Material properties used in the numerical analyses E/Ea G/Ea Material 1 1.0 0.2 Material 2 32.5 8.2 Material 3 25.0 0.5 The laminate configuration used in the present numerical study is described in Table 2. It is an unsymmetrical five-layer laminate with seemingly random layer thicknesses and material properties that vary drastically from layer to layer. This example provides a very rigorous test of a laminate theory’s ability to model complex laminates subjected to sharply varying loads. While results are presented here for only one laminate, it should be noted that the current model has been tested on many problems with a wide variety of lamination sequences and geometries, including sandwich beams. In each case, comparisons between predictions of the present model and the exact elasticity solution were at least as good as the comparisons presented here. p11 CIT 29 Table 2. Description of laminate sequence Laminate Layer Number Relative Thickness Material Number (Volume Fraction) 0.100 0.250 0.150 0.200 0.300 5 layer random MAWN— DJ—‘UJN— For a thick laminate with span-to-thickness ratio of four, the normalized maximum deflection and the normalized maximum inplane normal stress predicted using a single sublaminate are plotted in Figure 4 versus the number of elements used in the mesh. The results are normalized by the exact solution of Pagano (1969) which has been modified for beam analysis by specifying the appropriate constitutive equations for beams, Eq. (2), instead of those for cylindrical bending. When only four elements are used, the error in both cases is less than four percent. However, it appears that a nearly converged solution is attained using approximately ten elements. Thus, all subsequent finite element results are obtained using a uniform mesh of ten four-node elements (with full integration) for LZZT and ten two-node elements (with interdependent interpolation of Tessler and Doug, 1981) for FSDT. Note in Figure 4 that there is approximately five percent error in the predicted stress when ten elements are used. This error exists because the current model does not predict the exact result for this laminate when the aspect ratio, Uh, equals four. Thus, the error is primarily due to assumptions in the theory as opposed to finite element approximation errors. The accuracy of these predictions increases rapidly as the aspect ratio of the beam increases. In Figure 5, the predicted normalized center deflection versus the span-to—thickness (hi ant Spa “564 the. 30 ratio of the beam is shown. The deflection predictions of FSDT degrade Sharply for aspect ratios less than about 50, largely because the material properties of adjacent layers vary drastically. It should be pointed out that a shear correction factor of 5/6 (valid only for isotropic beams) was used to obtain the FSDT results, which accounts for some of the error. Predictions obtained using LZZT with one element (one sublaminate) through the thickness are accurate for all aspect ratios greater than about two, and no shear correction factor is needed. Due to the interpolation schemes used, there is no locking in either model as the beam thickness (and thus the thickness of each element) decreases. In Figures 6 and 7, the through-thickness distributions of inplane displacement and inplane normal stress, respectively, as predicted by FSDT, LZZT, and Elasticity, are plotted for a thick laminate having a span-to-thickness ratio of four. It can be seen that LZZT, with only one element through the thickness, does an excellent job of predicting the through-thickness distributions of both inplane displacements and stresses, even for this very thick laminate. FSDT is not able to capture these variations. These two plots highlight the importance of including the zig-zag through-thickness variation of inplane displacements (and hence, inplane strains) for laminates in which the material properties of adjacent layers vary considerably. In Figures 8 and 9, the utility of being able to use multiple LZZT elements through the thickness of a laminate is demonstrated for the prediction of transverse Shearing stresses and the through—thickness distribution of the transverse deflection in a thick laminate with span-to-thickness ratio of four. When five LZZT elements (one per physical layer) are used through the thickness, the predictions of transverse shearing stress (calculated using the constitutive equations, Eq. (2)) are indistinguishable from the exact solution, and the fo 31 variation of transverse deflection is very good. Because increasing the number of sublaminates in the model also increases the finite element discretization through-the- thickness of the laminate, a greater number of degrees of freedom is required for such refined analyses. It has been shown, however, that excellent results can usually be obtained using only one element through the thickness of a laminate. Further, because inplane normal stresses are predicted very accurately, transverse stresses can be obtained by the now standard approach of integrating the two-dimensional equations of equilibrium through the thickness (as in Averill and Yip, 1996a). 2.5 Summary The new theory and finite element model described show excellent promise for accurate, efficient, and convenient modeling of laminated and sandwich beams. The level of accuracy in the prediction of through-the-thickness variations of displacements, strains, and stresses can be varied by choosing the number of sublaminates used in the model. In most cases, only one sublaminate, or one finite element through-the-thickness, is needed to achieve the desired accuracy of both global deflections and local inplane stresses, even for very thick laminates. 32 1.05 1 (0”)nu ............ A ..................... A: r -A -------- 1 I ’ (0")cxact \. ‘ 0.95 : i 0.9 1 w... 1 1 ”exact 1 0.85 f 0.8 r f 0.75 3 1 ‘4 L 4A 4“ A l_ L k # L #4 A g 41 2 4 6 8 10 Number of elements in the mesh Figure 4. Predicted normalized maximum deflection and normalized maximum inplane normal stress versus number of elements in the mesh 1 _ F l 4; 5 $'—é.>fl " . 0.9 h 1 0.8 1 w I 0.7 f wexact » . 0.6 f 1 0.5 f T L I 0.4 f j 1 J 10. 100. 1000. 10000. L/h Figure 5. Normalized center deflection versus span-to-thickness ratio of a five-layer simply-supported beam subjected to a sinusoidal load 33 0.5-'f:1f"i MT 1 ' ‘ Elasticity 0'4 """""" FSDT l 0.3 . """ LZZT: 1 sublaminate 0.2 0.1 11 4s Fifi r -0.1 -0.2 -0.3 - -0.4 .05 -- r 1.051111014111051 1.5 ux10 Figure 6. Axial displacement versus normalized thickness coordinate at the end of a simply-supported beam subjected to a sinusoidal load P U! .1 Elasticity 0.4 ..... . ..... FSDT 0.3L —"‘" LZZT'Isublaminate r 0.2 0.1 a-I h II a FIN -0.1 -0.2 -0.3 -0.4 -0.5 h. V Y r .10. Figure 7. Axial stress versus normalized thickness coordinate at the midspan of a simply-supported beam subjected to a sinusoidal load 0.5 0.4 0.3 0.2 0.1 ' a-m -0.1 -0.2 .03- -0.4 -0.5 ' {x \ Elasticity “ l- \ ' \ \ ........... FSDT ; ,\ ---'- Lzznrsublamiuale fl ‘ e ———— LZZTZ'3sublaminates I — """ — LZZT: 5 sublaminates Y T 3‘1 D II & Figure 8. Transverse shear stress versus normalized thickness coordinate at 0.5 0.4 0.3: 0.2 0.1 - #IN -0.l - -0.2 -0.3 -0.4 -0.5 - . 6.1 Figure 9. x=U10 of a simply-supported beam subjected to a sinusoidal load V V *T a? fii ' v v v v v V7 Elasticity u" —-—-- Lzzrz-rsublamiuaie — - _ — LZZT: 3 sublaminates _ “““ — LZZT: 5 sublaminates A A A l A L A A l A 6.3 _6 6.4 6.5 6.6 wxlO 'II'ansverse deflection versus normalized thickness coordinate at midspan of a simply-supported beam subjected to a sinusoidal load dis orc C0115 Chapter 3 MINDLIN PLATE ELEMENT In the following discussion, we shall study the small deformation response of plates using Mindlin plate theory, a popular theory in the application of finite element methods for the bending of plate structures. This is a prologue to the more complex laminate theory to be discussed later. This chapter will illustrate the field and edge consistency concepts that are crucial in the development of robust quadrilateral plate elements without introducing the complex nature of the displacement fields of laminate theories. 3.1 Introduction The use of Mindlin plate bending theory for the development of plate finite elements is popular, as only a C0 continuity requirement is imposed on the field variables (transverse displacement, w and rotations, 6x and By). This allows the development of simple, low- order elements containing three or four nodes with only three bending degrees of freedom per node. Unfortunately, early studies showed that these simple low-order elements are ‘flexurally challenged’ (i.e., they ‘lock’) and exhibit violent stress oscillations in the thin regime. Prathap (1994), in his review paper, provided explanations to the locking behaviors exhibited by many such low order quadrilateral elements. Fried (1974) was the first to provide the necessary insights as to why the displacement approach failed in such a constrained problem, namely the vanishing of the transverse shear strains in the thin limit, 35 T: di; no To C01 CO! Par fOu 36 and he suggested the use of reduced/selective integration as a possible cure. Hughes, et al. (1977) employed this reduced/selective integration technique on the shear strain energy to produce the first ‘simple and efficient’ plate bending element that was free of locking. However, this plate element was found to have two zero energy modes that caused the element performance to deteriorate if the element was distorted and used in the general quadrilateral form. Recent studies by Tessler and Dong (1981), Tessler and Hughes (1983), Hughes and Tezduyar (1981) and Zienkiewicz and Xu (1993) tried to alleviate locking by attempting to satisfy the inconsistent kinematically derived strain field using an unequal order interpolation of field variables (e.g. quadratic for w, linear for 6,r and By ). The use of unequal order interpolation requires that elements have different degrees-of-freedom at different nodes. In the present case, the quadrilateral plate element would have eight nodes, with w defined at all eight nodes, 6,, and By defined only at the four corner nodes. To remove the four mid-side nodes, Tessler and Hughes (1983) imposed four edge constraints by requiring the tangential shear strains at the elemental boundaries to be constant along each of the four edges. This is, in fact, equivalent to satisfying the linear part of the tangential Kirchhoff shear strain constraint at each edge of the element. These four constraints are used to condense out the four mid-side nodes, resulting in an interdependency among the w field and the 0x and 9y fields, as opposed to the conventional models where the three fields have independent interpolation schemes. Zienkiewicz and Xu (1993) showed that their plate elements derived using mixed formulations do not lock. However, the plate element of Tessler and Hughes (1983) using the displacement approach needs to be supplemented by a residual energy balancing 37 technique using relaxation parameters to prevent locking in severe cases (Tessler, 1986). These relaxation parameters are obtained by tuning the element response for specific test cases. This energy balancing technique is now widely used. However, it is artificial in that it introduces an error that compensates for another error without removing the original error. Recently, in addition to the more familiar completeness and continuity requirements, Prathap and co-workers (1983, 1988, 1992, 1994) have performed extensive research on the importance of the field and edge-consistency requirements and variational correctness in the formulation to produce error-free, robust displacement-type finite elements. Prathap (1994) collectively called them the C -concepts, namely, continuity, completeness, consistency and (variational) correctness. He and his co-workers have developed robust finite elements based on the C-concepts, including four-node (Prathap and Somashekar, 1988) and nine-node (Naganarayana et al., 1992) quadrilateral plate bending elements. The assumed strain approach in which independent consistent shear strain fields, not derived directly from the interpolations for the displacements (Hughes and Tezduyar, 1981 and MacNeal, 1982), is one way to satisfy the field consistency requirement. The variational basis of the assumed strain methods was discussed by Simo and Hughes (1986). Substitute shear strains have also been used by Dvorkin and Bathe (1984); and Prathap and Somashekar (1988) to develop four-node elements as well as Huang and Hinton (1984); Jang and Pinsky (1987) in their development of eight- and nine-node Mindlin plate elements. Besides ensuring field consistency in their elements, Hughes and Tezduyar (1981) were the first to use tangential shear strain interpolations on each of the element edges to ensure edge-consistency in a general quadrilateral so that elements Tbs “he: and: 38 assembled in an arbitrary grid would not lock. Other recent robust plate bending elements include those by Crisfield (1983) (based on the classical plate theory), Bathe and Dvorkin (1986), and Hinton and Huang (1986), wherein mixed interpolations for the shear strains and tensorial transformation were used to achieve field and edge consistent elements. In the current paper, the main focus is to improve the performance of the four-node element MIN4 of Tessler and Hughes (1985) by identifying and eliminating the inconsistent terms in the shear strain fields using an approach similar to that of Prathap and Somashekar (1988). This will remove the need for ‘tunable’ relaxation parameters to prevent locking in the thin regime. After performing such modifications, the consistent shear strain fields of the new MIN4-CC element are, surprisingly, shown to be identical to those developed by Prathap and Somashekar (1988) in the element QUAD4-CC. The primary differences between the elements then appear to be the representation of consistent load vectors and the post-processing of transverse deflection, which is quadratic in the present element. The behavior of these and other elements are thoroughly illustrated by numerical examples. 3.2 Formulation The Mindlin theory for a linearly elastic isotropic plate is briefly summarized below. The displacement components are: ux(x,y,2) = -29,(x,y) uy(x. y, z) = —29y(x. y) (18) “Ax, y. z') = W(x, y) where w(x, y) and rotations, 6x(x, y) and 9y(x, y) are defined at the plate’s midsurface and the sign convention is as shown in Figure 10. All quantities are referred to a fixed Sll will fall. 39 system of rectangular, Cartesian coordinates. A general point in this system is denoted by either (x, y,z) or (x., x2, x3) , whichever is more convenient. Throughout, Latin and Greek indices range from 1 to 3 and 1 to 2, respectively, unless otherwise stated. From Eq. (18), the strain components will be: e = —z9 E = -Ze 8.. = xx x,x yy y.y «A. (19) sz = (wry—6y) 7x2 = (wax-ex) ny = ‘z(0x,y+0y,x) where a coma in the indices denotes partial differentiation. Assuming cu is negligible, and with the generalized plane stress assumptions, the stresses are given by: E on = _ 2[Z(9x,x+”9y,y)] l-u — E 9 9 on " _ 2[Z(U x,x+ ya” (20) 1—1) on = —G[z(0x,y+0y,x)] on = G(w,x-0x), 0y: = G(w,y-9y) where E and G are the elastic and shear modulus, respectively, and t) is the Poisson’s ratio. Moment resultants of the stresses 6 and (5xy are defined as: xx’ ny Mxx = Izcxxdz = -D(9x’x+uey’y) My), = JzO'yde = -D(09x’x+9y,y) (21) l-u M”, = Izoxydz = —DT(9x,y+9y,x) where M n and M yy are bending moments, M is the twisting moment, and D is the xy bending rigidity defined by .J’ 40 3 12(1 — uz) where h is the thickness of the plate. Finally, the transverse shear force resultants, Q x and Qy are defined as: Q, = [0,de = kZGh(w,x - ex) (23) Q), = joyzdz = kZGh(w,y-—9y) The factor k2 in Eq. (23) is a shear correction factor that attempts to take into account the actual non—uniform distribution of shear stress in the thickness direction. A value of k2 = 1:2/ 12 is used for all computations here. Figure 10. Sign convention for rotations and stress resultants 41 The equations of motion, the natural boundary conditions and the displacement—based finite element are obtained using the principle of virtual work: SW = BWW+ SW“, = O (24) where Win, and W are the internal and external work done, respectively. The virtual ext strain energy of the plate in indicial notation is given by: SWW = IoijfieijdV i, j = 1, 2, 3 (25) where the integration is carried out over the volume, V, and summation of repeated indices is implied. The virtual work due to external forces is: SW“, = — J(q5w+ma89a)d£2+ j (MBanB)59adI'— j Qanaawdr n r, r, (26) (M3 = 1,2 where q, the distributed transverse loading in the z-direction, and distributed moments, ma are given by: q = Ipfzdz‘i'lczzlflf/ZZ _ fwd +'.'(o I H: I ml '— IZP x Z 2 xZ h/2 XZ -h/2) (27) m2 = szfde + g(oyz|h/2 + OyzLh/z) p is the density of the plate, and If = fig. is the body force per unit mass. M 0,3 and Qa are the prescribed boundary moments and shear forces, respectively, 9 and F are the area and boundary domains of the plate, respectively. Using the decomposition rule, (M gunmen can be written as: (1145034989,I = Mnn59n+Mnsaes (28) where 56,I = Seana, 595 = 59—59,}; and 59 = seaga. dl: lh: 42 Now, let the moduli Dams be defined by: 01111: D, 01122 = ”D, 02211 = “D, D2222 = D 1—1) (29) 01212 = 01221 = D2112 = 02121 = TD and all other terms are zero. Then using Eqs. (18-23) in Eq. (25), the virtual work equation can be written as: j 050,759,, 559.1, Bdo — j kGh(w,a — cameada Q n + JkGh(w,a—9a)8w,adn = L1(59)+L2(5W) (30) n a,B,y,5 = 1,2 where the boundary and load conditions for the plate are written as two linear functionals: L,(6Q) = —Jma50a dQ- j (Mnn69n+Mnsaes) dI‘ Q ’2 (31) L2(5w) = — Iqow dQ— answ d1" (n, s: no sum) (2 r, F refers to the complete boundary. If F1 is the portion of the boundary where displacements are specified and F2 is the portion where surface tractions are specified, then I‘lul"2 = I‘ and I‘lnl‘2 = Q. The equilibrium equations are obtained as: MBovB" Q0: = ”ma (32) —Qa,a = q de “01 MOT. 43 and the following boundary conditions can be identified: Simply Clamped 3939:3314 supgloztwed-Zl Free Symmetric symmetric W = 0 W = 0 w = O Q = O Q = 0 93:0 MM=0 95:0 Mns=0 Mm=0 n: Mun: nn=0 Mun: 9n=0 Note that simply supported-1 (SS-1) is the appropriate simply-supported condition for Mindlin-type plate elements to deal with general boundary types such as curved- boundaries. Simply-supported-2 (SS-2) is based on thin-plate theory, and in cases where there is no danger of over-constraining such as in simply-supported, thin rectangular plates, SS-2 can be used to further eliminate degrees-of-freedom. 3.3 Finite Element Model The finite element model utilizes the interdependent interpolation concept of Tessler and Hughes (1983, 1985). The transverse deflection, w is initially approximated using quadratic serendipity interpolation functions, N 1., while all other degrees-of-freedom are approximated using linear Lagrange interpolation functions, Pi. This, however, will result in mid-side degrees-of-freedom associated with w (a consequence of the higher order interpolation functions used). The finite element therefore is not of a convenient form as depicted by Figure 11a. The mid-side nodes are condensed out to achieve a uniform four- node quadrilateral element. This is done by explicitly enforcing the linear part of the Kirchhoff constraints, YSZ’S = (W’s—95)”- = O (33) along each side. A brief summary is outlined below (refer to Appendix B for details). 0 9 c k (a) unconstrained element (b) constrained element 0 {w, 9x, 9),} O { W} Figure 11. Topology of the finite element model Independent bi-linear in—plane displacements and rotations and a Serendipity (eight- node) transverse deflection are initially assumed, i.e., 8 w = 2N1”; i=1 4 4 e,r = 213,9“. ey = 210,9”. (34) i=1 i=1 4 4 where Pi(§, n) and N i(§, n) are respectively bi-linear and Serendipity shape functions, and 9 By, and W, denote nodal degrees-of-freedom (refer to Figure 12). 0X ll xi’ 9: l mu Vat a: Figure 12. Quadrilateral element coordinate description 8V “'1 3.4 DUI 45 The four mid-side w degrees-of-freedom are condensed out by the use of four differential edge constraints: (w._.,.-e,,,.)|§’n=il = 0 (se 10,111) (35) where 1,, is the edge length and subscript k represents the corresponding edge as depicted in Figure 12. Each of these constraints is explicitly solved for wk + 4 and back-substituted into Eq. (34) (refer to Appendix B for details). This produces a four-node, quadratic, coupled transverse deflection: 4 w = 2 [Piwi + Naiflgi + Nnieni] (36) i: 1 where N§i and Nut” are given by: 1 I N§1='N5’ N§2=-NE,1’ N§3=-ZN7’ Ng4=—N§3 4 l 1 (37) N11] =ZN8, Nfl2=ZN6’ Nn3=—Nn2’ Nn4=_an and 1 1 N5=§(1—§2)(1—T]), N6: i(1_112)(1+§) l 1 (38) N7=§(l_§2)(l+n), N8: i(1_nz)(1_§) 91;: and 9111' in 59- (36) are rotations corresponding to the natural (local) coordinate system (see Figure 8.1). This approach effectively increases the order of the element without introducing any additional nodes or degrees-of-freedom. 3.4 Consistent Shear Strain Fields The interdependent interpolation scheme alleviates the shear locking problem but does not eliminate it totally. This stems from the fact that constrained strain fields, i.e. the 46 transverse shear strains, are still not field consistent. In a general quadrilateral, it is difficult to see the consistent form of the shear strains clearly because of the non-uniform mapping from a natural (local) coordinate system into a Cartesian (global) system. In order to see this clearly, the transverse shear strains need to be converted to the natural coordinate system. Prathap and Somashekar (1988) termed the transverse shear strains in the natural coordinate system as covariant shear strains. Examining the covariant shear strain component, Ygz’ using Eqs. (34, 37, 38): 752 = (ng- 9g) = Ppgwi + (Néi’é - 13199:: + Nni'éeni = “—BIQFGHI +6T1249113 +9114} (39) (l-Tl) 4 0+“) + {—egl—Bgz—Wl +W2}+ 4 {—9§3—9§4+w3—w4} The coefficients of the quadratic term (1 —n2) are clearly inconsistent since they contain only 6711' rotations alone and will therefore lead to spurious constraints and shear locking when thin plates are modeled. Therefore, eliminating the inconsistent terms and simplifying, a consistent shear strain, 752 is obtained. The consistent shear strain is mathematically represented as: 7g: = Ppgwi + (Navg — P096 (40) This shear strain is equivalent to interpolating the w degree-of-freedom using: W = Piwi+N§i6§i (41) instead of Eq. (36) with the interpolation strategy for the rest of the degrees-of-freedom remaining unchanged. In a similar manner, the consistent shear strain, '7,” can be derived from the dc SE 5q 511 C01 for 9011. Intel Sim Bode 47 inconsistent shear strain, 7112 and is given by: 7112 = Pi’nwi+ (Nni’n-Pi)eni (42) or with the w degree-of-freedom interpolated using: w = Piwi + Nniflni (43) This consistent field approach using the new assumed shear strain fields is similar to the selective/reduced integration concept except that it only removes the relevant inconsistent terms and therefore will not produce spurious zero energy modes. Prathap (1994) demonstrated through the use of the Hu-Washizu theorem that the kinematically derived strain field when replaced with a consistent constrained strain field is equivalent to a least- squares smoothing operation. Simo and Hughes (1986) also showed that the assumed strain method is variationally correct 3.5 Edge Consistency Even with the choice of field-consistent transverse shear strain fields, edge- consistency requirements must still be ensured to prevent locking in an element, especially for general quadrilateral cases. It is essential that the tangential shear strain, Y” on an edge be consistently matched or else there is a spurious constraint generated on the edge. This can be achieved if the transformation from the natural coordinate to the Cartesian coordinate systems are done in such a way that the edge constraints are preserved consistently. Normal jacobian transformations of two adjoining elements at their own integration points will redefine this consistency causing a mismatch of tangential shear strain at the common edge. Prathap and Somashekar (1988) suggested the use of element nodal coordinate transformations for the interpolation of the covariant—based shear strain St: N( C0] 48 fields to preserve the consistency definitions of tangential shear strain, 7“ along a common edge between two elements. The nodal covariant degrees-of-freedom are transformed to the Cartesian counterparts {7&2} = all 012 {VIZ} {VIZ} = b“ b12{Y§Z} (44) ynz “21022 sz 7y: b21b22 Y112 where [a] = [b]_I . Following Eq. (19), the Cartesian coordinate based shear strain, 7x2 using: can now be expanded as: Yx (g’n) = W’ —9 Z x x (45) = (wig... + w,.,n.,> - 9, Introducing the covariant rotations, 9: and on into Eq. (45) and using the transformation definition of Eq. (44), the Cartesian based shear strain, 7“ can be expressed in terms of the two covariant based shear strains, 7:: and 7112 ‘sz(§’ Tl) = (W,§ — 6§)§9x + (Win — en)n9x = 7&3, + 7,211., (46) = 75,sz ‘" Ynzbrz and now Eq. (46) can be modified using the definitions of the consistent transverse shear strain fields, '75: and 7112 from Eqs. (41 , 42) to give a consistent shear strain, E: 7,4611) = 7§zb11+7mb12 = (Pi,§wi + (Nip: -' Pi)9§i)b11+(P"nWi+(N l (47) mm " Pi)9ni)b12 Nodal transformation must be canied out for the nodal covariant rotations, 6:1. and 9111' to convert them back to the Cartesian coordinate definitions of the rotations, 6n- and Gyi. Thus, the shear strain, '7“ in terms of On. and fly,- can be written as: 49 4 7.45.11) = 2‘, {9...1(a,.),(N§,.§— Fab“ +(a21),(N,,,.,, —P,-)b121 i=1 (48) + eyi[(al2)i(N§i’§ " Pi)b11+(a22)i(Nni’n _ Pi)bl2] where by. are evaluated at the integration points and a l. j are computed at the nodal points. Similarly, the shear strain definition for 7y: will be: 4 7,413.11) = 2 {e,.-[(a,,),(N§,-.§ - Fab21 + (a21),\\ /J<45°\ j Loading: unit forces at free end / . . A (b) Trapezordal Theorctlcall out-of-plane: 0.4321 twist: 0.03406 / A45° é / / / J / 1 l (c) Parallelogram Figure 14. Straight cantilever beam pa tw 116.1 153 Plar. Table 3. Straight cantilever beam 52 Tip load direction MIN4-CC MIN4 QUAD4-CC QUAD4 AQR8 (a) rectangular shaped elements Out-of-plane shear 0.980 0.980 0.980 0.986 0.981 Twist 0.850 1.203 0.850 0.886 1.011 (b) trapezoidal shaped elements Out-of-plane shear 0.963 0.893 0.963 0.968 0.965 Twist 0.805 1.248 0.805 0.896 1.029 (c) parallelogram shaped elements Out-of-plane shear 0.978 0.975 0.978 0.977 0.980 Twist 0.778 1.117 0.778 0.890 1.159 All elements performed well with the regular mesh. The irregular-shaped elements which contain considerable distortions cause performance deterioration of all the elements except for AQR8 as can be seen from the results in Table 3. For all elements, both the parallelogram and trapezoidal meshes caused unsymmetrical results at the free ends with a twist load. This, however, is to be expected because of the unsymmetrical mesh. Values reported above are therefore the mean values. Curved Cantilever Beam. The curved cantilever beam as shown in Figure 15 is the next test. This will also be a test of the effect of slight irregularity. Out-of-plane shear load is applied at the free end to produce the out-of-plane deformations. The normalized out-of- plane deflection for the curved cantilever beam is reported in Table 4. 53 Geometric and Material Properties: inner radius = 4.12; outer radius = 4.32; are = 90° : thickness = 0.1; E = 1.0x107; u = 0.25; Loading: out-of-plane unit forces at tip Theoretical tip deflection: 0.5022 F ired Figure 15. Curved cantilever beam Table 4. Curved cantilever beam Tip load direction MIN4-CC MIN4 QUAD4-CC QUAD4 AQR8 Out-of-plane shear 0.944 1.023 0.944 0.951 0.956 Long Cantilever Test. This test was taken from Prathap and Somashekar (1988). Tip out-of—plane shear loads are used to produce a constant variation of shear force and a linear variation of bending moment. Different configurations are used to check for distortion sensitivity. These are as shown in Figure 16. // / 40 Geometric and Material Properties: % é \ length = 100; 0 50 ,00 /b 60 100 Width = 1.0; depth = 0.1; (a) (b) 10 E = 1.0x106; / / /l) 90 [00 A0 100 Loading: unit forces at free end (c) (4) Theoretical: “4000: (Mxx)| = 75 ; (Mu)2 = 25 Figure 16. Long cantilever beam 54 Normalized results obtained for the different configurations are identical to the four- node element of Prathap and Somashekar (1988). The regular mesh configuration produced results that agree with field-consistent 2-node linear shear flexible beam elements. Linear bending moment can also be reproduced at the centroids of each of the elements. However, in the distorted configurations (b, c and d), the elements can only produce constant bending moments. Prathap and Somashekar (1988) reported that such results indicate that the assemblage of two distorted elements acts only with the efficiency of a single linear beam element. The effect of distortion will therefore lower the efficiency of the two-element patch. Note, however, that even in extreme distortion i.e. collapsing the quadrilateral elements into triangular elements (configuration d), locking is not present. Table 5. Long cantilever beam test Configuration w (tip) MIX M xx .____=____=1 (element 1) (element 2) (a) 0.938 1.000 1.000 (b) 0.750 0.667 2.000 (C) 0.750 0.667 2.000 (d) 0.750 0.667 2.000 Rectangular Plate. A rectangular plate will be used to test the convergence characteristics of the new element under different boundary supports, aspect ratio and loading conditions (see Figure 17). The test results are detailed in Tables 6 and 7. Uniform convergence characteristics are clearly seen for all elements in all cases. In the case of a clamped plate with aspect ratio of five, models that use MIN4 converge more slowly than the other models. 55 Geometric and Material Properties: a = 2.0; b= 2.0 or 10.0; thickness = 0.0001; E = 1.7472x107 ; u = 0.3 ; Boundaries = simply supported or clamped 7 sym Mesh = N x N (on 1/4 of plate); fl—l | J Loadingzuniform load, q = 10_4 or I b j center load, P = 4.0)(10—4 Theoretical solutions: Boundary supports Aspect ratio Displacement at center of plate b/a Uniform pressure Center load Simple 1.0 4.062 1 1.60 Simple 5.0 12.97 16.96 Clamped 1.0 1.26 5.60 Clamped 5.0 2.56 7.23 Figure 17. Geometry, material properties, boundary conditions and theoretical solutions for the rectangular plate tests Table 6. Results for rectangular plate, simple supports: uniform load Normalized transverse defection at center (a) aspect ratio = 1.0 Mesh MIN4-CC MIN4 QUAD4-CC QUAD4 2 x 2 1.079 1.035 0.978 0.981 4x4 1.019 1.009 0.995 1.004 6 x 6 1.008 1.004 0.998 1.003 8 x 8 1.005 1.002 0.999 1.002 (b) aspect ratio = 5.0 Mesh MIN4-CC MINL QUAD4-CC QUAD4 2 x 2 0.800 1.098L 0.971 1.052 4 x 4 0.992 1.001 0.978 0.991 6 x 6 0.996 1.001 0.991 0.997 8 x 8 0.998 1.001 0.995 0.998 3 d For def Bur cen- Elfin €ch Very 56 Table 7. Results for rectangular plate, clamped supports: center load Normalized transverse defection at center (a) aspect ratio = 1.0 Mesh MIN4-CC MIN4 QUAD4-CC QUAD4 2 x 2 0.865 0.990 0.865 0.934 4 x 4 0.965 1.000 0.965 1.010 6 x 6 0.985 1.001 0.985 1.012 8 x 8 0.992 1.000 0.992 1.010 (b) aspect ratio = 5.0 Mesh MIN4-CC MIN4 QUAD4-CC QUAD4 2 x 2 0.318 0.634 0.318 0.519 4 x 4 0.825 0.699 0.825 0.863 6 x 6 0.911 0.795 0.911 0.940 8 x 8 0.946 0.860 0.946 0.972 To assess the robustness of the element, a square, simply-supported plate subjected to a doubly-sinusoidal load was analyzed for span-to—thickness ratios, a/h, from 10 to 105. For all cases, an 8x8 mesh was used in a quarter-model of the plate. The normalized center deflection and maximum normal stress results (with respect to the elasticity solution of Burton and Noor, 1994) are plotted in Figures 18 to 21. Stresses are evaluated at the centroid of each element, which is the optimal stress recovery point for the MIN4-CC element (see Barlow, 1976). There is no sign of numerical ill-conditioning even with extreme thinness. Predictions of center deflection using MIN4 deteriorate rapidly in the very thin regime for a plate aspect ratio of one. These results clearly show that the element is very robust for both moderately thick and thin plates. 57 1.005 , r ' fi'fifi ' f ' r v--1l v v 4 ffifi..| , , V . 7" 1 E f:— ................ ‘ .................... A .................... A l l, 0.995 j 0.99} { iv : : 0.985} «1 0.98: 5 E + MIN4-CC 3 “-975 .' ~~~A MIN4 1 : . — - -- — QUAD4-CC 3 0.97 r j ”1000. ‘ 10000 7100000 . a/h Figure 18. Normalized center deflection versus span-to-thickness ratio of a simply-supported square plate under sinusoidal loading (a/b = l) 1 V. f; .'._. '7’; ' '_" ‘ '_" 0.95 « W 0.9 . 0.85- . MIN4-CC i t 1 0.8 _ MIN4 , . QUAD4-CC . 10. 100. 1000. 10000. 100000. a / h Figure 19. Normalized center deflection versus span-to-thickness ratio of a simply-supported rectangular plate under sinusoidal loading (db = 5) 58 ] - - - 1 f f ’ 4 0'99 A .................... A .................... ‘ ..................... A4 - ........... - —————————— - ........... - 0.98 6xx 0.” - + MIN4-CC I —-I-— QUAD4-CC A I 0.95; ‘AU‘UAL - - - ...--..1 . ..--..1. 10. 100. 101/11,); 10000. 100000. a Figure 20. Normalized maximum normal stress versus span-to-thickness ratio of a simply-supported square plate under sinusoidal loading (a/b=1) 1 ~ w - + . . . - - 1 0.975 _ ____f____ 0.95 - 0.925 . an 0.9. 1 0375’ + MIN4-CC 111.... 10. T “”100. ‘ ‘ 11000 0‘7710000.‘ ‘ 100000 Figure 21. Normalized maximum normal stress versus span-to-thickness ratio of a simply-supported rectangular plate under sinusoidal loading (a/b = 5) 611 .11 59 Circular Plate. A circular plate with clamped edges is commonly used to analyze an element’s behavior in a moderately thick plate regime. Of particular interest is a solution due to a center point load. According to Reissner’s theory, a logarithmic singularity exists under the load. Figure 22 shows the coarse mesh of a symmetric quadrant used to analyze the problem. y 41 ‘ Geometric and Material Properties: radius, a = 5.0; 1 thickness, h = 0.1 or 2.0; Boundaries = clamped gym MCSh = 12, E = 1.0x106; u = 0.3 Loadingzuniform pressure, q = 1 or center load, P = 4 sym Theoretical: 1 Uniform _ i 2_ 2 2 _4_ 2 2_ 2 w- 64D[(a r) +l-‘Uh (a r)] Center load V). 8D 1 r] 2 "- kGha a P 2 2 2 a W — m[(a -r )—2r In;- Figure 22. Clamped circular plate (quadrant) Figures 23-26 show the non-dimensionalized deflection along an edge of the clamped circular plate. Results of center deflection show that the MIN4-CC element is capable of reproducing accurate results for a clamped circular plate under uniform or center loads. Note that the deflection curves obtained using the MIN4—CC and MIN4 elements are quadratic in nature whereas only linear variations are predicted by the QUAD4-CC element. However, this is only a result of the post-processing of deflection in the case of MIN4-CC, and this could be done for the other elements as well. -0.2 -0.4 1 £1 -0.6 P Reissner Theory -0.8 : + MIN4-CC [ A - MIN4 . —--- — QUAD4-CC -1 - . _ 0 1 2 3 4 5 Figure 23. Normalized deflection (W = wl61tD/(Pa4)) along line of symmetry (y=0) of a thin (a/h=50) clamped circular plate under center load 0'7 T I -0.5' l -1- . w I -151- ‘ . ——- ReissnerTheory 1 : + MIN4-CC , -2 ' WA MIN4 1 p ' _..._ QUAD4-CC . 0 1 2 3 4 5 Figure 24. Normalized deflection (W =wl61tD/(Pa4)) along line of symmetry (y=0) of a thick (a/h=2.5) clamped circular plate under center load 0 7 V 4 —0.2 f - -0.4 r - _ l . W 1 . -0.6 - j I 1 _0 8 ’_ Reissner Theory .‘ ' » MIN4-CC ‘ » “A” MIN4 4 .1 L - - -- - QUAD4-CC .‘ 0 1 2 3 4 5 Figure 25. Normalized deflection (W=w64D/(qa2)) along line of symmetry (y=0) of a thin (a/h=50) clamped circular plate under uniform load 0 . r . l 1 '005 P " , . , l W .1 . . l ,- ' .A f P /-/ .' T r 7, ' _. " — Reissner Theory 1 -1.5- ,. —+— MIN4-CC : Z ,-’-’ A" A MIN4 j .‘ ............. ' —.--— QUAD4-CC 4 0 1 2 3 4 5 x Figure 26. Normalized deflection (W = w640/(qa2)) along line of symmetry (y=0) of a thick (a/h=2.5) clamped circular plate under uniform load VT 10 [l 001 33.; the 1 CODd reSul Const 62 Morely’s Acute Skew Plate. The last test is to assess the behavior of the new element under skew distortions (see Figure 27). Of particular interest is the acute 30° skew plate simply-supported on all edges. An exact solution for the uniformly loaded condition has been provided by Morley (1963) which indicated a singularity at the vicinity of the obtuse vertices, with moments of opposite sign. Geometric and Material Properties: 3:100; thickness = 0.1; YA 55-] or ss-Z angle, B = 30° or 60° Boundaries: simply-supported (SS-1 or SS-2) ss-l 33.; Mesh = Nx N / /a E = 1.0x10"; 0 = 0.3 S - r 55-2 é Loading: uniform load, q = 1 J] Theoretical: S a qa‘ ‘3 w = 0.408—D—X10 (B = 30°) 4 —3 w = 2.505%X10 (B = 60°) Figure 27. Morley’s acute skew plate All solutions except for those obtained using MIN4 converged very slowly because of the singularity and showed a much stiffer solution for the acute 30° angle case compared to the solution of Morley (1963). Results obtained with MIN4-CC are comparable to the nine-node element of Naganarayana et al. (1992). The choice of using a combination of SS-l and 88-2 boundary conditions or purely SS-2 boundary conditions has an effect on the results for the 30° angle case. Solutions using the combination of SS-l and SS-2 boundary conditions yield a stiffer result with a coarse mesh but both choices of boundary conditions converged to about the same solution when the mesh was refined. The stiffer result using the combination type boundary conditions could be traced to an over- constraining of thin plate elements caused by SS-2. This phenomenon was also seen by 3.7 Cdg, dish appr 63 Rossow (1977) who reported that displacement based finite element models underestimated the deflection at the center of the plate by almost 20 percent even with a very refined mesh. Since the acute skewed plate is a numerically difficult problem, the accuracy of the results can be considered to be fairly good. Table 8. Normalized center deflection: uniform load (a) B = 30° (i) Combination of 88-1 & SS-2 b.c. (ii) Purely SS-l b.c. Mesh MIN4-CC MIN4 QUAD4-CC MIN4-CC MIN4 QUAD4-CC 2 x 2 0.811 1.637 0.487 1.037 2.030 0.622 4 x 4 0.756 1.118 0.699 0.956 1.187 0.877 8 x 8 0.774 1.021 0.762 0.856 1.052 0.841 16 x 16 0.818 0.999 0.815 0.841 1.018 0.838 32 x 32 0.863 0.995 0.862 0.876 1.008 0.875 (b) B = 60° (i) Combination of 88-1 & SS-2 b.c. (ii) Purely SS-l b.c. Mesh MIN4-CC MIN4 QUAD4-CC MIN4-CC MIN4 QUAD4-CC 2x2 1.185 1.362 0.711 1.359 1.619 0.815 4 x 4 1.014 1.074 0.925 1.041 1.11 1 0.948 8 x 8 0.987 1.028 0.965 0.992 1.038 0.970 16 x 16 0.990 1.009 0.985 0.991 1.016 0.985 32 x 32 0.986 1.003 0.985 0.986 1.007 0.985 3.7 Summary The MIN4-CC element contains consistent shear strain fields as well as consistent edge tangential shear strains, obviating locking even when the element is severely distorted. Even though an interdependent interpolation scheme was used for the element approximation of deflection, the present consistent strain fields were shown to be equivalent to those developed by Prathap and Somashekar (1988). Thus, the main difi 016 than CC Ihec the c 64 differences between these two elements are in the consistent load vectors and in the predicted transverse displacement distributions. Numerous numerical tests demonstrated that the MIN4-CC element is both robust and accurate. The assumed strain field approach that was successfully utilized to derive the MIN4- CC element will next be used to develop a robust LZZB element using our new laminate theory. It will ensure the satisfaction of both field and edge consistency requirements in the element and prevent element locking in the thin regime. Chapter 4 LAYERWISE PLATE ELEMENT - THICK PANELS The field and edge consistency concepts have been introduced in the previous chapter. It is now easier to introduce the new laminate theory for the bending of plate structures. In this chapter we will solely concentrate on the formulation of the new zig-zag sub-laminate theory for plates and the associated finite element model with direct application to thick sandwich panels. The extension of the finite element model to thin plate bending structures will be discussed in the next chapter which deals with the locking phenomenon. 4.1 Introduction The layerwise plate element derived using the new laminate theory could take the form of a four-node plate element. However, similar to the case of the layerwise beam element discussed in Chapter 2, it is more advantageous to cast the element in the form of an eight- node brick. This topology allows the laminate thickness to be conveniently subdivided and modeled by multiple finite elements (representing sublaminates). It is thus possible to increase the accuracy of the finite element model, as needed, to capture through-the- thickness gradients and transverse (interlaminar) stresses. It is also possible to simulate delaminations using the redundant node concept. Both the theory and the finite element model have been developed in two forms. The first form, called here LZZ3, is based on a high-order zig-zag theory and will be discussed 65 ht di [01 Var EEK 66 in detail in this chapter. The finite element model LZZ3 contains seven degrees-of- freedom per node: u, v, w, 0x, 0y, 1,, 1y . The first five degrees-of-freedom are the usual ones used in plate/shell finite elements. The remaining two degrees-of-freedom, Ix, 13y , are the transverse shear tractions (or interlaminar shear stresses, as the case may be). If only one sublaminate (element) is used through the thickness of a laminate, then these transverse shear degrees-of-freedom can be eliminated at the element level, leaving only the traditional five engineering degrees-of-freedom. The second form of the model, called here LZZl, is based on a first-order zig-zag theory (see Cho and Averill, 1996) and contains the usual five engineering degrees-of-freedom only. This model will not be discussed here. 4.2 Theory Formulation In the present theory, the laminate is composed of N T perfectly bonded layers stacked together in the thickness direction. The thickness and material stiffness properties may vary arbitrarily from layer to layer. These layers can be modeled as M sublaminates, with each sublaminate containing N m layers, where m is the sublaminate number (see Figure 28). Mathematically, this is represented as: NT= sz (50) de 101 111116 3.8311 “Shel: 67 z is a local (sublaminate) thickness coordinate Z,z with its origin at the bottom of the sublaminate T Z is the global (laminate) thickness coordinate ZM . ; v ° 1 ' ////?}ZW///////‘/ Z ‘ 4 z 2 m "2‘” Sublaminate 2 : Wm z2(2) Z Zr (2) I T z0(2) ----- Sn blaminate—I— - — — — Figure 28. Schematic of sublaminate and layer divisions In order to facilitate the development of the theory, all expressions in the following derivation pertain to the mth sublaminate. The sublaminate number designation is omitted for brevity. The constitutive relations for the kth layer of the mth sublaminate with respect to the laminate coordinate axes are given by: T k k k k- to ‘(k) C(11) C(12) C(13) 0 0 C(16) Fe ‘00 "x C(k) Can 0 0 Cut) x" Oyy 22 23 26 8yy (10 (k) 0.22 = C33 0 0 C36 822 (51) Ty: Cff,’ cf]: 0 7y. 13 xz Sym C151? 0 7.2 .1”; C(k) .7101. 66 The 13 anisotropic material constants in Eq. (51) can be expressed in terms of 9 independent constants. The layer strain-displacement relations using small strain assumptions are: = u, 8 = 12,}, 822 = W,Z 8 x Y)’ XX (52) M = w,y+v,z 112 = w,x+u,z ny — v,x+u,y where a coma in the indices denotes partial differentiation. 50! 81's 1:111 6115 dim 68 In each sublaminate, an independent displacement field is assumed in which the through-the-thickness variation of in-plane displacements is described by a cubic polynomial in the local thickness coordinate with a piecewise linear (or zig-zag) function superposed upon it. The transverse deflection is assumed to vary quadratically with the thickness coordinate. The assumed displacement fields are initially assumed in the following form: 3 n—l It. . ui")(x,y,z) = z z uk+ z (z—z,)§i k=0 i=1 3 n-l n k. . k=0 i=1 1...): 0(1 —11+w:(%)+%(ill -1—1 where the subscripts b and I refer to the bottom and top surfaces, respectively, of the sublaminate. h is the total thickness of sublaminate m. All quantities are referred to a fixed system of rectangular, Cartesian coordinates. A general point in this system is denoted by either (x,y,z) or (x1, x2, x3), whichever is more convenient. The use of the non-conforming bubble function assigned to field variable uz is to ensure the removal of the Poisson’s ratio stiffening effect whenever the full three- dimensional constitutive relations are used. Poisson’s ratio stiffening effect is the phenomenon that is observed in continuum elements when only one layer of linear elements is used through the depth to model a region of flexure. This stiffening effect originates from the bending energy terms and emerges from the inability of the approximation functions to interpolate the transverse normal strain, 8 linearly through 22’ the depth of the plate (see Prathap, 1985). th sh 61: am “’he 69 Since the sole purpose of the bubble function is to eliminate the Poisson’s ratio stiffening effect, terms involving W3 in the computation of transverse shear strain relations are ignored. This simplification will eventually lead to W3 being a nodeless variable (degree-of-freedom) and thus C-1 continuous in the resulting finite element model. It is then condensed out at the element level. It is possible to eliminate the degrees-of-freedom E,- and 1"], in Eq. (53) by enforcing the conditions of transverse shear stress continuity at each interface. The conditions of shear stress continuity at the jth interface, where z = 2]. , are: (1') _ TU“) 1:(J') _ TU“) (54) sz - xz yz - yz Making use of Eq. (53) along with the infinitesimal strain-displacement and linear elastic constitutive relations, Eq. (54) can be solved to give closed-form expressions for E, and 01' in terms of the remaining degrees of freedom. . . Bwb .. A 8w, aw, ‘11 = “11(V1'l'a—y )+azi"2+a3i"3+a4t(§; ‘a—y ) 8w, Bwb) .. aw, .. .. +b“(ul+0—y )+b21“2+b3i“3+b4i(a_x _a—x - = . +_ + . + . + ._ __ t 611 V1 8y Cztvz C3tl’3 C41 3y ay (55) . aw, .. . 8w, aw, +dh-(ul-1-5; )+d2iu2+d3iU3+d4i('a—x- —§; ) where Thes Each trans. 11min 70 (k) k—l k-l apk = ak[fp + apq]+bkL cm] q=l =1 (k) k—l k—l bpk - bk[fp + 2 dpq]+ak[q2 bP‘l] q=l =1 ( ) k-l k-l cpk = Ck fp + Z anJ-l-dkL cpq] q= =1 (56) (k) k-l k-l dpk = d,[fp + Z dpq)+c,£z hm] q=l :1 (k) 2 f(lk) _ 1 f3 — 32,, (k) (k) zk f2 - 22p f4 - Z and ,. l k 1 k k 1 k a. = — <50 Ak+1 . 1 (1+1) k k 1 k (6.. C§.’- Cfo’ ’Cfts’1—1 Q, R H Ak+1 (k+l) (k+1) (k+1)2 Ak+1=C44 C55 ‘(C45 ) These conditions therefore reduce the number of degrees of freedom by 2(Nm — 1) for each sublaminate, leaving us with ten degrees of freedom per sublaminate. Additional simplification of Eq. (53) can be achieved through satisfaction of the transverse shear traction boundary conditions at the top and bottom surfaces of the laminate. These conditions are: “’41 C01} 71 To) _ T T1N.) 1: xz z = 7-0 xb xz z = h xt (58) 1(1) _ T (N...) _ T yz z=zo yb yz z=h y! where 111,, In, tyb, 1),, are the applied shear tractions (or interlaminar shear stresses, as the case may be) at the bottom and top surfaces of the sublaminate. The displacement field is cast in its final form by introducing the surface variables: l A A 0:“51) =“o ”0:“;0 - zo Z0 (59) (N..) _ (N..) Eqs. (58, 59) can be solved either analytically or numerically (see Appendix C) so that the displacement field for each sublaminate can be expressed in terms of the operative degrees of freedom (all functions of the inplane coordinates, x and y, only): 12 (60) ab, vb, wb, 0x1), Byb, txb, tyb, u,, v,, wt, 0“, 0 xi, 12y, yt’ where the subscripts b and t refer to the bottom and top surfaces of the sublaminate, respectively, and Bwb aw, 9x0 = '37 9x! = 3; (61) awb aw, 0y, = 3'; 9,. = '87 The replacement of the partial derivatives of the deflection with the “new” rotations was done to allow the assumed displacement degrees-of-freedom to remain C0 continuous. The constraints in Eq. (61) will have to be enforced explicitly by a penalty formulation during the development of the finite element model. The displacement field in Eq. (53) for the mth sublaminate can now be represented in 0V. 72 terms of through-thickness shape functions as: ( ) 4(k) k ux" = ua (x, y)(k)(z)dz 1:1 zt-l NM (MU)? = 2 (:1 rxy(x, y,z)‘1’(k )(z):z) i=1 z,_ 01 = 1,2; k =1,2,4,..., 74 N, z, (Q90, = [I]: rxz(x, y, z)Ma(z)dz] i= 1 1 N" z,- (Qy)a = 2 £1 Tyz(x, y,z)M a(z)dz] i= 1 —1 N, z,- (Qz)a = [ I Ozz(x, y, z)Ma, Z(z)dz] i: 1 7-1-1 Nun 2‘ (67) (Qz) = 2(1 “a” y, Z)M3 z(z)dz] 121: (R,)ff’= [I 1,,(x y.z).,aagf> ds — j [(Qn)a8wa + 9,803] ds :2 32 01,0: 1,2; k = 1,2,4, ...,7 or refers to either top or bottom surface of the plate where q, the distributed transverse loading in the z-direction, and the distributed “generalized” moments, 21 and ma are given by: 901 = Ipsza dz+[ozzMa]g q = Ipsz3 dz (m1)? = jpbxoff’dznt 42"”): xza 01 = 1,2;k = 1,2,4, ...,7 (70) h 0 k k k (m2)( ’ = J‘pbwa,’ dz+[tyz\rf,’] a M 113’ M 23, Q, and Q" are the prescribed “generalized” boundary moments and shear forces and they are: (k) k k (M,,)f,’ = (My)? (114.02" = (M,,),, (k) _ (k) (k) _ (k) (M21), - (Myx). (M22)... - (Myy)a (71) (Q11)... = JrnzMadz Q" = IrnzM3dz or =1,2;k =1,2,4,...,7 s1 refers to the portion of the top and bottom boundary contour s of the sublaminate where displacements are specified and s2 is the portion where surface tractions are specified. Making use of Eqs. (64, 65, 69), the virtual work equation can be written as: To the gen “he 76 0 - 8(Mxx)0k) MM”)? 8(Mxy)g()+a(Myx)0f) - I - 8x + By + 3y 8): -(R,)f,f"—(R,)ff’184f,f" (72) (9,), 8(0),, . . A — ax + a; -(Qz)a]8wa+(Qz)6w3]dA+L(ug‘))+L(wa) a =1,2;k = l,2,4,...,7 where L028”), L(wa) are the boundary terms that result from integration by parts and 12:: ) is defined by Eq. (53). The boundary and load conditions for the mth sublaminate can be written as two linear functionals: . k A k k A k L(ug‘)) = — 1(1):, +m2);)5u;) dA— 1(MIB+M25); 013611;) ds A 32 L(w.) = - 114.5%. +0503] 44 - 11(2),». + (“2,6031 ds (73) A 52 01 =1,2;k = l,2,4,...,7 and the governing equations in indicial notation are given by: swa: (-Qp,0+Q3).. = q. ..(k). 2 (k) _ 2 (k) 5%. - 20:1(M10’1‘R0h " ‘20=1(’"B)a (74) 8W3: ’(Q3)a = ac 01,117 =1,2; k =1,2,4,...,7 To understand the boundary conditions associated with this higher-order laminate theory, the degrees-of-freedom in the displacement field of Eq. (62) can be recast into a more general form: u = Hanna) ”(0,1,(0 x (I: :10 :0 :0 I: 1.2.3; <1 = 1’2 (75) “y = anawna'l'fisa‘ysa where the new set of degrees-of-freedom will be represented as: 77 ( ( ( S S S -1 _1_2_3 a=1,2 (76) al.): (al. 1 u}. l. u}. l). = (14.9.1.1. and subscripts s and n refer to the tangential and normal directions, respectively, on the boundary surface. The shape functions, (DU) (1)“) ‘l’ig and WE}; are readily obtained 311’ nor from the original shape functions (see Appendix D) using transformation equations for Cartesian vectors, e.g. u luxux+l u n fly )7 (77) u, = —lnyux+lnxuy where In = cos(n, x) = dy/ds and l"y = cos(n, y) = —dx/ds are direction cosines of the outward normal of the boundary. Substituting Eq.(75) into the boundary and load terms associated with the iig‘) degrees-of-freedom of Eq. (73) will yield: L(ag<>) = —j[(ms)g>8rtw+(mn)g>5rtw] dA A (78) + 11(M.)§P517.o + (It/[01005171111] ds 52 where the new generalized moments are defined as: ’ l I l 11 (“1’51” = 1919910 b.1111 dz + [1.913. + 0.111 . (l) (79) 1 (M 5)}? = IKOxxnx + txyny)<1>f3;+ (“trxynJr + oyyny)‘1’ 1: 1,2,3; 01 =1,2; [3: s,n 1301] 612 de 935 685 COI [rat 78 Now, the following sets of boundary conditions can be identified: Simply Simply . Clamped supportwl 3:553:35“ Free Symmetric“ um = 0 uml = 0 um = 0 (May) = 0 (1145):!” = 0 u... = 0 (May) = 0 (M39) = 0 (M33) = 0 u... = 0 wa=0 wa=0 wa=0 Qna=0 Qna=0 9.. = 0 (M053) = 0 9.. = 0 (Mag?) = o (Mpg?) = 0 0.. = 0 (Mpg?) = 0 (M,)f,2> = 0 (Mpg?) = 0 9.. = 0 I... = 0 (M053) = o t... = 0 (M03) = 0 (M03?) = o I... = 0 (M553) = 0 (M333) = 0 (44353) = 0 I... = 0 *Symmetric conditions are only applicable to laminates with orthotropic materials that are stacked with their principal material directions at either 0° or 90° to the coordinate axes (cross-ply lami- nates are a special case of these specially orthotropic laminates). The essential boundary conditions are similar to the plate/shell theory formulation with the exception of the transverse “shear traction” (or interlaminar shear stresses, as the case may be) degrees-of-freedom. A closer look at the boundary conditions will, however, seem to yield an ambiguity in the boundary conditions associated with the “shear traction” degrees-of-freedom, 1: m and Tna- These degrees-of—freedom appear explicitly in both the essential and natural boundary conditions. In actual fact, there is no ambiguity. The essential boundary conditions associated with the “shear traction” degrees-of-freedom for the boundary types must always be complied with as dictated by the variational approach. In cases where the shear traction conditions are known, (e. g., zero shear traction is commonly encountered), it is therefore desirable (but not necessary) to specify the “shear traction” degrees-of-freedom to correspond to this shear traction condition. However, in ant de: “11 79 the event when the applied shear traction load is non-zero, it is imperative that the shear traction must be included as an applied load using either Eqs. (70, 71) or Eq. (79). Also, note that simply supported-1 (SS-l) is the appropriate simply-supported condition for the plate elements with both general laminate and curved boundary types, and are very similar to those of the Mindlin plate element. For cross-ply laminates or laminates with orthotropic materials that are stacked with their principal material directions at either 0° or 90° to the coordinate axes, SS-2 can be used to further eliminate degrees-of-freedom since 15 and In are directly related to the shear strains through the kinematic assumptions. While the finite element geometry could take the form of a four-node plate element with fourteen degrees-of-freedom per node, it is advantageous to use the topology of an eight-node brick element as shown in Figure 29a. This topology allows the substructuring of the laminate into sublaminates for mesh refinement in the thickness direction. This will increase accuracy and also permit the simulation of delamination using the redundant node concept. Furthermore, using the eight-node brick element topology allows this new 3-D ‘structural’ element to be coupled with conventional 3-D continuum elements, although compatibility would not be strictly satisfied in this case. The finite element model utilizes the interdependent interpolation concept of Tessler and Dong (1981), Tessler and Hughes (1983, 1985) and Tessler (1990). The transverse deflection, wb, w, is initially approximated using quadratic interpolation functions, N j, while all other degrees-of-freedom are approximated using linear Lagrange interpolation functions, P j. This unequal interpolation scheme is a consequence of trying to satisfy the Kirchhoff constraints in a consistent manner. This, however, will result in mid-side degrees-of-freedom associated with wb, w, (a consequence of the higher order 4.4 final 31X 1 511101 80 interpolation functions used). The finite element therefore is not of a convenient form as depicted by Figure 29b. The mid-side nodes are condensed out to achieve a uniform eight- node brick element. This is done by explicitly enforcing the linear part of the constraints of Eq. (53) along each side in a similar fashion as Tessler and Hughes (1983) (see Appendix D for details). This approach effectively increases the order of the element without introducing any additional nodes or degrees of freedom and alleviates the shear locking problem so an exact order of numerical integration may be used. a 8-node brick element §% ( ) (b ) unconstrained element (“1’ Vt) w,, ext’ eyt’ In: Ty!) (“1r vh’ W1» 9x0, eyb’ "xh’ Iyb) (c) constrained LZZJ element mid-side nodes condensed out 12:: Figure 29. Topology of the LZZ3 finite element model 7 (“1’ vl’ ”'1’ 0”, 0”, I“, tyt)7 2 (“b' Vb’ Wb’ 910’ eyb’ 1:rtb’tyb)2 4.4 Numerical Studies Several numerical experiments dealing with anisotropic plates are presented. Spectral analysis of the LZZ3’s stiffness matrix revealed six zero eigenvalues associated with the six rigid body modes. No spurious zero energy modes are present. This is to be expected since the stiffness matrix is integrated using a 3x3 Gauss quadrature rule “in the plane”. Numerical results are presented for bending of three simply-supported laminated ele C0 501 61'; €161 Ihre 81 square panels subjected to a double-sinusoidally varying transverse load (see Figure 30). The laminates are very thick, having a span-to-thickness ratio of four. Numerical examples such as these are an excellent test of a model’s ability to capture local effects. Physically, these examples are similar in many ways to local models of the region directly beneath an impact load. P—————————.———___..._.__’ double-sinusoidal """""""""""" transverse loading Note: Load is distributed over entire top surface Figure 30. Schematic of loading and boundary conditions for example problems Numerical results are obtained using a quarter-model with discretization of 8x8 elements in the plane of the panel and one or more elements through the thickness. Comparisons are made between predictions of the current model and an exact elasticity solution (Burton and Noor, 1994). The stress distributions across the thickness are evaluated at the center of the elements, since the center can be treated as the optimal location for stress recovery (Barlow, 1976). A Random Five-layer Panel. The first example is a random five—layer laminate. The lamination scheme is given in Table 9 and the material properties used are listed in Table 10. Three levels of thickness discretization were used. The first model employed a single element through-the-thickness (a sublaminate containing all five layers), the second used three elements through-the-thickness (one for the bottom layer (i.e. layer one), one for the 82 next to bottom layer and the third, for the remaining three layers) and the last model used five elements (one for each physical layer) through-the-thickness of the entire panel. Table 9. Lamination scheme for random five-layer panel Layer No. Material Thickness 1 l 0. 10 2 2 0.25 3 3 0. 15 4 1 0.20 5 3 0.30 Table 10. Material properties for random five-layer and sandwich panel Matl. 1 Matl. 2 Matl. 3 Core E” 1.0e6 33.0e6 25.0e6 5.0e4 522 25.0e6 21 .0e6 1.0e6 1.5e5 E 1.0e6 21.0e6 1.0e6 5.0e4 33 0 0.01 0.25 0.25 0.01 12 0 0.25 0.25 0.25 0.15 23 0 0.25 0.25 0.25 0.15 13 G 0.5e6 8.0e6 0.5e6 2.17e4 12 G 0.5e6 4.0e6 0.2e6 4.20e4 23 G 0.2e6 8.0e6 0.5e6 2.17e4 13 A Sandwich Panel. The second example is a sandwich panel with five composite face sheets on the top and bottom of a core material. The material properties are already listed in Table 10. The lamination scheme is given in Table 11. Note that the five-layer face sheet laminate is by itself a rather challenging laminate to analyze due to the variations of 83 material properties and layer thicknesses. Three levels of thickness discretization were used. The first model employed a single element through-the-thickness of the entire sandwich panel, the second used two elements through-the-thickness - subdividing the laminate equally, and the other used three elements through-the-thickness - one for the bottom face sheets (a sublaminate containing five layers), one for the core, and one for the top face sheets. Table 11. Lamination scheme for sandwich panel Layer No. Material Thickness l 1 0.010 2 2 0.025 3 3 0.015 4 1 0.020 5 3 0.030 6 Core 0.800 7 3 0.030 8 1 0.020 9 3 0.015 10 2 0.025 1 1 1 0.010 A Panel from TACOM ’s Composite Armored Vehicle. The third example utilized a lamination scheme that is similar to that present in the US Army TACOM’s Composite Armored Vehicle (CAV). The laminate can be divided into three sections, as described below: (1) Inner Shell. Four plies of S-2 Glass/Phenolic Fabric with stacking sequence { [902/02] }with thickness of 0.01 for the four plies stack. Thirty-seven plies of S- 2 Glass/855340 Epoxy Tows with stacking sequence {[0/ 90], [45/—45/0/9O]4 , [45/0/45], [90/0/—45/45]4 } with thickness of 0.021. 001 lay wit m0 The Tat emph 84 (2) Armor Core. One layer of EPDM rubber (thickness = 0.0625) and one layer of ceramic tile with inserts (thickness = 0.7). (3) Outer Shell. Twelve plies of S-2 Glass/855340 Epoxy Fabric. Stacking sequence is [0/90/45/-45/0/90]s with thickness of 0.01. The laminate consists of 55 layers and is nearly two inches thick. To allow comparisons between predictions of the current model and the exact elasticity solution, all layers are oriented at either 0. or 900 from the reference axis by replacing all 45° plies with 90° plies and the —45° plies with 0° plies. This is not a restriction of the current model, however, which is capable of modeling completely general lamination schemes. The material properties are listed in Table 12. Table 12. Material properties for TACOM’s composite armored vehicle (CAV) panel S-2 Glass/ S-2 Glass/ S-2 Glass/ EPDM Phenolic 8553-40 8553-40 Ceramic Fabric Tow Fabric R“b°°’ _ 1—__ = L i: E” 3.0e6 6.2e6 3.0e6 3.0e3 5.0e6 £22 3.0e6 1.0e6 3.0e6 3.0e3 5.0e6 £33 1.2e6 1.0e6 1.1e6 3.0e3 1.25e5 1) 0.13 0.29 0.13 0.45 0.15 12 u 0.18 0.37 0.18 0.45 0.15 23 u 0.18 0.29 0.18 0.45 0.15 13 612 1.0e6 0.3e6 1.0e6 1.0e3 2.5e6 G 4.6e5 0.3e6 3.9e5 1.0e3 8.5e3 23 013 4.6e5 0.3e6 3.9e5 1.0e3 8.5e3 Once again, three levels of thickness discretization were used. The first model employed a single element through the thickness of the entire panel, the second used four 85 elements through-the-thickness with each sublaminate representing the different sections - one for the inner shell, two for the armor core, and one for the outer shell. The third model used three elements t11rough-the~thickness with the two sublaminates in the armor core combined and the rest similar to the second model. The convergence characteristics of the different LZZ3 models were studied with different mesh density. The center deflection at the mid-surface was used for comparison even though the variation of the transverse deflection through-the-thickness for a thick laminate can be substantial. This is readily seen in Figures 32, 38 and 46. Three discretizations were used for all cases. The cases with one and three elements through-the- thickness (sublaminates) are as described above. The two elements case is obtained simply by subdividing the panel thickness into two equal parts. The convergence characteristics for all the cases are depicted in Table 13. Table 13 showed that for all cases even with a coarse mesh of 4x4, the results are converged with the exception of the one-sublaminate model in the CAV panel case. For the CAV panel, results clearly demonstrated that thickness refinement through the use of more sublaminates is needed for laminates with such complicated layups. A note of caution is that more thorough comparison between the different (sublaminate) models must, in general, take into account through-the-thickness variation of the displacements and stresses for thick panels. 86 Table 13. Results for square plate, simple supports: sinusoidal load Normalized transverse deflection at center (a) random five-layer panel Mesh l sublaminate 2 sublaminates 3 sublaminates 2 x 2 1.070 1.072 1.073 4x4 1.008 1.008 1.011 6 x 6 0.997 0.997 0.999 8 x 8 0.978 0.993 0.998 (b) sandwich panel Mesh 1 sublaminate 2 sublaminates 3 sublaminates 2x2 1.131 1.120 1.133 4 x 4 1.054 1.033 1.056 6x6 1.040 1.014 1.042 8 x 8 1.035 1.007 1.037 (c) CAV panel Mesh 1 sublaminate 2 sublaminates 3 sublaminates 2 x 2 0.863 1.052 1.057 4 x 4 0.822 0.999 1.003 6 x 6 0.815 0.989 0.993 8 x 8 0.812 0.986 0.989 Results of the analyses for all three examples are shown in Figures 31 to 51 for the 8x8 mesh discretization, where predictions by the LZZ3 models of the through-thickness variations of inplane displacement, ux, transverse deflection, u inplane stress 0'”, Z, inplane shear stress ’1: transverse shear stress In and transverse normal strain, ezz are xy’ compared to the variations predicted by three-dimensional elasticity. Even for these extremely thick laminates, it can be seen that good predictions of ux, uZ , on and 1,0, are obtained using only one element (sublaminate) through-the-thickness, with slight improvements in these predictions when the thickness is discretized using three or more elements. This is to be expected since the deflection, uz , and the inplane stresses, CH and 87 ’1' can be captured adequately by zig-zag theories which are equivalent to our LZZ3 xy ’ model when one sublaminate is used. In Figures 34, 41 and 49, the transverse shear stress is calculated using the constitutive relations with no shear correction factors. In each example, more than one element through-the-thickness is needed to capture the variation of sz adequately. The deviation from the exact solution is more drastic in the results of the sandwich and CAV panels (Figures 41 and 49) where these panels have material properties that are drastically different in adjacent plies. The predictions of sz could be improved for the case when one element is used through-the-thickness by recovering the stress using the three- dimensional equations of equilibrium. From Figures 36, 44 and 51, plots of the normal transverse strain ezz are shown. Results of an is not expected to be good for models using small number of elements through-the-thickness because of the simple approximation of the assumed displacement field, uz. However, one can clearly see that the normal transverse strain, ezz will approach the elasticity solution as more elements are used in the through-the-thickness direction. This is especially seen in Figure 36 for the random five-layer panel when one sublaminate is used to represent each layer. In general, however, accurate normal stress, 0' has to be 22 ’ recovered using the three-dimensional equation of equilibrium if the number of elements used in the through-thickness direction is kept low. 4.5 Summary The current LZZB model shows excellent promise for efficient and accurate analysis of thick laminated composites and sandwich panels. This zig-zag approach has the following desirable pr0perties: (i) interlaminar transverse shear stress continuity conditions are 88 satisfied; (ii) transverse shear tractions at the top and bottom surfaces of the plate are satisfied exactly; (iii) transverse normal tractions at the top and bottom will tend to the prescribed load value with an increase in the number of sub-layers; (iv) a small and fixed number of degrees-of-freedom per sublaminate are needed to accurately describe the kinematic behavior of complex sublaminate regions; (v) traditional engineering degrees- of-freedom (displacements and rotations) are used; and (vi) the eight-node brick topology permits the use of adaptive techniques for through-thickness discretization. The casting of the technical theory in terms of surface quantities also allow the model to represent the applied loads better. In most plate theories, surface tractions are applied at the reference surface which usually corresponds to the mid-surface. This is acceptable for thick plates but not for thin plates. For a thick plate as observed by the numerical experiments, the through-the-thickness distribution of in-plane displacements, stresses as well as the transverse stresses are non-symmetric with respect to the mid-surface. Results showed that the current LZZ3 model is much more efficient and potentially more accurate than 3D continuum-based models for analysis of laminated composites and sandwich panels. The LZZ3 model is also computationally competitive with traditional low-order plate elements. 89 Elasticity ........... LZZ3: 1 sublamuzate ————— LZZ3: 3 sublaminates 4 a» — — — — LZZ3: 5 sublaminates 0.1 " // 7: = 4 + 2 ~. 7. ° "- -o.1 L \" -o.2 - \ J -0.3 i '004 \ .05LAAA- ,Anm ‘4 A A; ‘ .. . 44444 -o.4 -o.2 o 0.2 0.4 0-6 -6 uxXIO Figure 31. Axial displacement versus normalized thickness coordinate of a 0.5 0.4 0.3- 0.2 0.1 3‘“! -0.1 - -0.2 - -0.3 -0.4 .05_:;. simply-supported square plate (random S-layer) subjected to sinusoidal load Y Elasticity LZZ3: 1 sublaminate -‘--- LZZ3:3subIaminates —--- LZZ3:5subIaminates A +4 A A A A l A A A L 3.2 3.3 A ‘ 3.4 ‘ A 3.5! -6 uzx10 Figure 32. Transverse deflection versus normalized thickness coordinate of a simply-supported square plate (random 5-layer) subjected to sinusoidal load 0.5 i T i I i ' l 7 f I I T Elasticity 0.4 . ---------- LZZ3: 1 sublaminate i ----- LZZ3: 3 sublaminates “-3 ' ——-— LZZ3:5sublamiuates 0.2 0.1 ' .a. = 4 h 4 5 o -0.1 - -o.2 - -0.3 ~ « -0.4 -o.5 - Figure 33. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (random 5-1ayer) subjected to sinusoidal load 05 _ vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv . 0.4 . a .~ * _ = 4 0.3 . h w 0.2 0.1 ' 5 0 Elasticity “ """"" LZZ3: 1 sublammate' -0.1 D ,' -' ————— LZZ3: 3 sublammate' s ‘ ,,,,,, :73"? --—- LZZ3:5sublammate° s .002 ' , ' 1 .03 _ ......................... -0.4 ............. \ -o.s - ................................. ' « ' 1 e2 - 1 “0.8 .006 '004 '002 0 1x2 Figure 34. Transverse shear stress versus normalized thickness coordinate of a simply-supported square plate (random S-Iayer) subjected to sinusoidal load a-m 0.5 0.4 - 0.3 0.2 0.1 -o.1 ~ -o.2 - -o.3 ~ -0.4 -o.s~ 91 r vvvvvvvvv T Elasticity LZZ3: I sublaminate — -- - - L223: 3 sublaminates - - — - LZZ3: 5sublaminates AAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAA A Figure 35. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (random 5-layer) subjected to sinusoidal load i J 0.5 - i .\ W 0.4 - a = 4 .- .‘. h .' \ 0.3 . '\ 0.2 " 0.1 - . Z 0 l \. 7. ,x _0.l . I .\ -o.2 - ‘\ Elasticity .0'3 P \ """""" LZZ3: I sublaminate . _0 4 -~-‘- LZZ3:3sublamiuates ' .- \ — — —- - LZZ3: 5 sublaminates -0.5 - A - . . - A . l . - . . - . . A 0 0.2 0.4 0.6 0.8 -6 23x10 Figure 36. 'Il'ansverse normal strain versus normalized thickness coordinate of a simply-supported square plate (random S-Iayer) subjected to sinusoidal load 92 0.5»' ' ' fi' - ~ - - . . - , - a f 04 xx 0'3“ Elasticity 02. ...... ”23:13“! 0 I ' —-—-— 1.23:2:ch . l 01. ——-— LZZ3:3sublamiuates _/ + . ./ z i.- 0 -o.1» a -0.2 7'_4 -0.3r -0.4 «0.5L 4 . A -l uxx10-6 Figure 37. Axial displacement versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.4 0.3 . Elasticity """""" LZZ3: I sublammate' 0.2 ' ----- LZZ3: 2 sublaminates 01L --—- LZZ3:3sublamiuates a-IN V 41.1 -o.2 -o.3 -o.4 .05 1 I' . A A. 1 _A. A A L —6 uzx10 Figure 38. 'h'ansverse deflection versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5 VVVVVVVVVVVVVV 1 ”K 0.4 0.3 P 0.2 " g = h 4 0.1 ' E 0 h -0.1 L Elasticity i _01 _ """""" LZZ3: 1 sublaminate -'- - - LZZ3:2sublamiuates '03,, “-—" LZZ3:3SflbWS 43.4 M'fi fl -0.5-.".-:-.-...4---7 ................ 1 -30 -20 -10 0 10 20 30 0 Figure 39. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5 a 0.49 0.465 5. 0.45 h — Elasticity 0.43 ........... ”23.13“”. I —-1 - . — ~ - LZZ3: 2 sublaminates - — — - LZZ3: 3 sublaminates 0.4-30 4+0A3‘310‘AA‘20AA‘A304 0' Figure 40. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load - top face sheet 0-5 ’ '..'__'_'_'_L—'—'—'--'-'-3--'——'-'—'-—'_—'_—'—_L—'._—I_l.L..Z_'._'..i._‘..'_;_'-'_' ‘ ' '- « 'EI:___ .r' 004 T1. I 0.3 . i J 0.2 ' 1 fl = 4 h 001 . Ill .4 II' 2 : ‘ it 0 I . W. .0 2 _ .’ """"" LZZ3: I sublammate’ i ' 1|; ----- LZZ3:2sublaminates .03 . 2 - - — - LZZ3: 3 sublammates .004 ”_‘_ -—- l\"-.. -o.5 a - . - - -‘.“.‘:Z‘if‘.“f.‘::f‘f'—f—T—"f:‘i‘fj:“Tf. f- . . . - . -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 1' xz Figure 41. 'Ii‘ansverse shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidalload 0.5i"ififivv""";'VYV'YYTYTV'f‘ \ 0.4 03L 0.2 * g: 4 i 0.1' 5 0 " L '04 Elasticity .0 2 b """"" LZZ3: 1 sublammate' & ’ -‘—-- LZZ3:2sublaminates 4.3L —--- LZZ3:3sublaminates 1 -0.4 \ MA -0.5 A n A . ‘. l A L - . . . A, -10 -5 0 5 10 Txy Figure 42. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 95 0.49 0.465 E 0045 h Elasticity 0,43 """""" LZZ3: 1 sublaminate ——< ----- LZZ3:2sublaminates —— -- LZZ3: 3sublalninates 0.4 AAAAAA - A A 4 A - . - . . ...... - n A - -10 -5 0 5 10 I” Figure 43. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load - top face sheet 0.5~'-.:—s' ' i I 0.4 i ‘C"-~ \ 03- “"'~--.\ 0.2- ‘ 0.1. E. 0 h -0.1- i —— Elasticity -0.2 - """""" LZZ3: 1 sublaminate - - - - - LZZ3: 2 sublaminates -0.3* ———- LZZ3:3sublamiuates . -0.4 . -o.5--LA.-.----.3.--.3 0 5 10 15 .5 czleo Figure 44. 'Ii-ansverse normal strain versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 96 0 5 _ ............ f- - - ................ ‘ Elasticity """""" .LZZ3:1subhunbuue --- - - LZZ3:3sublaminates — — - - LZZ3: 4sublaminates J 5 h 4&5 ' . -- - - . -454 I .-ss-l,- ............... kg‘ -7.5 -5 -2.5 0 2.5 5 7.5 —6 u x10 Figure 45. Axial displacement versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 0.5: '7?" 7' 'Tfi' ' 'fijfi' - 0.4r ' 0.3- 4 a 0.2- 5 -4 0.1P z - 0- h .001. \ . :' | 1 ”'2' 5' l 9 Elasticity .03_ ' ,' 1 """""" LZZ3:Isublanunate ;’ I , _._-_ LZZ3:3sublaminatcs -0.4- g l,’ ,- — —- — — LZZ3: 4 sublaminates -o.s--':---.A---.--------/-{- -- ---------A- 42.5 45 47.5 50 52.5 55 57.5 —6 uzx10 Figure 46. 'Ii'ansverse deflection versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 97 0.5 ' f ' ' 7 Y 004 F .3: 'l 0.3 ~ Elasticity a _ 4 """"" LZZ3: 1 sublanuuate' 0.2 ' f, " - - - - - LZZ3: 3 sublaminates * 0 l I} - — — - LZZ3: 4 sublaminates 4t1 4L2- 4L3» 4L4» 4L5--. ll-- 3-IN -20. Figure 47. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 0 f . - 1 ~ - . -0J. fl = . h 4 -o.2 - 5. h -0.3 * Elasticity I ~~~~~~~~~~ LZZ3: 1 sublaminate -' —---- LZZ3:3sublaminates 4’" ' I” J — - — — L223.- 4sublaminates .20 A ‘ A .10 A A ‘ A 0 L A A A 10 ‘ A A #720 4* 0 Figure 48. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load - outer shell 0.4 - v - 1 v v v v V r v v r V fiv f v V r 1 r v—v— ff? ...... 0.5; i . e . . / v A I ’ q . . . . 0.3 " —— Elasticity 2 4 ----------- LZZ3:1 sublammate' 0.2 , . h --—-— LZZ3:3sublanunates 01_ ---- LZZ3:4sublaminates FIN -o.1 - -o.2 .03 -o.4 - -o.s - Figure 49. 'h'ansverse shear stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 05_ . . . r as - - . . . . . . . . . . . . . . 0.4) T'- " 0.3r a _=4 0.2~ h 0.1. J 5 o- . 'oel’ Elasticity -0.2 i """ ' """ LZZ3: I sublaminate _03 b ----- LZZ3: 3 sublaminates — — — - LZZ3: 4sublaminates -o.4- J -2 0 2 4 6 1:, Figure 50. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 99 0.5 - 7 f? f a T 0.4 - 0.3 P a at 0.2 - 7s = 4 0.1 . E 0 / . II -0.1 - . _03 . Elasticity : """""" LZZ3: 1 sublaminate .03 i ----- LZZ3: 3sublanu'nates 1 -- —- LZZ3:4subla.minates -0.4 *- . -o.5 » 3 - - - - . - - - - . - a - A . l s 0 20 40 60 —6 5:1le0 Figure 51. Transverse normal strain versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load Chapter 5 LAYERWISE PLATE ELEMENT - THIN PANELS The layerwise plate element developed in Chapter 4, like many conventional displacement-based C0 -plate finite elements, is plagued by the shear locking phenomenon when utilized in the thin regime. This chapter will primarily deal with the locking problem. The assumed strain field approach described in Chapter 3 for a Mindlin C0 -plate element will be adapted for use with the layerwise C0 -plate element to eradicate locking in the thin regime. 5.1 Introduction In the previous chapter, a new technical theory and associated finite element model were developed for the analysis of thick laminated composite and sandwich panels. This theory, like most refined theories, is primarily intended for thick plate applications where the effects of the transverse shear and normal stresses cannot be ignored (i.e., classical laminate theory is no longer accurate). Although the refined model developed from the theory may be primarily intended for thick plate situations, it is critical that this model exhibit robustness even for plates in the thin regime. In addition, the eight-node topology of the new element permits refinement in the thickness direction. This through-the- thickness refinement using several elements can also lead to very large span-to-thickness element ratios regardless of thick or thin plate geometries. This chapter discusses the 100 101 techniques used to make the finite element model of the new sublaminate zig-zag theory robust for both thick and thin plate applications. This chapter will focus on the critical problems and techniques to resolve shear locking, Poisson’s ratio stiffening and stress oscillations that usually plague three- dimensional displacement-based finite elements developed using general higher-order displacement fields. The causes of the shear locking phenomenon in thick plate elements have always been either the presence of an inconsistent transverse shear strain field (see - Somashekar et al., 1987) or a lack of edge-consistency of the tangential transverse shear strains (Prathap and Somashekar, 1988; Mohan et al., 1994). The use of reduced and selective integration techniques (Zienkiewicz et al., 1971; Pawsey and Clough, 1971; Hughes et al., 1977, Averill and Reddy, 1992), hybrid and mixed methods (Lee and Pian, 1978; N car and Andersen, 1977; and Spilker and Munir, 1980) and modified shear strain methods (Hughes and Tezduyar, 1981; Hinton and Huang, 1986; and Crisfield, 1984) are but some of the many ways used to alleviate or remove the shear locking phenomenon encountered by plate finite elements in the thin regime. A detailed discussion of the various techniques was reviewed by Prathap (1994). Using the field and edge-consistency concepts proposed by Prathap (1994), Yip et al. (1996) developed a robust and accurate four-node C0 Mindlin-type plate element for the homogeneous and isotropic case. The fundamental concepts used to develop this four-node plate element will be utilized to derive a field- and edge-consistent eight-node brick element from the “inconsistent” element described in Chapter 4. 102 5.2 Formulation In the previous chapter, the assumed displacement fields in a sublaminate were represented in terms of through-thickness shape functions as follows: . k k u,(x. y. z) = uf.’(x.y>£2z(ug2a - a...) + (95,... - Gamma. + M.) + calms... - 1,1,1} +11 + 1116424113.. - at...) (2) (3) (101) ”91:30: — 054a)(<1>§a’ z + Ma) + (1):“, 2(Tg3a — 1:40,) }] + (iglhcbg; 2(“g1a + “@200 + (11224030, + 03:20,) + «>82, it... + it...) + Mam... - w...» + (1.1—“2mg gum + 11540,“ cg; 4953,, + 9:400 + «>233. .112... + it...) + Maw... — w...» On examining Eq. (101), the coefficients associated with the quadratic term (1 — 112) will clearly lead to an inconsistent shear strain field, 7§z since they contain only 6111‘ rotational degrees-of-freedom. This will cause spurious constraints and shear locking in the thin regime where 7&1 tends to zero. The inconsistent terms came from the 108 interdependent interpolation scheme for degreescof-freedom wa. Leaving out the inconsistent terms, a consistent ' can then be written in a concise form, £2 - - (I) 1&2 = P .4219ng + (’5':th N§i’§9§ia)Ma i= 1,...,4; l=1,2,3; a=1,2 (102) The inconsistent terms came from the interdependent interpolation scheme for degrees-of- freedom wa . The consistent transverse shear strain, E2 of Eq. (102) can be derived using the following interpolation for Wu: Wu = Piwia'l'Ngl-egia l: 1, ...,4, on = 1,2 (103) instead of Eq. (83) with the interpolation strategy for the rest of the degrees-of-freedom remaining unchanged. In a similar manner, the consistent shear strain, 3112 can be derived from the inconsistent shear strain, 7712 , identified with 'sz , and is given by: - _ - (1) I'll - Piufili)a¢lla-Z+ (Pi’nwia+Nni’nenia)Ma (104) with the wa degrees-of-freedom, now interpolated using: ya = Piwia+Nni0nia i: 1, ...,4, a = 1,2 (105) The other two consistent transverse shear strains, igz and :Ynz obtained from the inconsistent covariant based shear strains, 7&2 and y," identified with sz of Eq. (99) will be of a similar form to Eqs. (102, 104) and they are: A _ (l) _ Ptugiqufia, Z + (Ppgwia 'l' Ngi,§e§ia)Ma ~< W N l _ — 1 (I) 106 i112 - Piuilib‘llnmz + (Pi’llwia + Nfli’flenia)Ma ( ) i= 1,...,4; l=1,2,3; 01:1,2 Strictly speaking this set of new shear strain fields, {3:2, inz’ 35,2, :Ynz} is still not field 109 consistent. They are only field consistent with respect to the in-plane coordinates. However, to make them totally field consistent (which may not be necessary at all as reported by Robbins et al., 1991) is simple. Any term associated with the thickness coordinate, 2 that is of an order higher than one in the set of shape functions, {(1)22 z, <14]: 2, W122“, ‘10:; 2} must be removed. But this would severely compromise the accuracy of the model in applications to thick plates. As will be demonstrated through numerical results, such a strict enforcement of field consistency is not necessary for the present model. This consistent field approach using the new assumed shear strain fields is very similar to the selective/reduced integration concept except that it only removes the inconsistent terms. Straightforward selective/reduced integration will not work as it removes a lot more terms than necessary in the strain fields and introduces spurious zero energy modes in the process. 5.5 Consistent Penalty Constraints In the new laminate theory and finite element formulation, four constraints were introduced to produce a C0 element. The constraints correspond to the top and bottom rotations and are given as follows: Cm = (Wu—ex)“ Cw = (my—9y)Cl a = 1,2 (107) These constraints are enforced by the penalty method in the finite element formulation. Note that the form of Eq. (107) is very similar to the form of the shear strain fields obtained from Mindlin plate theory in Chapter 3. These penalty constraints therefore need to be checked for consistency just like the above transverse shear strain fields. Using the 110 relation from Eq. (86), the constraints in Eq. (107) can be recast into their covariant counterparts as: can: = “11Cxa+012Cya 18 C (0) not = “21cm "' “22Cya Following the above procedures, the covariant based penalty components, Cga , will be examined: C501 = (W’§ — 9901 (1 - 2) = +{-9nla+en2a'en3a+en4a} (109) (1- +11) 1 4’ 4n){-9§1a _ 9€201 - ”’10: + W201} + (_T—{_e§3a — 9&4“ + w3a_ w4a} Again, the coefficients of the quadratic term (1 —n2) are inconsistent as they are functions of the rotational degrees-of-freedom 0,"- alone and therefore must be eliminated to prevent any spurious constraints. Eliminating the inconsistent terms and simplifying, the consistent penalty component, Cga can be obtained by interpolating the wa degrees- of-freedom using exactly the same form as in Eq. (102) above instead of Eq. (83) with the interpolation strategy for the rotational degrees-of-freedom remaining unchanged. Cg“ can then be written in a concise form, 6.3,, = Ppgwia + (N§,,§—P,.)9§ia (110) In a similar manner, the penalty constraint component, Cm can be derived from the mconmstent component, cm! as: Cnz = mem+ (NW, —P,.)en,.a (111) with the wa degrees-of—freedom interpolated using Eq. (105) above. These derived fields are totally field consistent since the thickness coordinate. 2, does not explicitly appear in 111 any of the equations. 5.6 Consistent 'li'ansverse Normal Strain Field A problem that is normally encountered with finite elements that utilize the three- dimensional constitutive relations is the rapid deterioration of the computed transverse normal stress as the span-to-thickness ratio increases. This problem is especially severe in modeling of isotropic plates under flexure. Predicted bending deflections are typically stiffer by about 10%, depending on the value of Poisson’s ratio, when one element is used through-the-thickness. This stiffening effect, however, is not apparent if (1) reduced (plane-stress) constitutive relations are used or (ii) when two or more elements are used through-the-thickness or (iii) for laminated plates. Prathap (1994) termed the effect as Poisson’s ratio stiffening. This phenomenon was attributed to an interpolation inconsistency that prevents the finite element model from simulating a state of zero transverse normal stress in the presence of general non-zero bending strains (see Robins et al., 1991). Addition of the bubble function using a nodeless variable W3 in the displacement field, uz , as seen in Eq. (80) eliminates the stiffening effect. This, however, does not ensure the field—consistency in the definition of the transverse normal strain, ezz. Besides the normal Kirchhoff constraints, i.e. the vanishing of the transverse shear strains, there is another conStraint related to the vanishing of the transverse normal stress, 82s = "(C313xx "' C328yy + C367xyVC33 (112) In the thin limit, the transverse normal strain must be consistent with the combined Poisson effects from the dominating in-plane strains. Notice that 8zz exhibits a higher order in-plane interpolation (quadratic) and a lower order transverse interpolation (linear) than 8 8),), and ny' In general, these interpolation inconsistencies will prevent Eq. xx’ 112 (112) from being satisfied exactly and a deleterious numerical constraint will be imposed. Therefore, these inconsistencies in the transverse normal strain field must be eliminated. Again, using the same argument as in the transverse shear strains, only the in-plane interpolation inconsistencies need to be removed as the contribution to the strain energy by terms involving the thickness coordinate z is minimal. The in—plane normal strains, an and e the in-plane shear strain, 'ny as well as the Y)” transverse normal strain, Eu are expanded in terms of the covariant based degrees-of- freedom in a similar fashion as in Eq. (101). Omitting all the tedious algebraic processes, the consistent transverse normal strain, ézz , is found to be: ézz = {(1-§)[W1+ W4 + (6111 —6714)]a + (l + §)[w2 + 1123 + (0,]2 - 9713)]11 (113) +(1+11)[w3+w4+(0&4—0g3)]a}Ma,z/8+W3M3,z or = 1,2 or in a concise form: ézz(§, Tl) = (Piwia + 11750:“, + 1171".0§im)Mm,z + W3M3,z (114) i = 1,...,4,0l =1,2 where the interpolation functions are: P] = (2—é—n1/8; fig. =(1-n1/8; 1V... = (1—§)/8; F2 = (2+E3-n1/8; figs = -(1-n)/8; 7% = (1+§)/3; (115) P3 = (2+§+n)/8; ITI§3 = -(1+n)/8; Kim = -(1+§)/8; 134 = (2-§+n)/8; 175,4 = (1 +n)/8; 171,4 = —(1-§)/8; This is different from what is normally obtained using a l-point or reduced quadrature rule to evaluate the original transverse normal strain, 82:: 113 8ZZ = uZ’Z (Piwia+N§19§ia+Nnienia)Ma’z+W3M39z (116) i=1,...,4, or =1,2 Again, the “consistent” strain field, in is only in-plane interpolation consistent and not transverse interpolation consistent (linear for E cubic for the other three strains, 8 e Z: ’ xx, yy and ny ). It can be made totally field consistent even with respect to the thickness coordinate, 2 by removing all terms associated with quadratic 2 terms or higher in the contributing a y and ny strain fields when evaluating the strain energy contribution xx’ 8y of the transverse normal strain, ézz. For the same reasons given previously for the transverse shear energy terms, this last step is not necessary as the contribution of the transverse normal strain field to the overall strain energy is very minimal in the thin regime. 5.7 Edge Consistency Even with the choice of field-consistent transverse shear strain fields, edge- consistency requirements must still be ensured to prevent locking in an element, especially for general quadrilateral cases. The consistent matching of the tangential shear strain, 752 , on any common inter-element edge is crucial. Any deviation will give rise to spurious constraints on the edges when the elements are non-rectangular. Normal jacobian transformations of two adjoining elements at their own integration points will redefine this consistency causing a mismatch of tangential shear strain at the common edge. Prathap and Somashekar (1988) showed that using element nodal coordinate transformations for the interpolation of the covariant-based shear strain fields and, in this case, the covariant- based penalty constraint fields as well, will preserve the consistency definitions of 114 tangential shear strain, 752 , as well as the penalty constraint fields, C along a common sa’ edge between two elements. So for the inconsistent Cartesian based transverse shear strain, 7x2, the set of inconsistent transverse shear strain fields, {7&2’7112} of Eq. (93) will now each be replaced with their own consistent definitions from Eqs. (103, 104). 7xz(§’ Tl) = Y§2b11+7mb12 .. (1) = (Piué21¢ga, z + (Pi’éwia + Néi’gegia)Ma)b” (117) - (l) + (PiugliLCDna’ z + (Pi’nwia + Nni’nenia)Ma)b12 i = 1,...,4; l=1,2,3; a =1,2 Following the procedure of Prathap and Somashekar (1988) (see also Yip et al. (1996)), nodal transformation must be carried out to convert the set of nodal covariant degrees-of- freedom, {limp rim-a} (see Eq. (92)) back to their Cartesian coordinate definitions, 17115,) of Eq. (81). Using Eqs. (86, 117), the Cartesian based shear strain, 7x2, after the nodal transformation will become: 4 _ 1 1 sz(§’ T1) = Z {“ia[(all)ipi¢(§02, zbll + (“2011095102, zb12] i=1 1 1 +vioj(an)ip,.‘§c:,1;ll +(a22).P.<1>( ’ 1912] l 1 1111.2 3 3 +Txia[(all)ipi¢(§ozzbll +(021)1P‘q’( ) 1712] l ' Tia. Z (3) (3) (118) + Tyiaual2)ipi(b§(1,zbll + (a22)ipi¢na,zb12] (2) (2) + exia[(all)i(Piga- Z + N§,.§Ma)b1, + (a22)1.(1>,.<1>,mM + Nanawlz] Where by. are evaluated at the integration points and al.j are computed at the nodal points. Similarly, 115 4 7W0; 11) = 2 {uia[(all)ipiq’éla),zb21 + (a21)iPi‘P$|l(:,zb22] i=1 1 1 +Via[(012).-P1‘P(§oz,zb21+(022)'P"¥( ) [’22] l 1 1111.17. (3) (3) + "xia[(a11).-P iwéa. zb21 4‘ (“21)1P1‘Pna. zb22] (119) (3) (3 + TyimI(“12)ilpiwtm. zb21 + (“22);P1‘Pnoin22] 2 2 + 6...I.—b221 (2) (2) + Oyia[(al2)i(Pi‘P§a, Z + Ngjnga)b21 + (a22)i(Pl‘PT|a, Z ‘l" Nni,nMa)b22] + Wialpi,§b21+ Pi’anZIMa} = 1, 2 This is repeated for the Cartesian coordinate based penalty constraints, C m which are written as: Cxa(§’ Tl) = (W’x—ex)a (120) = ((W’§ 9x+w9nn9x)—ex)a a = 1, 2 Introducing the covariant rotations, 9g and 9n into Eq. (120) and using the transformation definition of Eq. (86), the Cartesian based penalty constraints, Cm can now be expressed simply in terms of the two covariant based penalty constraints, Cs,“ and Cna Cxa(§’n) (w9g-e§)a§9x+(wvn—en)an9x C§a§9x+ Cnan’x (121) = Cgabll + Cnablz The above equation is then replaced with the consistent definitions of the penalty constraint components, Cg“ and CW from Eqs. (110, 111) to give a field consistent C m : an“; ll)= Egab11+ Enabn =(Pi'§wia + (Ngi’§_Pi)0§ia)bll + (Pi’nwia+ (Nni’n -Pi)0nia)b12 (122) i = 1,...,4; a =1,2 116 Once again, nodal transformation must be carried out for the nodal covariant rotations, 0b- and 9111' to convert them back to the Cartesian coordinate definitions of the rotations, 0 xi and 0”- . Thus, the penalty constraint component, Cm expressed in terms of 0“. and 0),,- can be written explicitly as: 4 Exa(§in) = 2 {exia[(all)i(N§i’§_Pi)bll + (021)1(Nni’n ‘Pi)b12] " =‘ (123) + eyia[(012)i(N§i’§ —Pi)bll + (“22).‘(Nni’n ‘ Pi)b12] +wialpi’fibll +Pi,nb12]} a = 1,2 The penalty constraint definitions for Cm are of a similar form to the above equation and are given by: 4 Eyrfigt 11) = z {9x)a[(al 1 )i(N§1'1§ - Pi)b21+(021)i(NTli’Tl — Pi)b22] i: l (124) + 9,,a[(a,2),(N§,,§ — P,.)b21 + (a22)i(Nm-,n — P0522] + wia[Pi,§b21+Pi,nb22]} 01 = 1, 2 Likewise, the consistent normal strain, ézz after transforming the covariant rotations into the Cartesian rotations will be given by: 4 £2,205.11) = 2 {9,.,1(a1,),IT/§,+ (anti/“,1 i= 1 + 9y1a[(012),-1V§,- + (022)1Nni] wiaPi }M (125) (1’2 + W3M3’z The consistent definitions of the Cartesian shear strain fields, 7x, and 7),, , normal strain field, ézz as well as the penalty constraints, Em and CW derived in this manner will eliminate locking problems due to spurious constraints and as such an exact order of integration may be used to evaluate the shear strain energy. 117 5.8 Numerical Studies First, a spectral analysis of the new, field consistent layerwise finite element model was performed to check for spurious deformation modes. An unconstrained element yields six zero eigenvalues corresponding to the six independent rigid body modes (three translations and three rotations). The new finite element was subjected to patch tests proposed by MacNeal and Harder (1985) to determine its suitability for arbitrarily distorted element shapes. Though the new finite element has an eight-node brick topology (it could also take the form of a four-node plate element except that it will lose the ability for stacking in the thickness direction; see Chapter 4 for discussion), it is based on plate kinematics and as such, it is more appropriate to use the patch tests designed for plate elements. Membrane and Kirchhofi Patch test. This is an accuracy test of the element to reproduce the pure stretching and flexural modes. The patch is shown in Figure 53. Elements are of arbitrary shape patched together to form a rectangular exterior boundary. As such, boundary conditions corresponding to constant bending curvatures are easy to apply. The applied displacement boundary conditions and the theoretical solution are also shown in Figure 53. The patch is modeled with a single layer of elements through the thickness. The present element passed both patch tests. It is able to reproduce both the in-plane strains and stresses exactly for the membrane patch test as well as the constant bending moments and the surface stresses exactly for the Kirchhoff patch tests for this arbitrarily shaped patch. Circular Plate. A circular plate will be used to check the performance of the LZZB 118 model for a non-uniform mesh. This problem was discussed in Section 3.6. A circular plate with clamped edges for two load and span-to-thickness ratios will be analyzed. Figure 54 shows the mesh of a symmetric quadrant used to analyze the problem. Two through-thickness discretizations are used: one and two sublaminates. Their results will be compared with the exact solutions based on Reissner-Mindlin theory (see Hughes et al., 1977). YT b I Node x y 4 3 1 0.04 0.02 a 2 0.18 0.03 2 3 0.16 0.08 ’ x 4 0.08 0.08 —> (a) Membrane Plate Patch Test Boundary conditions: u = 10_3(x+ y/2) v = 10'3(y+x/2) Theoretical solution: —3 e”: cyy=1xy =10 , on[ = o”, = 1333. ‘ny = 400 (b) Bending Plate Patch Test Boundary conditions: w = -10-3(x2+xy+y2) 0 X —1o’3(x/2 + y) e, —1o‘3(x + y/2) Theoretical solution: Bending moments per unit length: Mn = M” = 1.111x10'7, Mxy = 3.333x10'8 Surface stresses: on = a” = 20.667, Oxy = 10.200 Geometric and material properties: a = 0.12, b = 0.24, thickness = 0.001, E = 1x10‘5 , u = 0.25 Figure 53. Patch test for plates 119 y A Geometric and Material Properties: 1 radius, a = 5.0; thickness, II = 0.1 or 2.0; Boundaries = clamped Mesh =12; E = 1.0x106; W". 1) = 0.3 Loading:uniforrn pressure, q = 1 or center load, P = 4 1: Theoretical: ; Uniform sym 2 4 = _q—[ 2_ 2 _ 2 2_ 2] l w 640“ r)+l_uh(a r) a Centerload w = -—P—[(a2-r2)-2rzlne— 80 ln-C] 161“) r kGha2 a 513 where D = _— 12(1—v2) Figure 54. Clamped circular plate (quadrant) Figures 55 and 56 show the non-dimensionalized deflection (normalized to the thin plate solution) along an edge of the clamped circular plate for the LZZ3 models under the two different load conditions as well as span-to-thickness ratios. For the thin clamped plate under center or uniform loads, both models are capable of reproducing accurate results. For a very thick plate of span-to-thickness ratio 2.5, results correlated well with the exact solution away from the predicted logarithmic singularity in the point load case but not very well (error of 10%) in the uniform load case. The finite element solutions also showed variation of the deflection in the through-thickness direction for thick plate cases. Through-thickness deflection can differ by as much as 20% for the center point load case and 4% for the uniform load case (see Figures 55 and 56 for the deflection of the top and bottom surfaces of the plate). In view of the fact that we are trying to model a very thick plate, results obtained by both models are deemed to be reasonable. Also these results reinforce the results of the earlier patch tests in that the LZZ3 models can perform well for non-uniform mesh. 120 Reissner Theory 1 sublaminate ‘ 2 sublaminates j Figure 55. Normalized deflection (W = w161tD/(Pa4)) along line of symmetry (y=0) of a clamped circular plate under center load Reissner Theory + 1 sublaminate . “HA“ 2sublaminates ‘ . tap (loading) surface j 0 l 2 3 4 5 x Figure 56. Normalized deflection (171 = w640/(qa2)) along line of symmetry (y=0) of a clamped circular plate under uniform load 121 Numerical results are again presented for the bending of two simply-supported laminated square panels (of Chapter 4) and an isotropic plate subjected to a double- sinusoidally varying transverse load. The material property of the isotropic panel is the same as in the patch tests. The laminated panels chosen for the tests are the sandwich and TACOM’s composite armored vehicle (CAV) panels. The detailed lamination sequence, geometric and material properties for both panels are listed in Section 4.4. This time, the laminates are thin panels having a span-to-thickness ratio of ten and one thousand. The latter case will be an excellent test of a model’s ability to accurately simulate bending in very thin laminates. Numerical results are obtained using a quarter-model with various discretizations of the plane of the panel and either one or more elements through-the-thickness. Comparisons are made between predictions of the current models, a finite element model based on first-order shear deformation (Mindlin) theory and an exact elasticity solution (Burton and Noor, 1994). The Mindlin-type element selected is a nine-node Lagrangian element with selectively reduced integration and appropriate shear correction factors, (see Whitney, 1973; Chatterjee and Kulkami, 1979). The simply-supported, square plate under double-sinusoidally varying transverse load will be used to test the convergence characteristics of LZZ3 models for different mesh densities and the two span-to-thickness ratios. The center deflection at the mid-surface was used as a comparison measure. Three thickness discretizations were used for all cases. For the laminated panel, the cases with one and three elements (sublaminates) through- the-thickness are as described in Chapter 4. The two elements case is obtained simply by subdividing the panel into two equal parts. The results of all the cases are shown in Table 122 14, where predicted deflections are normalized by the exact solution. Table 14. Results for square plate, simple supports: sinusoidal load Normalized transverse deflection at center Normalized transverse deflection at center (a) isotropic panel (i) span-to—thickness ratio = 10 (ii) span-to-thickness ratio = 1000 Mesh 1 sublaminate 2 sublaminates 3 sublaminates l sublaminate 2 sublaminates 3 sublaminates 2 x 2 1.071 1.073 1.071 1.061 1.067 1.067 4x4 1.022 1.021 1.018 1.015 1.016 1.017 6x6 1.013 1.011 1.008 1.007 1.007 1.007 8x8 1.010 1.008 1.005 1.004 1.004 1.004 (b) sandwich panel (i) span-to-thickness ratio = 10 (ii) span-to-thickness ratio = 1000 Mesh 1 sublaminate 2 sublaminates 3 sublaminates 1 sublaminate 2 sublaminates 3 sublaminates 2 x 2 1.095 1.094 1.095 1.058 1.065 1.059 4 x 4 1.026 1.024 1.026 1.015 1.015 1.012 6x6 1.014 1.011 1.014 1.006 1.007 1.007 8x8 1.009 1.006 1.009 1.004 1.004 1.004 (c) CAV panel (1) span-to-thickness ratio = 10 (ii) span-to-thickness ratio = 1000 Mesh 1 sublaminate 2 sublaminates 3 sublaminates 1 sublaminate 2 sublaminates 3 sublaminates 2 x 2 1.052 1.072 1.079 1.060 1.060 1.060 4 x 4 0.993 1.012 1.020 1.014 1.014 1.014 6 x 6 0.975 0.999 1.009 1.006 1.006 1.006 8 x 8 0.971 0.997 1.005 1.002 1.003 1.003 In all cases, converged results are obtained even for a coarse mesh of 4x4. Note that for the span-to-thickness ratio of ten which is moderately thick, more thorough comparison between the different (sublaminate) models must, in general, take into account the through-the-thickness variation of the deflection (see for example, Figure 61). Plots of the normalized center deflection for span-to-thickness ratios from four to 104 are shown in Figures 57 to 59 for two of the LZZ3 models as well as the nine-node Lagrangian FSDT element. The first LZZ3 model selected employed a single element through-the-thickness of the entire panel, and the other used three elements through-the- 123 thickness. The results clearly demonstrated the superiority of the LZZ3 models as compared to finite elements based on first—order shear deformation theories. This is especially obvious in the moderately thick TACOM’s CAV panel which has a very complicated layup scheme. The LZZ3 models exhibit no shear locking or Poisson’s ratio stiffening phenomenon when applied to plates in the thin regime. However, when the element span-to-thickness ratio exceeds 103 , there is a need to increase the computational precision (e.g. double precision to quad. precision) due to numerical ill-conditioning of the stiffness terms. This is especially true in modeling thin plates with finite elements possessing full three-dimensional capability. There is, therefore a limit imposed on the span-to—thickness ratio beyond which the resulting system of equations will be numerically unsolvable unless some form of “relaxation” parameters are used. This numerical ill-conditioning is in no way related to the shear locking phenomenon attributed to field- and edge-inconsistencies but is strictly related to the machine finite word length and it is termed “machine locking” by Briassoulis (1993). The study also showed that the deflection results are relatively insensitive to the choice of value of penalty parameters used to enforce the constraints of Eq. (41). A low value, however, will yield additional zero eigenvalues during the eigenvalue test. The through-the-thickness variations of displacements and stresses for the two laminated panels as predicted by the two LZZ3 models as well as the nine-node Lagrangian FSDT model are shown in Figures 60 to 88. Plots show comparison among the three finite element models and the three-dimensional elasticity solution for through- thickness variations of inplane displacement, u x , transverse deflection, uz , inplane stress 6 inplane shear stress I transverse shear stress In and transverse normal strain, xx’ xy’ 124 cu. For both the symmetric (sandwich panel) and unsymmetrical (CAV) layups with span-to-thickness ratio of a thousand, the FSDT as well as the LZZ3 models yield accurate results for the inplane displacement, u x (Figures 75 and 82) but not for the transverse shear stress I xz (Figures 79 and 86) with the exception of the three elements through-the- thickness LZZ3 model. In general, accurate inplane displacements will lead to accurate inplane stresses (see Figures 77-78, 80, 84-85, 87). However, for moderately thick laminates with span-to-thickness ratio of ten, the FSDT model is inadequate for both the sandwich and CAV panels. In the case of IJrz , predictions using the FSDT model with Whitney’s shear correction factors are very poor, as expected. As for the LZZ3 models, more than one element through-the-thickness is needed to capture the variation of IJrz in each of these rather complicated examples. The transverse shear stress shown in Figures 64, 72, 79 and 86, is computed using the three-dimensional constitutive relations with no shear correction factors. Predictions of the normal transverse strain, cu improve as more elements are used in the through-the-thickness direction. Nevertheless, as mentioned earlier, it is more common to recover the normal stress, 0'ZZ using the three-dimensional equations of equilibrium. 5.9 Summary The present LZZ3 model using the assumed strain field approach has been shown to be robust and accurate for application to thick and thin panels. It passes the membrane and Kirchhoff patch tests and has no spurious energy modes. This new finite element model is therefore a viable alternative to the elements derived using discrete layerwise or zig-zag theories as well as the conventional continuum elements. 125 1.035 E . —e— FSDT i A LZZ3:Isublamiuate 1-03 r — E1- - LZZ3: 3 sublaminates 1.025 u 1.02 Z ”z’exact 1.015 :- 1.01: I lows: ——————— g_ —r H h “"5" ~ ~ _ _ a: b 4 1) A A1 1 - - -..---s - A 2A..-.\/ . - -----y - - .....y 1 10. 100. 1000. 10000. a/h Figure 57. Normalized center deflection versus span-to-thickness ratio of a simply-supported square plate (isotropic) subjected to sinusoidal load 1.] if T fivfi T + FSDT : WA LZZ3:Isublaminate 1.08 . - El - LZZ3:3sublamiuates 1.06 L ”Z uz’exact l 04 L 1.02 ‘ fl- - q‘ n " nfig,“ vn w-s f" "fii 1 L . . - .4-..1 . g2 AAAHY - . -...-.<.{\ r . A_....y 1 10. 100. 1000. 10000. a/h Figure 58. Normalized center deflection versus span-to-tbickness ratio of a simply-supported square plate (sandwich) subjected to sinusoidal load 1.15:r “H": ° . A LZZ3:Isublatninate 1 C --El~- LZZ3:3sublaminates I 1.05 r . u » 1 z 1 _ r .............. _ n “”46" u El’ 1 z’exact 1 0.95 :1 g 0.9 E r 0.85 : . . 1 ’ l A /J\A AL A A A AAJ A A AAA! A A A A_L ‘ l 10. 100. 1000. 10000. a/h 126 Figure 59. Normalized center deflection versus span-to-thickness ratio of a simply-supported square plate (CAV) subjected to sinusoidal load v v v v 1 Y v v 1 1 v v v *v 1 v Y r 0.4 0.3 ~ 0.2 - 0.1 - 5 o h L -o.1 E, ,. i0 -0 2 _ """""" FSDT ’ ————— LZZ3: I sublaminate _03 1 — - — — LZZ3: 3 sublaminates 41.4 41.5 - - ----------------- -6 .4 Figure 60. Axial displacement versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 127 1L5 r ' ' " ' ' ' ' ' ' ' ' 1' si\" ' : 1 _— E’aSfiCi‘y \' ' / L 0'4 ---------- FSDT 1 , I 0.3 . — ' — ' - LZZ3: I sublaminate 1 /,/’ ' - — - - LZZ3: 3 sublaminates ' //’ A 214 A A 216 I 212 ‘ uzx10_6 A 210 ‘ Figure 61. 'D'ansverse deflection versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 115) ..- **I* * « IIIII f """ .Lj;;;.;' 4 0.4 0.3 r 1 0.2L 2 = h 10 0.1 r 5 0 -0.1 * Elasticity 1 ........... FSDT 412 ----- .Lzzssrsabuuuauuc ‘ .0.3L —’—— LZZ3:3sublaminates .0'4 A5 415- i‘ff'TILL- - .- ;- . -- --» ------------------ - -75 -50 -25 0 25 50 75 100 0x}? Figure 62. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 128 0.5 vvvvvvvvv 4 - - - ~ 4 - - 44 . ................... 0.49 . 7. - 10 0.465 / i h """""" FSDT ''''' LZZ3:1sublaminate 0.43 .4 — - — - LZZ3: 3 sublaminates 0.4 AAAAAAAAAAAAAAAAAAA 4- -75 -50 -25 0 25 50 75 100 0’” Figure 63. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load - top face sheet 0.3 . 0.2 - ‘ g = 10 1 1 E y 0.1 ’ 1 .1. o ' -0.1 T Elasticity """""" FSDT '0-2’ '. -'-‘- LZZ3:Isublaminate ; _03 _ 1‘ — - - - LZZ3: 3 sublaminates j -o-4 ..................... % ........... -05 .;..:..:.;.;I.:-.:'-:'. ...... :‘.-..-..: AAAAAAAAAAAAAAAAAA -1 5 -1 25 -1 -0.75 -0.5 -0 25 0 TIZ Figure 64. Transverse shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 129 0.5-{4' * 0.4 \ 0.3» 0.2- 2:10 0.1? E. 0 '0'1‘ Elasticity . ---------- FSDT .0'2 — ' - ‘ — LZZ3: I sublaminate .03. -——— LZZ3:3sublaminates 4 -0.4 \ 057‘1 -40 -20 0 20 40 Ix] Figure 65. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5 - - . fl 444444 - : - - - - 0.49 ’ a 7! — 10 0.465 E, 0.45 h : Elasticity ___q 0.43 ........... FSDT : ———— LZZ3: 1 sublaminate — — - - LZZ3: 3 sublaminates 0.4 A n . A_ A A n - A A A : A A A A a A A A A n A .40 -20 0 20 40 In Figure 66. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load - top face sheet 130 05’ Eta 0.4 '1, “‘\ i \. l "'3' ‘. —- Elasticity E ‘ 0.2 ~ ‘. —‘—‘— LZZ3:Isublattttnatc 1 ‘ — - - - LZZ3: 3 sublaminates 1 0.1 - 1). : z ‘ ' 1'. 0 1. : -0.1 r ‘- g 1 \ 1 -0.2 ~ 1. i 7... = 10 _03 . ‘- i ‘ 1 -0.4 '1. 4 4,5 31 --0----5----m---A15-- .5 Esz10 Figure 67. Transverse normal strain versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5—TV """ ;,'\'-_. """" r f f T f - r Y T 0.4 ’ 0.3 . Elasticity ....... FSDT 0.2 . ''''' LZZ3: 1 sublaminate ‘\ —— —- LZZ3:3sublamiuates 0e] 7 '. \ 5 0 - h -0.1 > Figure 68. Axial displacement versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 0 .5 _.‘ - . - . . - . - . - ......... ‘ . . e 0 4 ‘1‘ i '1 o ' " l . I . l 0'3 ‘-\ Elasticity / 0.2 , 1 ---- LZZ3: I sublaminate ,’ . - — - - LZZ3: 3 sublaminates // 0.1 ' '. / z ‘ ’l — 0 " 1 / / h t \ -0.1 - .' 1‘ ' 1 -0.2 L 1’ -6 g a 1 uz = 1523x10 (FSDT) , _ = 10 .03 . , ’1 h I -0.4 i .1. [I] -0.5--’---. ------------ I / ----- -.A 930 940 950 960 970 111x10-6 Figure 69. lransverse deflection versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 005 I. ''''''''''''''''' 0.4 r 03 t Elasticity ........... FSDT 0.2 ' """" LZZ3: I sublaminate — — - — L223.- 3 sublaminates / / l. I ‘75 Figure 70. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 132 o ............................... v [I ..... l : , / -0.1 -0.2 Z a _ r. 7. - 1° -0.3 / 7. .. I . ' --" _‘ Elasticity _ -0.4 1 . j _ _ _: ........... FSDT j, - ’ ''''' LZZ3: I sublaminate ' 1:11;} —— -— LZZ3: 3sublamiuates ‘0 5 AAAAAAA 1/ A A A 1 A" AAAAAAAAAAAAAAAAA L -75 -50 -25 0 25 50 75 Oxx Figure 71. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load - outer shell 0.4 r 0.31 Elasficily FSDT 0.2 r —" ' - LZZ3: I sublaminate 0 1 - - — - LZZ3: 3 sublaminates R'IN -0.l - -0.2 r -0.3 ~ -0.4 . -0.5 . - Figure 72. 'h'ansverse shear stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 34" 0.5 0.4 0.3 0.2 0.1 ~ -0.1 - -0.2 -0.3- -0.4 -0.5 133 V 1 V V le fi‘lfi Elasticity FSDT T " ' — ‘ - LZZ3: I sublaminate - — - — LZZ3: 3 sublaminates ‘ .10. Figure 73. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 0.5 v v 77 vvvvvvvvvvvvvvvvvvvvvvvvvvv . \\ 0.4 - ‘. ‘ \ \ . \ 0.3 ’ \. \ \ fl = 0.2 r -\ \ h 10 w 0.1 ~ \ \ . \ \ 5 o» ‘r / h -o.1 r -0.2 - Elasticity 4 _03_ ““ LZZ3:1sublaminate - — — - LZZ3: 3 sublaminates -0.4 - -0.5 * .............................. ‘ 20 30 40 50 60 —6 euxlo Figure 74. Transverse normal strain versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load I34 0.5 Kr'fi'flf'hn "f'r'fi"""""- 0.4 0.3’ 1 0.2' g: 3 h 10 0.1r . 5 0 h '0-1F Elasticity ........... FSDT '0'2' ----- LZZJ:Isublaminate _03 _ —- —- LZZ3:3sublaminatcs . '0e4 \ -0.5*.........L-.Ln.s ...A.ssan.s.....: -0.075 -0.05 -0.025 0 0.025 0.05 0.075 a I Figure 75. Axial displacement versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 05*? 1 f f f f T f : 'u‘ I | I 0.4 : i 1: 0.3»: g i I I 0.2“ 2:103 i l- h - l 0.1» : E 2 ° T: : .01. 9 :- ' Elasticity . : -: ........... FSDT . h '0'2 g - - - - - LZZ3: 1 sublaminate : { -0.3-E —--- LZZ3:3sublaminatcs i :. ’ . l -0.4 : i 4: -o.s-5-.---A.E--A.L-A..-!.-'- 51.9 51.95 52 52.05 u Z Figure 76. 'lransverse deflection versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 135 0.5V ' , ' ' ' r ' V V L V f - r f ' ' f '7" 0.4 l 0.3? ‘ .- J 0.2 7', = 103 0.1b * 5 0 '0-1' Elasticity ........ .. FSDT ‘ '02 —"‘- LZZJ:IsubIamiuate .03. -—-— LZZ3:3subIaminates . -0.4 / I -o.sh"_.:.rr.nn....-s...4 -l -0.5 0 0.5 l 0' Figure 77. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5..-7..--c--,v-,vaa 0.49 0.465 3:103 [ 5045 h h . Elasticity 0.43— ___________ FSDT "—‘— lZZ3:13ublalmuate' -~ -- L223: 3sublamiuates 0'4-1“‘k-o.s%flioku‘o.sé 1 xx Figure 78. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load - top face sheet 05- """ 1:; 1.? 1:.;.;;;;;_;:_<::_v . ~ . v . - M“ 0.4 .‘ l 0.3 . ,' 0.2 - (_l = 103 . h E 001 h J z o ‘ h .001 p Emmi” ........... FSDT '0‘2 ' u ----- LZZ3:1sublaminate 3 -0.3 - l r — -- LZZ3: 3 sublaminates * -l75 -150 425 -100 -75 -50 -25 0 Figure 79. 'D'ansverse shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5 * k i ' ' - ' , - . 9 0.4 1 0.3 ’ 0.2L 2 = 3 h 10 0.1 ~ . 5 0 h -0.1 ’ Elasticity """""" FSDT '0-2' “"'— LZZ3:Isublammat' e 4 .03_ ——-- LZZ3:3sublaminates -0.4 \ -o.s-.-Ahr....r an-..c.zai- -0.4 -0.2 0 0.2 0.4 I” Figure 80. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load 0.5_-\-I.--.V.V-H.flfl H--_-s-fvws. 0.4 '33 0.3L 'x.\ . 0.2r ‘-\_ 2: 3 J \\ h 10 091b ‘\' E 0 '\ \ h \ -0.1 \.\ Elasticity x. -0.2- —'--- LZZ3:1sublantinatc ‘-\. .03. --—- LZZ3:3sublamtnates \,\.\ ~0.4 “\.\3 -0.5-.-. his-.-u.c-c- U--.Jns+_u..-‘: -0.015 -0.01 -0.005 0 0.005 0.01 0.015 8 22 Figure 81. 'Il‘ansverse normal strain versus normalized thickness coordinate of a simply-supported square plate (sandwich) subjected to sinusoidal load sssssss v v v .,....,v v . oesr\ 1 d 0.4 r _ Elasticity 03 ........... FSDT 0.2 “‘—‘— LZZ3:Isublaminatc — - - - LZZ3: 3 sublaminates 0.1 Y FIN é / r .0.1 -0.2 - -0.3 - -0.4 -0.5~L.-.--.A-.s--s-.+---.-.n-.* -0.2 -0.1 0.1 0.2 0.3 10 #lfi T E: Figure 82. Axial displacement versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 05 -' ' * E ' a- fl . ~l. 1 l 0.4. , l‘ 0.3* 3 l : a _ l I 0.2. ,..—10 I E. 0.1- : i Z a . T. 0’ l i‘ . ' u. -0.1 I : .o.2- —— Elasticity . :. ----- FSDT I I .0.3- ----- LZZ3:Isublaminate :. — - — — LZZ3: 3 sublaminates ' I -0.4- I : 051.1.1-.a-11u.11...119..-:4 170.1 170.2 170.3 170.4 170.5 uz Figure 83. Transverse deflection versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 0.5_ ............ .,ss../.fi,.d 0.4’ - . Elasticity 0'3 """""" FSDT 0.2. —“"‘ LZZ3:1sublaminatc 0 1 — — -- LZZ3: 3sublaminates 5 o. / . h w [E 1 .02L 1 7 e f I j7 g = 103 '0o3' / 1 h -0.4- / ,’ / 0511-111..-//.-L-LT.-. ...... -300000 -200000 -100000 0 100000 200000 0’ xx Figure 84. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 139 o .................. - ........ c -0.1 I; a 3 / 1 ii = 1“ / 1 -0.2 / E. I -0.3 [I l / l Elasticity _0.4 f ........... FSDT _ 7 ""‘ LZZ3:1sublaminat¢ [ - - - - LZZ3: 3 sublaminates I _05 ..... WJ/l ....... -300000 -200000 -l00000 0 100000 200000 0 Figure 85. Inplane normal stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load - outer shell P u- 1 0.4 ' 0.3 . 0.2 ' 0.1 - 3'IN Elasticity ........... FSDT *‘—“ LZZ3:1sublaminatc 4 - - — — LZZ3: 3 sublaminates -0.l* /' -0.2r ‘x. -0.3+ °"\.\, -0.4~ 0‘. ‘. ‘. ‘- ‘. ‘. ~ ' - . I l A A A A I A A A L A A L A a j A A A A l A %L 4 -250 -200 -150 -100 -50 0 Figure 86. Transverse shear stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load 0.5; ........ C ............ 0.4- 0.3- g _ 103 0,2. h 0.l~ .2- 0b \ 4 h -0.1 . -0.2* . Elasticity '04“ ------ FSDT .04_ —-—-- LZZ3:1sublamiaate ' —--— LZZ3:3sublantiuatcs -0.5.1\ 450000 400000 -50000 0 50000 100000 1: xy Figure 87. Inplane shear stress versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load vvvvvvvvvvvvvvvvvvvvvv 0.5 0.4 0.3 0.2 ~ 0.1 - Elasticity """" LZZ3: 1 sublaminate - - - - LZZ3: 3 sublaminates a-IN -0.1 . -0.2 - '004 -0 5 . . AAAAAAAAAAAA Figure 88. Transverse normal strain versus normalized thickness coordinate of a simply-supported square plate (CAV) subjected to sinusoidal load Chapter 6 NON -LINEAR MODEL In this chapter, some aspects of the formulation and solution of geometrically non- linear problems will be treated. The study will concentrate on small strain/moderate rotation problems, typically encountered in slender, flexible structures. The solution will be based on a Total Lagrangian formulation, and the procedure adopted to trace out the complete load-deflection path is the Newton-Raphson method. 6.1 Introduction Most of the early work on geometric non-linearity related primarily to the linear buckling problem (see Crisfield, 1994 for his review on geometric non-linear problems). For genuine geometric non-linearity, ‘incremental’ procedures (or forward Euler) were originally adopted by Turner et al. (1960) using the ‘geometric stiffness matrix’ in conjunction with an updating of coordinates. For non-linear analysis, a number of ‘alternative’ strain measures are used. The common strains are the rotated engineering, Green’s, Almansi’s and rotated log strains. These strain measures are then used in the virtual work equation to develop the finite element model resulting in either the ‘Total Lagrangian’, ‘Updated Lagrangian’ or ‘Co- rotational’ formulations. These formulations are identical from a continuum mechanics point of view, the only difference being the choice of reference configuration. However, in 141 142 discretized form, applied to a structural element based on some approximating theory, the performances of the various formulations are usually quite different (see Mattiasson et al., 1986). Composite structures are designed such that their in-plane stiffnesses are almost always much greater than the transverse stiffnesses. It is therefore pertinent to consider only the von Karman non-linear strains (a special case of the Green’s strains) which will account for moderately large deflections and small strains. The inclusion of the von Karman non-linear strains will allow geometric non-linear analysis of laminated composite plates using the Total Lagrangian approach. The non-linear algebraic equations of the non-linear model are solved by iterative methods, which will be discussed in Section 7.3. In iterative methods, the expressions involving the dependent variables and their derivatives are computed using their values from the previous iteration, so that the integrals can be evaluated (by numerical integration methods). As in linear analysis, problems occur with shear locking and, again as with linear analysis, these can be improved by various methods such as the use of (selective) reduced integration or the assumed strain field approach as discussed in Chapter 3. We will once again adopt the approach of the assumed strain field approach for the layerwise plate element to prevent shear locking. 6.2 Theory Formulation The displacement field given in Eq. (62) with respect to the initial configuration is reproduced here for convenience: 143 a“) = 17(k)(X, Y)ff’(2) x a (n) __ ..(k) (k) k = 1, ...,7 uy _ ua (X, Y)‘Pa (Z) a = 1,2 (126) 14:") = wa(X, Y)Ma(Z)+w3M3(Z) with the shape functions (by), ‘l’g‘l, M a and M 3 defined in Appendix C. In general, the upper-case X, Y, Z and the lower-case x, y, 2 will relate to the initial and current coordinates respectively. In the present work, this distinction is not really necessary since the Total Lagrangian procedure always refers measures to the initial undeformed configuration. The stress measure that is conjugate to the Green’s strain tensor, E is the second Piola-Kirchhoff stress tensor, .5 and they are given as: 1.5 = %(F‘T-E—l); 12,, = %[%§%_8"] (127) and _ poaXi 8X1. 128 '7 — pdxksk’dxl ( ) p T S = B905") 4.5-(5"); S Using the second Piola-Kirchhoff stress tensor, .5 and the Green’s strain tensor, E as a suitable pair of energy conjugate stress and strain measures, the virtual work equation is given by: j SUM,j dVo = j tifiui dFo+ j pbifiui dVo (129) V F V 0 0 0 V 0 and F0 are the volume and surface area, respectively, in an arbitrary known reference configuration. The vectors, Q and g are the body forces per unit volume and the surface tractions, respectively in the initial undeformed state. If the in-plane strain terms, (e.g., (aux/8.302 and (Buy/8X)2 ), are small compared 144 to one, then they can be considered negligible, and the Green’s strain tensor, E of Eq. (127) degenerates to the von Karman strains. They can be expressed as: E Bux 1 BuZ 2 xx - 327 +502?) w au u - y z 2Eyz 5—2- +37 2E),z Variation of Eq. (130) leads to: dfiux 8Eix:=:;§ dfiu SEW = 37 y afiuz SEZZ = -—Z' dSu ZSEYZ = 7 y afiux 25Exz = a—Z dfiux 26Exy = 5'? 2 duz ) Ezz=a—Z Bu _ x icy-37+ au ___y 1(a_“z 'ay+2ay Bu auz I first (130) Bu __y+ 8X 2E (auz auz $53? ) duz dfiuz + (a7 a7 (BuZ dfiuz ) ) +5??? (131) Bfiuy duz dfiuz a")? ( asaz duz +_.__ ax aY ( +3737 ) l The constitutive relations relating the Lagrangian-Green’s strain and second Piola- Kirchhoff stress for the case of small strains and large rotations are of a similar form to Eq. (51). The proof is outlined in Appendix E. 145 F - PS w.) Cfi’j’c‘é’c‘g’ 0 0 Cfi’g’ -E '(k) x" C“) C“) 0 0 C“) "’ Syy 22 23 26 E” (k) (k) 53, = C33 0 0 C36 Ezz (132) SYZ C2? C3? 0 2Eyz S 2E XZ Sym C(Slg) O 12 SM, 2E , _ ._ (k) _ n. where the index k is to denote the layer number. If the non-linear (von Karman) strains, Eq. (130), with the constitutive relations, Eq. (132), are to be included in the virtual work equation, Eq. (129), the procedure discussed in Section 4.3 does not change, except that the resulting algebraic equations will be non- linear. The Principle of Virtual Work yields: 0- Mail" My)? ems)? await" ’I" 3X+8Y+BY+6X - (Raff) - (R,)ff’16a;"’ (9,), ale), ' ax + 31’ - (Qz)a]5wa + (flaws (133) a BWB a BWB a aw, ‘[a—x[(”x)asa—x ] * 517W mm; ] +a—x[(Nyx>as'a—y ] 3w + a—aY-[(ny)aB-a7p]]wa]d.4 + may) (1,3 = 1,2;k =1,2,4,...,7;l =1,...,7 where Lag”) are the boundary terms. The stress resultants are as defined in Eqs. (66, 67) except that the stresses are replaced by the second Piola-Kirchhoff stresses and the reference configuration is the initial undeformed configuration. The additional stress resultants that are not defined by Eqs. (66, 67) are: 146 N i (Ix/x)”, = 2 053m, Y, Z)Ma(Z)MB(Z)dZ] i=1 0 N i (Ny)aB = 2 OS”(X, Y, Z)Ma(Z)MB(Z)dZ] (134) i=10 N i (”Mala = 2 styor, Y, Z)Ma(Z)MB(Z)dZ] = (Nyx)aB i=1 0 The boundary and loading conditions for the plate are: mag“) = — jun, + m2);"’5ag" dA — j (Ml B + M23);k)n38ilflk) as A 32 L(wa) = — I[qa5wa + 315%] dA — j [(Qn)a8wa + Qnaw3] as A 52 (135) —j [(Nn)anB,n + (an)anB,s]5wa dS s 0t,|32=1,2;k = l,2,4,...,7 where the definitions of Eqs. (70, 71) were used after again replacing the stresses with the second Piola-Kirchhoff stresses. The subscripts n and 5 refer to the outward normal and tangential directions respectively and they are with respect to the initial undeformed configuration. Note that the following expressions were used to obtain the final form of Eq. (135): N" = 11513, + 2nylnxlny + Nylfiy NS = (Ny' Nx)lnxlny + ny(13X +13y) .3. = l 5.3—-1 9— (136) ax nx an "Y 83 a _ a a 3): — lny 524-1,” .a—S where [M and I"), are the direction cosines of the outward normal of the boundary in the initial configuration. 147 The finite element model for the layerwise plate theory is obtained by substituting the element approximations, Eq. (137), into the virtual work equation, Eq. (133). u = Piui v = Pivi w = PlWi+in9xi+Nyi9yi 6x = Pie”. 0y = Piflyi (137) 1,, = Pit“. ry = Pity. l = 1, .. ,4 Pi are the bi-linear Lagrangian functions and N N yi are defined in Appendix B. xi’ The finite element model can then be expressed in a matrix form, [Ke]{de} - {Fe} = 0 (138) where [K e] , {d‘} and {F e } are the elemental stiffness matrix, nodal degrees-of- freedom vector and generalized force vector, respectively (see Appendix F for details). Similar to the linear layerwise finite element model, the inconsistent transverse shear strain, normal strain and penalty constraint fields must be replaced by their consistent counterparts, Eqs. (118—119, 123-125) in the non-linear model to eliminate the shear locking and Poisson’s ratio stiffening effects. 6.3 Newton-Raphson Iteration Method The non-linear algebraic equations of the non-linear model are solved by iterative methods. The element stiffness matrix in Eq. (138) is non-linear and unsymmetric when the non-linear von Karman strains are included. Therefore, the assembled non-linear equations will be non-linear and unsymmetric. The assembled non-linear equations must be solved, after imposing the boundary conditions, by an iterative method, which seeks an approximate solution to the algebraic equations by linearization. The two most common 148 methods are the direct iteration and the Newton-Raphson iteration methods. The latter method is probably the most widely used because of the symmetry of the tangent stiffness matrix obtained. The Newton-Raphson iterative scheme is described in many books (see for example, Ochoa and Reddy (1992); Crisfield (1994)). The method is based on a Taylor series expansion of the algebraic equations, Eq. (139), about the known solution, {d}i at the ith iteration. The global finite element equations can be written as: [K({d})]{d}-{F} = 0 (139) For an approximate solution {d}i : {R} = [K]{d}-{F}¢0 (140) where {R} is called the residual and is generally not equal to zero unless convergence has been attained. Expanding {R} about {d}i , we obtain {R(d”‘)} = {R(di)}+[3§§i]'<{d}i+l—{d}") (141) 32{R} l'+1 l' 2 +2’—,[-a——{d},] ({d} —{d}) + Since our goal is to minimize the residual {R(d‘)} , we can write Eq. (141) as: 0~ {R(d") )i}+[KT ]{Ad"}+0({Ad"})2 (142) where [K T] is the tangent stiffness matrix (or geometric stiffness matrix), .- _ .- _ 3{R} ‘ [KT] - [KTM )1 - [am] (143) (the stiffness coefficients of [K T] are listed in Appendix F). If the finite element model is derived using a variational principle, it can be shown that the tangent stiffness matrix, [K T] , will always be symmetric even if the direct stiffness matrix, [K] , is symmetric for 149 non-linear problems. Neglecting quadratic terms and higher, the truncated Taylor series can be rewritten as: {Adi} = -[K"T]"{R(di)} . _1 . . (144) = -[K1(d')] (([K(d‘)]{d'} - {F})) and the total solution at the (i + 1 )th iteration is given by: {di+ 1} = {d‘} +{Ad’} (145) The iterative solution of Eqs. (144, 145) is continued until a predetermined criterion and tolerance, 8 (normally chosen in the range, 10.5 S 8 S 10—1 ) is satisfied. Two typical error critieria are: . , , . _ ({Adi}T{Adi})l/2 Displacement cntenon. IIdII2 _ ({di}7{di})1/2 < (146) - - - . _ ({R’}T{R’})“2 Force (orResrdual)cr1ter10n. IIRII2 — ({F}T{F})1 /2 S There are many other convergence criteria such as the energy norms, but the above two criteria often are the most useful. The method described above yields the solution for a given load, but it is often desirable to know the load-deflection path. Also, these methods may be inefficient, even non-convergent, for highly non-linear problems. It is therefore desirable to use incremental strategies to solve most non-linear structural problems. Eq. (139) can be rewritten as: [K({d})] = MFO} (147) where {F 0} is a reference load and A. is a scalar parameter. The methodology does not change except that the load is applied in increments. For the kth increment, we can rewrite Eq. (144) as: 150 {Adi} = -[KT(dl)1;‘<<[K(dl)1{dl}-x,.{F,})) (148) where {d)‘;+'} = {61%;} +{Adi} (149) A, = 1,,_,+A)., and {d3} is the solution at the end of the last increment. In general, it is good to have three to six iterations per increment. This is controlled by the load increment AM. Note that at the beginning of each iteration the tangent stiffness matrix and residual vector must be updated using the latest available solution, {di} or {dfi} . Rather than updating the tangent stiffness during every iteration, it is common to make the approximation: [Kidk = [thc (150) and as such update [K T] only at the beginning of each increment. Since the tangent stiffness matrix is kept constant and only the residual vector is updated during the iterations within each increment, there will be computational savings. Such an approach is called the (incremental) modified Newton-Raphson method. 6.4 Numerical Results Geometric non-linear responses of composite plates can be very different from isotropic, metallic plates. The responses of composite plates are highly dependent on the lamination scheme and boundary conditions and can be very significant even at small loads and deflections. Two laminates were selected to assess the accuracy of our model. The first is a single-layer orthotropic plate and the second is a cross-ply (0/90/90/0) _LI-rlm 151 laminate. The cases selected are attractive because test results exist for the non-linear behavior of square plates with all four edges either simply supported or clamped under uniform transverse loading (see Zaghloul and Kennedy, 1975; cited in Ochoa and Reddy (1992)). The geometries and material properties for the 0° ply that were used in the experiments are as follows: (i) Orthotropic plate: a=b=12,h=0.138 E, = 3x106, E2 = E3 = 1.28x106 6 (151) (ii) Cross-ply laminate (0/90/90/0): a=b=12,h=0.096 E, = E3 = 1.8282x106, E2 = 1.8315x106 6 (152) v,2 = v23 = v,3 = 0.2395 A quarter plate model will be used to analyze the two cases. Different meshes were used to determine the convergence characteristics of the new layerwise zig-zag element, LZ23 for the two cases. In Figures 89 and 92, the load-deflection curves of different inplane discretizations are shown. A mesh of 4x4x1 for both cases yielded converged results in both cases. Through-thickness refinement using four sublaminates (each representing a physical layer) reinforced the fact that a mesh of 4x4x1 is adequate to model both cases. The load—deflection curves for the through-thickness refinement for both cases are shown in Figures 90 and 93. In Figures 91 and 94, the present theory was compared with solutions obtained by the 152 classical laminated plate (CLPT) and the first-order shear deformation (FSDT) theories in addition to the test results obtained by Zaghloul and Kennedy (1975). A mesh of 8x8xl was used for the comparison. For the simply-supported orthotropic plate, LZZT as well as the shear deformation element yield results that are in excellent agreement with the experimental data (see Figures 89 to 91). The classical laminate bending element under- predicts deflection. The load-deflection curves obtained from the quarter-plate model for the clamped, cross-ply laminate are shown in Figures 92 to 94. The finite element solutions of the LZZ3 model under-predict the deflection (almost 15% at the maximum load of 2) when compared with the experimental data. This trend is also seen in the shear deformation element (see Figure 94). In fact, the load-deflection curves by the two different theories are almost identical. Since mesh refinements of the LZZ3 model showed that the results are converged. This therefore lead us to the same conclusion as Ochoa and Reddy (1992). They attributed the poor prediction of their shear deformation element to be due to inaccuracies in the representation of the material stiffnesses as well as the actual boundary conditions of the experiment. 6.5 Summary The prediction of the load-deflection path by the non-linear LZZ3 model does not differ much from the shear deformation element of Ochoa and Reddy (1992). This is to be expected since the span-to-thickness ratios of the two experiments are large. In Chapter 5, we have demonstrated that the first-order shear deformation element perform adequately for large span-to-thickness ratios except in cases of laminates with drastically different materials in adjacent layers or large number of plies. Therefore, in order to evaluate the 153 true performance of the non-linear LZZ3 model, one must also compare the in-plane displacements as well as the (Cauchy) stresses. However, because of lack of experimental data we are unable to do so. 154 05,-.--fs---,----,-ss D — Experimental f‘: I < . —e~— 4x4x1 A .-/' 0.4 ' A 8x8x1 ,..! . — a— I6x16x1 , linear ,3, , 0.3 r A "/ 1 ”max ,. non-linear 1 1'5; ,... v 0'2 ' ,;. ‘ /€‘/’: + 0.1 ’ 4., I ; / 1 ' I 0 - T 0 0.5 1 1.5 2 load, 2» Figure 89. Load-deflection curves for the bending of a simply-supported, square, orthotropic plate under uniform transverse loading - effect of inplane mesh refinement 0.5 ' ' i ' ' v v v v v 1 f v v V r ‘ J — .A. Expenme' ntal 0 4 ’ —9— 4x4x1 ' 4:: / 4. O b ' . . '/ . "A " 4x4x4 _-/ , 1 —EJ- — 8x8x1 . ._./.-,‘ J C — «a— 8x8x4 1W” "/ i 0.3 L -/ “ b ' é ' w 4 max .r‘" non-linear _ .,— . h »/ t‘ ”’5 i v- ‘ 0.2 _ / 1/ L a, .A i 1 a 0.1 - .,, . l '4 4 / J i ‘ J o A A g A A A A A A a A A A A r A A A A 1 0 0.5 1 1.5 2 load, 1. Figure 90. Load-deflection curves for the bending of a simply-supported, square, orthotropic plate under uniform transverse loading - effect of through-thickness mesh refinement 155 0.5 F 9‘ b ‘ . : —- Expenmental f. (- " . b 0.4 - ‘9_ CL" I, _ [:1 —-El— LZZ3 , ff . linear 0.3 - a . ,...” w "" r" . . W“ e. a non-linear - A. 0 .4”.- ,.M"" a 4 0.2 - .,.- ,... , 1 .' .02 . I , m" ‘ 3 g Q . L ._ - e v , L '43; V 1 I u r ~/¢ < 001 . 1,;/% 1 0 . A A A A A A A A A A AA A A AA 1 A A A . s 0 0.5 l 1.5 2 load, 2» Figure 91. Load-deflection curves for the bending of a simply-supported, square, orthotropic plate under uniform transverse loading - comparison among different theories and experimental results 0.5 f a V V V Y ' ‘ f V T if ' ' V V ' " r T h . 1 . _— Experimental 0 4 ’ —9— 4x4x1 A, . - '71 . . --¢.-- ‘ 8x8x1 A .,v’ , ’ —- E1- — 16x16x1 / . . ._ I»; ’ . 0,3 - linear ,-” . . , . I w I max 1 ‘3? 1 L / 0-2 t 4" non-linear t A .3. ‘ i .4.» .— 5;- " .c. .— ) ’3 ,A i‘ i“ 001 P ’. 5 r ’3- " ) /’~’ 1 0 . A A A A J A A A A L A A A A L A A A A l 0 0.5 1 1.5 2 load, 1 Figure 92. Load-deflection curves for the bending of a clamped, square, symmetric cross-ply (0/ 90/ 90/ 0) plate under uniform transverse loading - effect of inplane mesh refinement 156 0.5 a f a e r a fi . I —— Experimental I . —e— 4x4x1 fl . 0'4 f A 4x4x4 _ 4’0 J ’ —EJ- — 8x8x1 .,Iwfi’ - I — «tr — 8x8x4 , f ”’ 0.3 ‘ linear ,2” ”max . 7; ' i I, 2 0.2 . non-linear J A. - -19.: - 2 .3. ‘ ’ J 0.1 - , , .. ~ J /: J 1 0 ’ A A g 1 1 A A A A n A A A A n J 0.5 1 1.5 2 load, 1 Figure 93. Load-deflection curves for the bending of a clamped, square, symmetric cross-ply (0/90/90/0) plate under uniform transverse loading - effect of through-thickness mesh refinement 0.5 fl . - . f , . , J —-—-— Expenme’ ntal " : 0.4L CL” ” .293 . A FSDT ééf -, A ’ _’B" LZZ3 .’_:.)‘-' “J J . ., . 0.3 ’ linear e __./Q J w : ’. 7'9 J max C ’. / :‘fir non-linear A. ,J 0_2 . [.91 ¢ 6 ' . r I e, . : :‘2' ’4 ...... ' 'A“ _xéqu—“A J /, flex-W'A’gé" } 0'1 T /fl, “Aer-J" " 1 ; i". 1 0 0.5 1 1.5 2 load, 1 Figure 94. Load-deflection curves for the bending of a clamped, square, symmetric cross-ply ( 0/ 90/ 90/ 0) plate under uniform transverse loading - comparison among different theories and experimental results Chapter 7 CONCLUSIONS Two main objectives were accomplished —- the development of a new technical theory and a new, robust C0 plate finite element for sandwich and composite panels. The new technical theory has some common features of other approaches of laminate theories, in particular the discrete layerwise and the zig-zag theories. However, in the formulation presented here, the use of surface quantities to represent the degrees-of-freedom as well as the satisfaction of the non-homogeneous transverse shear traction at the top and bottom of the sublaminate is a significant departure from the other two types of theories. A review of the most relevant aspects of those differences are summarized in the next section. 7.1 Review Our technical theory as well as the discrete layerwise and zig-zag theories are assumed displacement approaches. Discrete layerwise theories do not generally satisfy the interlaminar continuity of transverse shear stresses explicitly (see Li and Liu (1996) for an exception to this) while both the new theory and the zig-zag theories do satisfy these conditions exactly. However, the discrete layerwise theories will achieve this in the limit with increasing number of sub-divisions. The concept of stacking elements through the thickness is readily achieved by elements derived from discrete layerwise theories as the only condition that must be met is 157 158 the continuity of displacements. This is because there is no specific requirement for the continuity of interlaminar transverse shear stresses. The zig-zag theories which have their reference planes at the mid-planes cannot be stacked. The elements using zig-zag theories must meet both the displacement continuity as well as the continuity of transverse stresses. This cannot be done since the zig-zag theories satisfy only homogeneous shear traction (traction-free) conditions at the top and bottom surfaces. Any stacking will only result in ‘zero’ or discontinuous shear traction between two sublaminates. In the present theory, the non-homogeneous shear traction conditions are satisfied exactly. This is achieved through the shear traction degrees-of-freedom that can be left unknown at the interface where two sublaminates are stacked on top of one another. The new theory allows the development of an eight-node brick finite element model because it uses surface degrees-of-freedom. Essentially, the theory is still based on plate kinematics and uses an approach that is in many ways opposite to what is normally done, that is, to degenerate a three-dimensional body into a two-dimensional surface. Casting of the technical theory in terms of surface quantities, however, has the advantage of allowing the model to represent exactly how the surface tractions are applied. In the case of the zig- zag theories, surface tractions are applied at the reference surface which corresponds to the mid-surface. This is acceptable for thin plates but not for thick plates. As observed in our numerical experiments (see also Bogdanovich (1991)), in a thick plate with symmetric laminate scheme, the through-thickness distribution of in—plane displacements, stresses and transverse stresses are non-symmetric with respect to the mid-surface. This is due to the presence of transverse normal stress in the plate. The assumption that transverse displacement is independent of the through-the- 159 thickness coordinate must be used with caution. The maximum deviation of transverse displacement from the constant value is in the range of i12 — 18% for laminated as well as sandwich panels of up to 55 layers and having span-to-thickness ratio of 4 as shown by our numerical computations. This deviation narrowed to less than 2% for span-to- thickness ratio of 10. Nevertheless, this assumption generally is acceptable in many practical problems if the span-to-thickness ratio is large and the material properties are not drastically different as in sandwich panels. However, neglecting the transverse normal stress through the assumption of inextensible normal strain will result in the inability to predict failure modes initiated by this stress component. The choice of a quadratic displacement field in the thickness direction for our theory overcomes this problem and is adequate as revealed by comparisons between the exact solutions (Pagano, 1969; and Burton and Noor, 1994) and our finite element computations. The layerwise theories result in elements with a large number of degrees-of-freedom. The finite element model derived using the new theory only has a small and fixed number of degrees-of-freedom while retaining the desired accuracy just as in the zig-zag theories. This is very attractive and efficient for large computational models. Furthermore, our theory will yield traditional engineering (displacement and rotation) degrees-of-freedom. These engineering degrees-of-freedom will allow easier imposition of boundary and loading conditions. In summary, it has been demonstrated through numerous numerical tests that the finite element model derived using this new theory is a viable alternative to the elements derived using discrete layerwise or zig-zag theories as well as the conventional continuum elements. It passes the patch test and is accurate and robust. 160 7.2 Areas of Future Research The approach to using surface quantities to ‘regenerate’ plate finite elements into three-dimensional elements is new and attractive. Some areas of current research include extending the theory to the generalized shell theory as well as to heat transfer theories. There is therefore a lot of potential for future research to expand this concept. As for the derivation of finite element models, there is a need to consider things such as material non-linearity and thermal effects, plasticity (in metal forming applications) and damage (for delamination problems). The development of a six-node wedge element using this theory is also needed. The reason is that for some geometric boundary types, the wedge element types sometimes are preferred over the eight-node brick element type. An important area of future research is to eliminate the restriction of the existing finite element formulations to flat plates. This is because the existing finite element has the restriction that the thickness coordinate must be straight. The curved shell element under development will alleviate this problem but yet still does not allow modeling of tapered structures or plates with variable thickness. Therefore, the formulation of a more general topology for the plate element will create greater flexibility in the model. Accurate transverse normal stress is important for analyzing crack problems or problems that have failure modes associated with this stress component. The model accurately predicts displacements and stresses other than the transverse normal stress. Another area of possible research is therefore to improve the transverse normal stress prediction in the model. The new finite element model has potential for use in a lot of applications. The element can be used for explicit and implicit dynamic analysis, structural optimization as 161 well as metal forming processes. In addition, this element will be able to perform damage analysis with the availability of new adaptive techniques like the superposition method (Fish, 1992; Kim et al., 1991) to improve the quality of finite element calculations in regions of unacceptable errors. BIBLIOGRAPHY BIBLIOGRAPHY Actis R.L. and Szabo B.A., 1993, “Hierarchical Models for Bi-directional Composites”, Finite Elements in Analysis and Design, Vol. 13, No. 2/3, pp. 149-168. Aminpour M.A., 1990, “A Four-Node Assumed-Stress Hybrid Shell Element with Rotational Degrees of Freedom”, NASA Contractor Report 4279. Averill R.C. and Reddy J.N., 1992, “An Assessment of Four-Noded Plate Finite Elements Based on a Generalized Third-Order Theory”, International Journal for Numerical Methods in Engineering, Vol. 33, pp. 1553-1572. Averill R.C., 1994, “Static and Dynamic Response of Moderately Thick Laminated Beams with Damage”, Composites Engineering, Vol. 4, pp. 381-395. Averill R.C. and Yip Y.C., 1996a, “Development of Simple, Robust Finite Elements Based on Refined Theories for Thick Laminated Beams”, to appear in Computers and Structures. Averill R.C. and Yip Y.C., 1996b, “An Efficient Thick Beam Theory and Finite Element Model with Zig-Zag Sublaminate Approximations”, to appear in AIAA Journal. Babuska I., Szabo B.A. and Chang R.L., 1992, “Hierarchial Models for Laminated Composites”, International Journal for Numerical Methods in Engineering, Vol. 33, No. 3. PP. 503-535. Barlow J., 1976, “Optimal Stress Locations in Finite Element Models”, International Journal for Numerical Methods in Engineering, Vol. 10, pp. 243-251. Barbero EJ. and Reddy J.N., 1990, “An Accurate Determination of Stresses in Thick Laminates Using a Generalized Plate Theory”, International Journal for Numerical Methods in Engineering, Vol. 29, pp. 1-14. Barbero E.J., Reddy J .N. and Teply IL, 1990, “General TWO-Dimensional Theory of Laminated Cylindrical Shells”, AIAA Journal, Vol 28, No. 3, pp. 544-553. 162 163 Barbero E.J., 1991, “3-D Finite Element Modeling of Laminated Composites by Layerwise Constant Shear Elements”, in Enhancing Analysis Techniques For Composite Materials, NDE-Vol. 10, Winter ASME Conference 1991, Eds. Schwer L., Reddy IN. and Mal A., pp. 229-237. Bathe KI. and Dvorkin EN, 1986, “A Four-Node Plate Bending Element Based on Mindlin/Reissner Plate Theory and a Mixed Interpolation”, International Journal for Numerical Methods in Engineering, Vol. 21, pp. 367-383. Belytschiko T. and Hsieh E.J., 1973, “Non-Linear Transient Finite Element Analysis with Convected Co-ordinates”, International Journal for Numerical Methods in Engineering, Vol. 7, pp. 255-271. Belytschiko T. and Glaum L.W., 1979, “Applications of Higher Order Corotational Stretch Theories to Nonlinear Finite Element Analysis”, Computers and Structures, Vol. 10. pp. 175-182. Bogdanovich A.E., 1991, “Spline Function Aided Analysis of Inhomogeneous Materials and Structures”, in Local Mechanics Concepts for Composite Material Systems, IUTAM Symposium, Blacksburg, VA., Eds. Reddy J.N., Reifsnider K.L., Springer Verlag, pp. 355-382. Briassoulis D., 1993, “The Performance of a Reformulated Four-Node Plate Bending Element in Moderately Thick to Very Thin Plate Applications”, Computers and Structures, Vol. 47, No. 1, pp. 125-141. Burton W.S. and Noor A.K., 1994, “Three-Dimensional Solutions for Thermomechanical Stresses in Sandwich Panels and Shells”, Journal of Engineering Mechanics, ASCE Vol. 20, 10, pp. 2044-2071. Cho M. and Parmerter R.R., 1993, “Efficient Higher Order Composite Plate Theory for General Lamination Configurations”, AIAA Journal, Vol. 31, pp. 1299-1306. Cho Y.B. and Averill R.C., 1996, “An Improved Theory and Finite Element Model for Laminated Beams Using First Order Zig-Zag Sublaminate Approximations”, submitted for publication. Chattetjee SN. and Kulkarni S.V., 1979, “Shear Correction Factors for Laminated Plates”, AIAA Journal, Vol. 17, No. 5, pp. 498-499. Crisfield M.A., 1983, “A Four-Noded Thin-Plate Bending Element Using Shear Constraints - A Modified Version of Lyons’ Element”, Computer Methods in Applied Mechanics and Engineering, Vol. 38, pp. 93-120. 164 Crisfield M.A., 1990, “A Consistent Co-Rotational Formulations for Non-Linear, Three-Dimensional, Beam-Elements”, Computer Methods in Applied Mechanics and Engineering, Vol. 81, pp. 131-150. Crisfield M.A., 1994, Non-Linear Finite Element Analysis of Solids and Structures, Vol. 1: Essentials, John Wiley and Sons. DiSciuva M., 1985, “Development of an Anisotropic, Multilayered, Shear-Deformable Rectangular Plate Element”, Computers and Structures, Vol. 21, pp. 789—796. DiSciuva M., 1987, “An Improved Shear-Deformation Theory for Moderately Thick Multilayered Anisotropic Shells and Plates”, Journal of Applied Mechanics, Vol. 54, pp. 589-596. DiSciuva M., 1986, “Bending, Vibration and Buckling of Simply Supported Thick Multilayered Orthotropic Plates: An Evaluation of a New Displacement Model”, Journal of Sound and Vibration, Vol. 105, pp. 425-442. DiSciuva M., 1993, “A General Quadrilateral Multilayered Plate Element with Continuous Interlaminar Stresses”, Computers and Structures, Vol. 47, pp. 91-105. Dvorkin E.N. and Bathe K.J., 1984, “A Continuum Mechanics Based Four-Node Shell Element for General Non-Linear Analysis”, Engineering Computing, Vol. 1, pp. 77-88. Epstein M. and Glockner PG, 1977, “Nonlinear Analysis of Multilayered Shells”, International Journal for Numerical Methods in Engineering, Vol. 13, pp. 1081-1089. Epstein M. and Huttelmaier HR, 1983, “A Finite Element Formulation for Multilayered and Thick Plates”, Computers and Structures, Vol. 16, pp. 645-650. Fish J, 1992, “The s-Version of the Finite Element Method”, Computers and Structures, Vol. 43, NO. 3, pp. 539-547. Flanagan G., 1994, “A General Sublaminate Analysis Method for Determining Strain Energy Release Rates in Composites”, AIAA-94-1356-CP, presented at the AIAA 1994 Conference, pp. 381-389. Fried 1., 1974a, “Finite Element Analysis of Incompressible Materials by Residual Energy Balancing”, International Journal of Solids and Structures, Vol. 10, pp. 993-1002. Fried 1., 1974b, “Residual Energy Balancing Technique in the Generation of Plate Bending Finite Elements”, Computers and Structures, Vol. 23, pp. 409-431. Friedman Z. and Kosmatka J .B., 1993, “An Improved Two-Node Timoshenko Beam Finite Element”, Computers and Structures, Vol. 47, pp. 473-481. 165 Hinton E. and Huang H.C., 1986, “A Family of Quadrilateral Mindlin Plate Elements with Substitute Shear Strain Fields”, Computers and Structures, Vol. 23, pp. 409-431. Huang H.C. and Hinton E., 1984, “A Nine-Node Lagrangian Mindlin Plate Element”, Engineering Computing, Vol. 1, pp. 369-379. Hughes T.J.R., Taylor R.L. and Kanoknukulchal W., 1977, “A Simple and Efficient Finite Element for Plate Bending”, International Journal for Numerical Methods in Engineering, Vol. 11, pp. 1529-1543. Hughes T.J.R. and Tezduyar TE, 1981, “Finite Elements Based Upon Mindlin Plate Theory, with Particular Reference to the Four-Node Bi-linear Isoparametric Element”, Journal of Applied Mechanics, ASME, Vol. 48, pp. 587-596. Jang J. and Pinsky P.M., 1987, “An Assumed Co-variant Strain Based 9-Noded Shell Element”, International Journal for Numerical Methods in Engineering, Vol. 24, pp. 2389-241 1. Jiang J. and Olson M.D., 1994, “Large Elastic—Plastic Deformations of Slender Beams: Co-Rotational Theory vs. von-Kérman Theory”, Computational Mechanics, Vol. 159 pp- 117‘128. Kim Y.H. and Lee S.W., 1988, “A Solid Element Formulation for Large Deflection Analysis of Composite Shell Structures”, Computers and Structures, Vol. 30, No. 1/2, pp. 269-274. Kim Y.H., Levit I. and Stanley G., 1991, “A Finite Element Adaptive Mesh Refinement Technique that Avoids Multipoint Constraints and Transition Zones”, CED- Vol. 4, Iterative Equation Solvers for Structural Mechanics Problems, ASME, pp. 27-35 Kong J. and Cheung Y.K., 1995, “Three-Dimensional Finite Element Analysis of Thick Laminated Plates”, Computers and Structures, Vol. 57, No. 6, pp. 1051-1062. Ling-Hui H., 1994, “A Linear Theory of Laminated Shells Accounting for Continuity of Displacements and Transverse Shear Stresses at Layer Interfaces”, International Journal of Solids and Structures, Vol. 31, pp. 613-627. Lee CY. and Liu D., 1992, “Layer Reduction Technique for Composite Laminate Analysis”, Computers and Structures, Vol. 44, No. 6, pp. 1305-1315. Lee SW. and Pian T.H.H., 1978, “Improvement of Plate and Shell Finite Elements by Mixed Formulations”, AIAA Journal, Vol. 16, pp. 29-34. Li X. and Liu D., 1996, “A Laminate Theory Based on Global-Local Superposition”, to appear in Communication of Numerical Methods in Engineering. 166 Li X. and Liu D., 1996, “A Generalized ZigZag Theory for Laminate Plate Analysis”, to appear in AIAA Journal. Liu D., 1995, “2in Theory for Composite Laminates”, AIAA Journal, Vol. 33, No. 69 pp. 1163‘1165. Lo K.H., Christensen R.M. and Wu E.M., 1977, “A High-Order Theory of Plate Deformation - Part 2: Laminated Plates”, Journal of Applied Mechanics, Vol. 44, pp. 669- 676. Lu X. and Liu D., 1992, “An Interlaminar Shear Stress Continuity Theory for Both Thin and Thick Composite Laminates”, Journal of Applied Mechanics, Vol. 59, pp. 502- 509. Lynn RP. and Arya SK, 1974, “Finite Elements Formulated by the Weighted Discrete Least Squares Method”, International Journal for Numerical Methods in Engineering, Vol. 8. pp. 71-90. MacNeal R.H., 1982, “Derivation of Element Stiffness Matrices by Assumed Strain Distributions”, Nuclear Engineering and Design, Vol. 70, pp. 3-12. MacNeal RH. and Harder R.L., 1985, “A Proposed Standard Set of Problems to Test Finite Element Accuracy”, Finite Elements in Analysis and Design, Vol. 1, pp. 3-20. Mattiasson K., Bengtsson A. and Samuelsson A., 1986, “On the Accuracy and Efficiency of Numerical Algorithms for Geometrically Nonlinear Structural Analysis”, in Finite Element Methods for Nonlinear Problems, Eds. Bergan P.G., Bathe K]. and Wunderlich W., Springer-Verlag Berlin, Heidelberg Mohan P.R., Naganarayana B.P. and Prathap G., 1994,“Consistent and Variationally Correct Finite Elements for Higher-Order Laminated Plate Theory”, Composite Structures, Vol. 29, pp. 445-456. Morley L.S.D., 1963, Skew Plates and Structures, Pergamon, Oxford. Murakami, H., 1986, “Laminated Composite Plate Theory with Improved In-Plane Responses”, Journal of Applied Mechanics, Vol. 53, pp. 661-666. Naganarayana B.P., Prathap G., Dattaguru B. and Ramamurty TS, 1992, “A Field- Consistent and Variationally Correct Representation of Transverse Shear Strains in the Nine-Noded Plate Element”, Computer Methods in Applied Mechanics and Engineering, Vol. 97, pp. 355-374. Noor AK. and Andersen C.M., 1977, “Mixed Isoparametric Finite Element Models of Laminated Composite Shells”, Computer Methods in Applied Mechanics and Engineering, Vol. 11, pp. 255-280. 167 Ochoa 0.0. and Reddy J .N., 1992, Finite Element Analysis of Composite Laminates, Kluwer Academic Publishers, Netherlands. Pagano N.J., 1969, “Exact Solutions for Composite Laminates in Cylindrical Bending”, Journal of Composite Materials, Vol. 3, pp. 398-411. Pawsey SE. and Clough R.W., 1971,“Improved Numerical Integration of Thick Shell Finite Elements”, Intemational Journal for Numerical Methods in Engineering, Vol. 3, pp. 545-586. Peng X. and Crisfield M.A., 1992, “A Consistent Co-Rotational Formulation for Shells Using the Constant Stress/Constant Moment Triangle”, International Journal of Engineering Science, Vol. 35, pp 1829-1847. Piskunov V.G., Verijenko V.E., Adali S. and Summers E.B., 1993, “A Higher-Order Theory for the Analysis of Laminated Plates and Shells with Shear and Normal Deformation”, International Journal of Engineering Science, Vol. 31, N o. 6, pp. 967-988. Piskunov V.G., Verijenko V.E., Adali S. and Tabakov RY., 1994, “Transverse Shear and Normal Deformation Higher-Order Theory for the Solution of Dynamic Problems of Laminated Plates and Shells”, International Journal of Solids and Structures, Vol. 31, No. 24, pp. 3345-3374. Prathap G. and Viswanath S., 1983, “An Optimally Integrated 4-Noded Quadrilateral Plate Bending Element”, International Journal for Numerical Methods in Engineering, Vol. 19, PP. 831-840. Prathap G., 1985, “The Poor Bending Response of the Four-Node Plane Stress Quadrilateral”, International Joumal for Numerical Methods in Engineering, Vol. 21, pp. 825-835. Prathap G. and Somashekar BR. 1988, “Field- and Edge-Consistency Synthesis of a 4-Noded Quadrilateral Plate Bending Element”, International Journal for Numerical Methods in Engineering, Vol 26, pp. 1693-1708. Prathap G., 1994, “The Displacement-Type Finite Element Approach - From Art to Science”, Progress in Aerospace Science, Vol. 10, pp. 295-405. Ramm E. and Matzenmiller A., 1986, “Large Deformation Shell Analysis Based on Degenerate Concept”, in Finite Element Methods for Plate and Shell Structures, Volume I .' Element Technology, T.J.R. Hughes and E. Hinton, eds., Pineridge Press Int., Swansea, UK, pp. 365-393. Reddy J .N., 1980, “A Penalty Plate—Bending Element for the Analysis of Laminated Anisotropic Composite Plates”, International Journal for Numerical Methods in Engineering, Vol. 15, pp. 1187-1206. 168 Reddy J .N., 1984, “A Simple Higher-Order Theory for Laminated Composite Plates”, Journal of Applied. Mechanics, Vol. 51, pp. 745-752. Reddy J.N., 1987, “A Generalization of Two-Dimensional Theories of Laminated Composite Plates”, Communications in Applied Numerical Methods, Vol. 3, pp. 173-180. Reddy J.N., 1993, An Introduction to the Finite Element Method, Second Edition, McGraw-Hill, Inc., New York. Reddy J.N. and Robbins DH, 1994, “Theories and Computational Models for Composite Laminates”, Applied Mechanics Review, Vol. 47, No. 6, Part 1, pp. 147-169. Robbins D.H., Reddy J.N. and Krishna-Murty A.V., 1991, “On the Modeling of Delamination in Thick Composites”, in Enhancing Analysis Techniques For Composite Materials, NDE-Vol. 10, Winter ASME Conference 1991, Eds. Schwer L., Reddy J .N . and Mal A., pp. 133-149. Robbins DH. and Reddy J .N., 1993, “Modeling of Thick Composites Using a Layerwise Laminate Theory”, International Journal for Numerical Methods in Engineering, Vol. 36, pp. 665-677. Rossow M., 1977, “Efficient Co Finite Element Solution of Simply-Supported Plates of Polygonal Shape”, Journal of Applied Mechanics, ASME, Vol. 44, pp. 347-349. Simo J.C. and Hughes T.J.R., 1986, “On the Variational Foundations of Assumed Strain Methods”, Journal of Applied Mechanics, Vol. 53, pp. 51-54. Somashekar B.R., Prathap G. and Babu CR, 1987, “A Field-Consistent, Four-Noded, Laminated, Anisotropic Plate/Shell Element”, Computers and Structures, Vol. 25, No. 3, pp. 345-353. Spilker R.L. and Munir N.I., 1980, “The Hybrid-Stress Model for Thin Plates”, International Journal for Numerical Methods in Engineering, Vol. 15, pp. 1239-1260. Srinivas S., 1973, “Refined Analysis of Composite Laminates”, Journal of Sound and Vibration, Vol. 30, pp. 495-507. Sun CT. and Whitney J .M., 1973, “Theories for the Dynamic Response of Laminated Plates”, AIAA Journal, Vol. 11, pp. 178-183. Tessler A. and Dong S.B., 1981, “On a Hierarchy of Conforming Timoshenko Beam Elements”, Computers and Structures, Vol. 14, pp. 335-344. 169 Tessler A. and Hughes T.J.R., 1983, “An Improved Treatment of Transverse Shear in the Mindlin-Type Four-Node Quadrilateral Element”, Computer Methods in Applied Mechanics and Engineering, Vol. 39, pp. 311-335. Tessler A. and Hughes T.J.R., 1985, “A Three-Node Mindlin Plate Element with Improved Transverse Shear”, Computer Methods in Applied Mechanics and Engineering, Vol. 50. pp. 71-101. Tessler A., 1986, “Shear-Deformable Bending Elements with Penalty Relaxation,” in Finite Element Methods for Plate and Shell Structures, Volume 1: Element Technology, T.J.R. Hughes and E. Hinton, eds., Pineridge Press Int., Swansea, UK, pp. 266-290. Tessler A., 1991, “A Higher-Order Plate Theory with Ideal Finite Element Suitability”, Computer Methods in Applied Mechanics and Engineering, Vol. 85, pp. 183-205. Toledano A. and Murakami H., 1987, “A Composite Plate Theory for Arbitrary Laminate Configurations”, Journal of Applied Mechanics, Vol. 54, pp. 181-189. Turner M.J., Dill E.H., Martin H.C., Melosh R.J., “Large Deflection of Structures Subject to Heating and External Load”, Journal of Aeronautical Science, Vol. 27, pp. 97- 106. Verijenko V.E., 1993, “Nonlinear Analysis of Laminated Composite Plates and Shells Including the Effects of Shear and Normal Deformation”, Composite Structures, Vol. 25, pp. 173-185. Wempner G, 1969, “Finite Elements, Finite Rotations and Small Strains of Flexible Shells”, International Journal of Solids and Structures, Vol. 5, pp. 117-153. Whitney J.M., 1973, “Shear Correction Factors for Orthotropic Laminates Under Static Load”, Journal of Applied Mechanics, ASME, Vol. 40, No.1, pp. 302-304. Xavier P.B., Lee K.H., and Chew CH, 1993, “An Improved Zig-Zag Model for the Bending of Laminated Composite Shells”, Composite Structures, Vol. 26, pp. 123-138. Yang P.C., Norris CH, and Stavsky, Y., 1966, “Elastic Wave Propagation in Heterogeneous Plates”, International Journal of Solids and Structures, Vol. 2, pp. 665- 684. Yip Y.C., Averill R.C. and Tessler A., 1996, “A Simple, Efficient and Robust C0 Mindlin-Type, Four-Node Plate Element”, submitted for publication. 170 Yu G., Tzeng G.Y., Chaturvedi S., Adeli H. and Zhang S.Q., 1995, “A Finite Element Approach to Global-Local Modeling in Composite Laminate Analysis”, Computers and Structures, Vol. 57, No. 6, pp. 1035-1044. Zaghloul S.A. and Kennedy J.B., 1975, “Nonlinear Behavior of Symmetrically Laminated Plates”, Journal of Applied Mechanics, Vol. 42, pp. 234-236. Zienkiewicz O.C., Taylor R.C. and Too J .M., 1971, “Reduced Integration Techniques in General Analysis of Plates and Shells”, International Journal for Numerical Methods in Engineering, Vol. 3, pp. 275-290. Zienkiewicz O.C., Owen D.R.J., and Lee K.N., 1974, “Least Square Finite Element for Elasto-Static Problems. Use of ‘Reduced’ Integration”, International Journal for Numerical Methods in Engineering, Vol. 8, pp. 341-358. Zienkiewicz CC. and Xu 2., 1993, “Linked Interpolation for Reissner-Mindlin Plate Elements: Part I - A Simple Quadrilateral”, International Journal for Numerical Methods in Engineering, Vol. 36, pp. 3043-3056. Zinno R. and Barbero E.J., 1994, “A Three-Dimensional Layerwise Constant Shear Element for General Anisotropic Shell-Type Structures”, International Journal for Numerical Methods in Engineering, Vol. 37, pp. 2445-2470. Zinno R. and Barbero E.J., 1995, “Total Lagrangian Formulation for Laminated Composite Plates Analysed by Three-Dimensional Finite Elements with Two- Dimensional Kinematic Constraints”, Computers and Structures, Vol. 57, No. 3, pp. 455- 466. APPENDICES APPENDIX A APPENDIX A _ 1 13 - -9 Pi " '5' P2 = g p3 ' 5 - —a ” p4=C+b'Z-‘ p5=_§ p6=—p4 1 ~ - “ (A-l) p7 = -(_l) = _‘1 p - —9—+hg C55 p3 E 9 C(sls) C 23 21 -2 P10 = ‘= P11 = +17: c (N...) c C55 21,-1343, cl. = hard-5+6, bi = aa,+&,a—&, (A.2) d 8 (1,2;—a, e _ 2a + 21,21 1 _ a,- (1) ‘ _ J (N...) C55 C55 Nm-l N,,—l a = a[1+ Z 21,] 13 = d[2hm+ Z 1),] ‘=' “ (A3) 171 3"! at 172 N,,,—1 —hm+ch,3n+ 2 (hm—2,)(&,a—e,) i=1 —1 —ch,3n — N2 (hm —z,)(&,a—e,) i=1 Nm—l hi, - 13h; - 2, (hm — z,)(&.13 — 13,-) (A4) 1': 1 1M _Ni(hm_ (l) C i=1 (1%th + N2 (hm fish-{(i— C55" i=1 E) El 2) L ¥__/ APPENDIX B 173 APPENDIX B lntmlmndcnflntemlation This appendix is relevant to both Chapters 2 and 3 which deal with the use of interdependent interpolation functions for the MIN4-CC and LZZ3 finite element models. The details of the derivation of these interdependent interpolation functions will be shown here. The interdependent interpolation scheme enforces the linear part of a constraint of the form: 9) = O (B-l) (W’s - . along each element side. In the case of the LZZ3 model, it may be for either the top or bottom surface. The constraint of Eq. (B.1) can be associated to the vanishing of the transverse shear strain in Chapter 3 for the Mindlin plate theory. However, in the laminate theory of Chapter 4, the constraint is not equal to the vanishing of the transverse shear strain. It does not have any physical meaning. Its sole purpose is to introduce a new variable, 0,, or 9), so that the new laminate theory becomes a C0 instead of a C l theory. The discussion that follows is valid for both the four-node Mindlin as well as the top and bottom surfaces of the eight-node LZZ3 finite element models. The subscript on that usually differentiates the two surfaces in the case of the LZZ3 model will be omitted here for clarity: I- 'I 174 co .rn 3\\‘ /£r— “4* V I3 5 n ll b34 7n " 9 b23 ‘5 9 12 a ——>> n 8 bu ’4 Y “41 I “12 _ . #x ————>—> 0). “a" ‘i'J‘r bi; = ””1 Figure B.l Quadrilateral element coordinate description The xy coordinate and the 5,1] coordinate systems are refer to as the “global” and “local” coordinate systems, respectively. Independent bi-linear rotations and a Serendipity (eight-node) transverse deflection are initially assumed, i.e., 8 W = 2Nl(xry)wi i=1 4 4 mo.» = Zia-(awe..- e,(x,y> = Zia-(awe,- (13.2) i=1 i=1 4 4 x = 2Pi(§,n)xi y = 2P1(§r ")Y, i=1 i=1 where Pi(§, n) , and N ,.(fi, 11) are respectively bi-linear and Serendipity shape functions, and 0 By, and w,» denote nodal degrees-of-freedom. Because of mapping of the global xi’ coordinates x and y into the local §, 1] coordinates as depicted in Eq. (B.2), then Pi(§’ Tl)=Pi(xaPa(§, 11): yaPa(§’ 11)) and N,(§.n)=1iI.-(xaPa(§,n),yaPa(§,n)) 175 The bilinear functions are given by: P.- = ,1,“ + §.§)(1 mm) (13.3) and the biquadratic serendipity functions by: N, = (1,10+§,lé)(lHum—(1-§2)(l+n,~11)-(1—112)(1+§.~§)]§3n.-2 (3.4) 1 1 +§(1—§2)(1 +nm)(1—§3)n?+§(l—n2)(l +§,«§)(1—n?)§3 where (5.11) <—: {—(l, 1)] . For side 1-2, let s denote the coordinate running along an element edge as shown in Figure B]. The constraint along this side (11 = —1) is: 9) ’s— s |n=—1 (3.5) (w Because of the anisoparametn'c interpolation, different degrees-of—freedom will appear at the corners and mid-side nodes (see for example Figure 29). To achieve a uniform four- node configuration, the four mid-edge w degrees-of-freedom are then condensed out by the use of four differential edge constraints: (mm—93’s”g Mil = 0 (s6 [0, 1,1) (3.6) where [k is the edge length. The differential constraint equation is written for each element edge. Each of these constraints is then explicitly solved for wk + 4 and back- substituted into Eq. (B.2). So for the edge 1-2, the differential constraint will be given by: (W’ss—esrs)ln =—1 = 0 (B7) Computing the following edge derivatives: 176 §> I — Wl‘n = —l = N‘wiln = _l gut-nu»,+§§<§+1)w2+(1-§2)w5 ,5={(.,..)..G..).-2.W.}i l = (wl + W2 — 2w5)—.§ 3 (B8) E) n=-l A W’SS n=—1 where ds = jsl lag and js is the edge jacobian. n : Now, the tangential edge rotation, also known as covariant rotation when normalized by its length (refer to Figure BI) is: 9 = 0 n (B9) and will be approximated using bilinear functions similar to the Cartesian rotations, 0,C and 0),. The derivatives of the tangential edge rotation along side 1-2 will be: e, = s|n=_1 = Piesi|n=_l 1 3.10 = §[(1-§)951+(1+§)9521 ( ) A l 1 63’s = §[052_esl]7 S n = -1 Substituting Eqs. (B8) and (BIG) into Eq. (B.7), the mid-side w5 degree-of-freedom is found to be: 1 1 . where j12 = js 1 is used to denote the line jacobian corresponding to edge 1-2. The n = - line jacobian can be obtained in the following manner: 177 i = x|11=-l = Pixi|n=—l 3) = y|n=_l : Piyi|n=—l i=lIx-x] y=l[y—y] ’s 2 2 1 ’s 2 2 l . _ J» 2 . 2 J,|n=_l — (x.,) +(y.,) (3.12) x —x 2 — 2 =12 ., 1» y) 2 2 =11/2 In general, the subscripts is one of the variables E, or n . Curve s generally corresponds to one of the edges of an element of reference on which parameter s is defined. The other mid-side wk + 4 degrees-of—freedom can be obtained by cyclic changes of indices 1 -—> 2 —-> 3 —> 4, etc. from Eq. (B.l 1) and can be easily verified to be: 1 1 . W6 = §(W2+W3)‘Z(9s3“952)123 1 l . l l . “’8 = §(W4+W1)‘g(951‘954)141 wk + 4 in terms of covariant rotations is readily obtained from Eqs. (B.1 1, B.l3): 1 1 l 1 1 1 (3.14) W7 = §(w3 + W4) — 2(6g3 — 9&1) 1 l by making the appropriate substitution of the corresponding edge transformation. Eqs. (B.ll and BB) or Eq. (B.l4) can also be cast in terms of Cartesian rotations by using simple transformation equation for Cartesian tensors: 178 9 10+l0, n nx x ny y (3.15) 85 = —ln).04+lm,0y where lmr = cos(n, x) = dy/ds and I"), = cos(n, y) = —dx/ds are direction cosines of the outward normal of the boundary. The following expression for the wk 4, 4 degrees- of-freedom are then obtained: l1(2cz-xl)—§1 i= 1 182 Nm—l T15 = 2 b4i i=1 M.,—1 1+ 2 a4,- i=1 "I :3 II NM-l -(1) «(1) -(1) 1,, = Q45 — 2 (hue... —a,,-Q4s) i=1 Nm—l “(1) AU) *(1) Tl,” = -Q55 - 2 (al.-Q55 -b1,-Q45) i=1 Nm—l T25: 1+ 2 (141 i=1 Nm-l T27 = 2 C41 i=1 Nm-l -(l) «(1) -(1) T29 = Q44 — 2 (due... -c,.-C4s) i=1 Nm-l -(l) -(l) .(1) T2,“ = Q45 - Z (61.-Q55 -d1,-Q45) i=1 T31 = ‘1 Nm—l T35 = 21v.“ 2 (Zn/{20% i=1 NM—l 737 = 2 (W.,—10041 i=1 NM—l “(1) «(1) T39 = _ZNMQ44 - Z (ZNm—ZixduQM - i=1 T16 = _T15 T18 = ‘T17 A(Nm) T1, 10 = "Q45 A(Nm) T1,12 = QSS T26 = ‘725 T = —T 28 27 (C3) A(N,,.) T210 = Q44 AN...) T112 = ‘Q45 T32 = 1 Nm-l T36 = ‘ Z (sz—zgd“ i=1 T38 = ‘T37 “(1) €11Q45) 183 N,,-1 (1) (I) “(1) T3,“ = ZNMQ45 — 2 (ZN... --Z )(611st -d11Q45) i=1 —1 T45 = Ni (ZNM- 1’41 T46 = ‘T45 i= 1 ~,,_1 ~,.-1 T47 = ZNm+ 2 (ZNM‘Z1W41 T48 = ' 2 (ZNA-Zilam' i=lN1=1 (C.3) ..(1) am (1) T49 = ZN" Q45 - N": (ZN... -Z)(b11Q44-a11-Q45) i=1 1 1 ..(1) (1) ( ) T411 = _ZN». Q55 - N2 (ZN... -Z1)(011Q55- b11'Q45) i=1 (1) CJJJ) all other T11 =0 and Q45 = A” k The displacement field of Eq. (62) after introducing the interpolation functions (Eq. (B.1)) can be written as: u = aJJ’P .(x. y)f,"’(z) u = 1190111,, y)\PJ"’(z) (C4) (w1ai’1(x.y) + amt/.111 y) + 9,..N..-(x. y))M.,(z) k II where indices 01 = 1, 2; i = 1, ...,4; k = 1, 2, ..., 7 and id = 4(a— l) +i .Indices on identify whether the degree-of-freedom is at the top or bottom, k indicates the degree- of-freedom and ia, the local element node number. u(.k ) w- 0 ia’ 111’ xia and eyia are the nodal degrees-of—freedom. Note that wia = 1253’: 1212513 where 813 is the Kronecker delta. Similarly, we can write 9x101 = 1753-: 12MB“ and Gym = 4(5) 1.21105]:5 “id: The shape functions, (by), ‘1’? and M a are defined as follows: 184 oJ," = 1+rJ," of,” = o o?” = —z+I“J,"’ all other of) = Ff,“ 1112) = 1.132) A w v I - 0 ‘1’?) = —z+A(15) all other ‘1’? = A3) M1 = (1-1;) ((2.5) ZN" Z M _ _ 2 ZN," M3 _ 4(1 “5.) ZN... ZN... where N,,,—1 (k) _ 2 3 . Fa - rum:Z +’21sz + Z (Z—zi)(rl,kad2i+r2,kad3i+r3,kac2i+r4,kac3i)’ i=1 NM—l 2 3 b b . a ’3,kaZ +’4,ka '*' 2 (Z—zi)(rl,ka 2i+r2,ka 31+r3,kaa2i+r4,kaa3i)’ i=1 and koc = 2(k—1)+a;witha =1,2;k = l,2,...,7. > A R- V II APPENDIX D 185 APPENDIX D Won: The shape functions, &7)b21/b“) (1)112 = (¢&J)+&2)b22/b12) (1)13 = ((1)114) + ¢&5)b22/b12) “’19.? = (4’8?) + ¢&”b22/b12) ‘l’ég = (W&2)+‘P&”b“/b21) ‘Péfl = (‘P&5)+‘P&4)b“/b2]) will = (‘P&7)+‘P&6)b“/b21) W112 = (W842) + Whinbiz/bzz) W11? = (W&SJ+‘P&4)b12/b22) “11130? = (TQM‘I’h‘SJbiz/bzz) multiplied with the factors bl2 and b2] respectively. (D.1) Note that bl2 and 1221 may both be equal to zero, but (D1181 and ‘ngl are always pre- APPENDIX E 186 APPENDIX E _l_‘I',,l‘AillLl [AA 11!. 1-1 ‘nAt 111‘ as fmll ° . are Botatians It is important to find out how the second Piola-Kirchhoff stress and strain tensors are related to the Cauchy stress and strain tensors for the case of small strains/large rotations. Mattiasson et al. (1986) showed that there is a definite relation between the quantities. The proof is reproduced here for completeness. The relation between the reference position vector, X , the current position vector, x, and the displacement vector, y is given by: 5=X+u; xi=Xi+ui (13.1) The Lagrangian-Green strain, E can be expressed as: l l where the comma denotes differentiation with respect to the initial configuration, and 8,]. is the linear Lagrangian strain tensor. Consider a small part of a continuum in a fixed Cartesian coordinate system, X ., i = l, 2, 3 , with base vectors 2,. Applying a finite rotation to a vector, d2; through an orthogonal rotation tensor, 13 , yields: dX‘ = ads; IdXI = ldfl = d5 (13.3) $ All! Here X i , i = 1, 2, 3 define another Cartesian coordinate system with base vectors, 1,- , which initially coinciding with X,- , but rotate with the body. 187 * Now, giving the body an infinitesimal deformation so that the vector (1X turns into d5: d5 = 011121 = dx:Z: (E.4) The following relation is obtained: dx + d1; = dx" + d1.“ (13.5) where d u is the relative displacement vector between two neighboring particles. Let Idxl = ds. The difference ds2 — d52 in terms of the Lagrangian strain tensor, E , is: 2 2 i i 1 is —dS = 2 124.1121,de = 2 EijdXidXJ (E.6) Here, it must be emphasized that the displacement gradient dui/de is of finite magnitude, while duI/dx; « 1 by virtue of the infinitesimal deformation assumption. This implies that the quadratic part of the strain tensor E; can be neglected. One can therefore write E-- = 81.. Noting that dXi = dx:,we find Eij = E1 l] j tj , and finally This result shows that for a case of small strains, but large rotations, the Lagrangian strain tensor components are equal to the linear Lagrangian strain components in a system co- rotating with the particle. Using the polar decomposition theorem, the deformation gradient tensor, If can be multiplicatively decomposed into a rotation tensor, 8 and a stretch tensor, Q as: ax. E = 131] where Fij = 337’- (E8) I For infinitesimal strain, the stretch tensor can be expressed as: 188 U z 1 + e4 (13.9) where e is a small number and 41:13 = 1:1. The deformation gradient can thus be approximated by: F = R (E.10) Noting that ratio of densities, p/ p z 1 in the case of infinitesimal strain, the second Piola-Kirchhoff stress tensor can then be written as: p T S = 3005") ~9-(E‘J) (E.11) z B . g . ET The component forms of the Cauchy stress tensor, 9' are: 9 = OijZiZj = 013-21; (E12) A* o o I 0 I where [1 18 a base vector in the co-rotating CarteSian coordinate system. * The matrices of tensor components 011' and 6U are thus related by the normal transformation laws: 7' t [0] = [R] [0 ][R1 (E-13) where [R] is the matrix of the components of R in the fixed system. Combining Eqs. (EU and E.13), we get: [S1 = 1R1<1R171€11R1>1R17 = [6"] (E14) The above equation shows an all important result that for small strains, the second Piola- Kirchhoff stress components are equal to the Cauchy stress components in a system co- rotating with the particle. This implies that a constitutive relation, formulated for the case of infinitesimal strain 189 and rotation, and relating Cauchy stress and engineering strain or rates of these quantities, can be used unaltered to relate Lagrangian strain and second Piola-Kirchhoff stress. For instance, constitutive relations given in the form: (E. 15) 11 n .26- f(§) or C: Q II can be written as: S = f(E) or S (5.16) 11 '0 151 APPENDIX F 190 APPENDIX F 1 ’ v ’ i The finite element model of Eq. (138) when written in the expanded form is: 4 . 11011 1K‘211K1311K141 11051 1K‘611K'711K181 1.111 {f‘} 1K2'11K221 1K2311K241 1K2511K2611K271 [W] (.121 {f2} 111/311111321 1x33] 1x34] 1x351 1106111671 1108] M3} {f3} 1161111691 1K431 1x44] 1x45] 1x46] 1x47] [W]. W} ,=. {f4} , (1.,) 1K“) 1x5211K531 1x541 1K5511K5611K5711K581 {451 {f5} [KM] [Km] 1K6311K641 1K6511K661 1K6711K681 {d6} {f6} 1x711 1x72] 1x73] 1x74] 1W] 1x76] 1W] [W] {cm {f7} _1K8'11K8211K8311K8411K8511K8611K8711K881 {d8} J {f8} 1 where 1dll=1u1l§ {d2}={V‘-}; 1231:1111; {d4}={9,1}; . l= 1,2,...8 (F.2) 11:51:19,}; 1.1614311; {d7}={T,1}; {d8}={W3}; A special note is that {W3} will be condensed out at the element level since it is a nodeless variable. Expressing Eq. (F. 1) in indicial notation, we obtain for the ith iteration: 01,13 = l,2,...,8 (Kfl’dff = (F‘ff , k = 1,2, ...8 (if a, 13:8) (R3) ” 1 (de=8) Here summation of k, B is implied. The tangent stiffness, [K T] of Eq. (140) can then be obtained from Eq. (F.3) by: 191 a a (17 i 0111311=1.2,...,8 T1514“: [(87.1% Kim 1 . F4 ( 1k) (JmadP+ 3d? 61",] j,k,m ={1.2,...8(ifa,[3,y¢8) (.) l (lfa,B,‘Y=8) The stiffness components of [K] and [K T] are obtained from the left hand side of the virtual work equation of Eq. (126) while using the constitutive relations, Eq. (129), to relate the second Piola-Kirchhoff stresses and the Green’s strain for large rotations, small strains. 1 5,615,. dVO = j €415,512}. dVo V V 0 0 (F5) i,j,k =1,2,...,6 Here, for convenience, contracted notation has been used instead of the normal tensorial notation, S] = Sxx; 52 = Syy; S3 = Szz; S4 = Syz; S5 = 52:2; S6 = Sxy; (F6) E1 = E”; E2 = E”; E3 = E22; E4 = 2E”; E5 = 2E“; E6 = 2Exy and Cik are the elastic coefficients defined in Eq. (129). The direct stiffness, [K] is therefore obtained by adding the individual stiffness contributions of each “strain” product from Eq. (ES). This is represented mathematically as: 13 6 6 B (1 (1 Kij = 2 Z Cmnlkij lmn (F7) m = 1n = 1 where [kgfihn are defined in Appendix G and Cm" are the elastic coefficients. The tangent stiffness, [K T] , is obtained using a similar approach and are listed in Appendix H. APPENDIX G 192 APPENDIX G i ' n ' rm n st in To derive the direct stiffness matrix, it is easier to look at the individual strain products of Eq. (F.5). Using Eqs. (123, 127-129, 134) and simplifying, we obtain, for example: an (I) 1 13,85, dVo = 81.“, (j p, ,p 1., @9811!) av)”, V. V. 6“” 1P ¢P 5 N 5 N 5 M + “it: 15 1x a1 13x 13+ xj,X 14+ yj,X 15] 13(“2.x) V . 0 + Ma[Pi, 1151.3 + ”1.111514 + Nyi, X8k5](uz, x) x P. <1) 0 ( "x 6 (GI) 1 (I) + iMfllpja X813 '1' ij’ X5l4 + NYj, X515](uz, X))] deujB i,j=1,2,...,4 k,l=1,2,...,7 01,13 = 1,2 to = 4(a—1)+i = l,2,...,8 jB = 4(B—1)+j=1,2,...,8 where 511‘ is the Kronecker delta. In general, Eq. (G. 1) can be expressed in the form: k,l=1,2,...,8 m,n =1,2,...,6 kl 1 id _ 1,2,...8(ifk¢8) JEnaEm dVo = 8dla[kia,jfl]mndjfi _ 1 (if k = 8) (6'2) V 0 '13- 1,2,...8(ifl¢8) J ' 1 (ifl=8) where nodal degrees-of-freedom, di-B are as defined in Eqs. (F2, F3) and [kg 143],” is the stiffness contribution from the strain products, E m, En. Comparing Eqs. (GI) and 193 (6.2), we obtain: kl 1km. 113111 = 1P1, XPj, xq’th‘I’él) dV V 1 + J [2P1 X¢J(1k)[Pj' ,5,3 + ij,x514 + N”. X515]Mfl(uz’ X) (6.3) '1’ [P1, x5k3 "' in, 2151.4 + Nyi, x5k51Ma(“z, X) l x (P13 xd’él) + 5MB[PJE x813 + ij, X514 + Nyj’ X515](uz, 20)] dV where the first and second integrands represent the linear and non-linear contributions to the overall stiffness matrix. The other stiffness components can be derived in a similar manner and they are as follows: kl ”101113112 = [P1, XPj, rq’lkaél) dV v 0 1 + Hip,- xcpgcqu’ r513 +ij,r514+Nyj, 1'5151M13(“z, r) v, ' (G.4) ‘1' [P1321513 + Na, x5114 ‘1' Nyi, X5k5]Ma(uz, x) 1 x (Pi. 4311141),» §M5[Pj’ ,8,3 + 1v,r 1., 1514 + N”, ,5,5](uz, 0)] dV kl — — — 1km, 113113 = [P1, xd’thMB, 2119,53 + ”11514 + ””5151 (W V0 ‘1' 1 [M 01M 13, 2“” i, x5k3 + ”111,115“ + Nyi, x5k5] (G.5) v 0 x [131.63 + 71745,, + Nyj515](uz, X)] dV Notice that the consistent normal strain field, 8 Eq. (125), has been utilized and the 22’ following interpolation functions have been introduced: 194 N”. = [(all)iN§i+(a21)iNni] _ __ _ (no sum on i) (G.6) N”. = [(a12)l.N§i + (a22)iNm.] and definitions of fig? Nni are given in Eq. (115). k8 ”1011113 = 1P1,x¢¢(zk)M3,z (IV V" ((3.7) + I[MaM3,Z[Pi, x513 + ”1.121514 'J' Nyi, x5115] (“2. x” dV v kl 1km, 1.5116 = 1 Pi, xogcnpj, yogi) +Pj, xwgn dV V 0 1 + l[§MB(Pi. X911“ '1' Malpi, x5k3 + in, x5k4 + Nyi, X5k5](uz, x» v 0 (G8) x (1P1, x513 "' ij, x514 '*' Nyj, x515](“z, y) + [P1, 1'513 ‘1' Nu, 1'514 + Nyj, r515](“z,x)) ‘1‘ Ma[Pi, x5k3 'J’ in, x514 + Nyi, x5ksllpj, 1%” "' P1, x‘i’é’Wz. x)] dV kl [kia, 113121 = 1P1, YPj,X\¥t(1k)¢él) ‘1‘, v 0 l + I [51”, 3111mm, + ~.,-,.6..+ Nyj,x515]MB(“z,x) V . 0 (G.9) + Malpi, r5k3 + in, r514 + Nyi, 1'56““; 1’) 1 x (Pi. xd’fsl) + §[Pj, x513 + ij’ x514 + Nyj, X5151“; 19)] dV 195 kl ”101113122 = JPi,YPj,Y\P1(1k)\Pél) dV v 0 l V, ' (G.10) + [P1, r513 + er,y514 + Nyi, y5kslMa(“z. r) 1 x (PI, 31,!) + 5143113., ,5,, + N, 1., ,8,4 + N”, Y815](uz, ,0] dV kl — — — 1km, 113123 = [P1, Y‘Pg‘)MB, 21])1513 + ”11514 + ””5151 (W V0 1’ 11541111111215.1513 + Na, r514 + Nyi, Y8k5] (6'1 1) V0 x [131.13,3 + N285” + Nyjagnuz, Y)] dV k8 ”101.1123 = [P1, yq’hk)M3,z (IV V" ((3.12) + IlMaM3,Z[Pi, r5113 + in, r5k4 '1' Nyi, Y5k5] (“a Y)] dV v 0 kl [1,4 1.4126 = (P1, figure, yogupj, xwgn dV V 1 + 1 [51111391, 1"?th + Ma[Pi, Y5k3 + in, Y5k4 + Nyi, Y8k5](uz, 1'» v0 ((3.13) X ([Pj, x513 '1' Nu, x514 ‘1' 1vyj,X515](uz. r) '1' [P1, Y513 + ”1.1, r514 + Nyj, Y515](“z,x)) '1‘ Malpi, r513 "’ in, 1'514 ‘J' Nyi, Y6k5][Pj, y‘pél) + P}, x‘l’é”1(uz, 1'” “JV kl - — — [k,a, 1.413, = j P j, chgwa, 2110,61, + 11114.5,4 + Nyifiks] dV V 0 1 _ _ " (3.14 '1'- J[§MBMQ’Z[Pl5k3+NXi5k4+Ny16k5] ( ) x [P1311613 + Nx1,x514 + Nyj, x515](“z, X)] dV 196 81 [k1,113131 = [Pfix‘bénMade v 0 l + l[§MBM3.Z[Pj, X513 +ij,x514+Ny,-,x5151(uz,11)] (W Va kl — — - [kia,jB]32 = [P1, Y‘PénMazlprsu + N115“ ‘J' NyiSkS] dV v 0 1 _ _ _ 81 [k1.1'13]32 = [1943501431 dV V 0 1 + J‘[§MBM3, Z[PJ'.Y813 + ij, 1’514 + Nyj, 13151042, y)] 611’ V0 kl — — — 1km, 113133 = [Mot 2MB, Z[P1‘5k3 + ”11514 + NyiakS] V0 x [131.8,3 + ij8,4 + Nyjsgj dV k8 — — — ”101.1133 = 1M3,2Ma, 21131-5113 +Nx15k4+Ny18k51 dV v 0 81 - — — [k1, 11,133 = 1 M3, 2MB, 291513 + 111145,4 + Nyjag] dV V 0 ”$131133 = [(M3,z)2 (IV V 0 (G. 15) (G. 16) (G. 17) (G.18) (G. 19) (G20) (G.21) 197 kl - — — [km ”3],, == A'111,Z[Pj,l’(bél)+ P j, x‘I'éJJHP1513 + N11511: + N 1.815] dV V y 0 1 _ _ _ V x([Pj,x513+ij,x514+Nyj,x515](“z.r) 81 [1,, 1.13136 = (M31110, yogi) +PJ.’ MU] (W V 1 G23 + [ [5M3 ZMB([PJ.’ Xian-AN“; X5,4+Ny 1., 115151041, ,) J l v, ’ + [P1, r513 +ij, y514 +Nyj, Y515](“z.x))] (W The consistent form of the transverse shear strains, Eqs. (118, 119) must be used to evaluate the stiffness components for the non-linear model. For convenience, the transverse shear strains will be expressed in a concise form as: — (k) «(1‘), - _ (k) ,.(k) 7.. = inaum, y. — Lyman, ((3.24) where the interpolation functions, Liz, LS; are readily obtained from Eqs. (118, 119) and are given as: (l) (1) (1) ina = [(011)1P1¢§a,zb11+(“21)1P1¢ b12] “(1.2 (2) (1) (1) ina = [(a12)1pi¢§a,zbll +(022)1Pi¢na.zb12] (3) ina = [Pvgbn +P1’nb12] (4) _ (2) (2) ina " [(011)1(Pi¢§a,z + Néi’fiMalbll + (021)1(Pi¢na,z + Nni’nMalbIZ] (6'25) (5) (2) (2) ina = [(a12)1'(Pi(b§a, z + NfiirgMalbll '1' (“22)1(Pi(bna,z + Nni'nMa)b12] (6) (3) (3) ina = [(011)11’1‘1’501, zbll + (“21)1Piq’na. 21,12] 7 3 3 . L5,“), = [(012)1P1¢(§d,zb11 + (022)1P1‘Pg‘1izb12] (no sum on i) 198 (l) (1) (1) Lyia = [(011)1P1‘Pga, zb21 +(021)1P "1' [’22] I “(1.2 (2) (1) (1) Lyia = [(012)1P1‘Pga.zb21+(022)1P1‘Pna,zb22] (3) Lyia = wia[Pi’§b2l +Pi’nb22] 4 2 2 Lg“; = [(011)1(P1‘PéoZz+N§11§Malb21+(0201031161,;+Nn1’nMa)b22] (G26) 5 2 2 L( ) = [(012)1(Piqléa),z+N§12§Ma)b21 +(022)1(Pi‘y( ) +Nni’nMalb22] M m. z 6 3 3 L( ) = [(011)1P1‘1’éo1),1b21+(021)-P"¥( ) b22] yia 1 1 na,z L3; = [(012)1P1‘P22,1b21+(022)1P1‘1’$.3J,1b221 (no sum on 1) 1ka 1.5144 = I 4243),, dV (G.27) Va [kg 1.11145 = 14’:ngng (6.28) Vo [kt-1,111.. = [lei-‘11.]... (6.29) [£341,155 = [42491, dV ((3.30) V kl [kia,jB]61 = IP13 X¢Bl)[Pi, y‘Pg‘) + Pi’x‘llakll dV V 1 + HIP!) JADE) + 5541311)}, X513 + ij, X514 ‘1' Nyj, 115151041, 11)) V0 x Ma([Pi, x5113 "' ”111,115“ + Nyi, 115/15““; 1') ‘1' [P1, Y5k3 "' in, r5114 'J' Nyi, r5k5](“z, x» 1 + iMfllpj’ X5” + N“? X5“ + Nyj. X5k5][Pi. Yq’hk) ‘1' P1, x‘Péf’Muz, x)] (W (G.31) 199 kl ”(01.113162 = 1P1. ”10m. y¢ékJ+P1,x‘I’&"’l W V kl (kid, jB]63 1 + j [(10]; $104,511,315, ,8,, + ij, ,8,4 + NH, Y815](uz, ,)) V0 x Ma([Pi, x5113 + ”111,115“ + Nyi, x5ks](“z, y) "’ [P13115113 + Na, 1’514 + Nyi, r5k5](“z, 21)) (G.32) 1 + §M111P1 1513 + N 514 + N 5115191. Yog") + P1,X‘P&k)](uz, V] W xj,Y yj,Y = 1 MB, Z[P,., ,(pgc) + P, X‘I’&’<)}[P 113,, + Fig.8“ + Ivy-515] W V 0 + [[MGMB, 21131-513 + NUS” + Nyjfifi] Vo X (1P1, x5k3 + ”111,115“ + Nyi, x5k5](“z, y) + “’1, 1’513 ‘J‘ N111, y5k4 ‘1' Nyi, r5ksl(“z.x))] (IV (6.33) k8 [k1a,1]63 = 1M3,2[P1, Y1J,") +P1.X‘P((xk)] dV v 0 + l[M3. ZMa([P1, 1151.3 ‘1’ Nx1,x514 + Nyi, 11815101,, y) (G34) V0 + [P1, r5113 + Nx1,y5k4 + Nyi, 75115105, x) )] dV 200 kl [k,a, 1.5166 = 1 [P1, yoga + Pi, fighupj, yog) + P 1., MO] (1V V 1 + j [5143(1104 Yg<> + P, X81510] V 0 + Ma([Pi, x5113 '1' ”111,115“ "’ Nyi, x5151“; r) + [P1, 1151.3 + in, r5k4 + Nyi, r515](“z,x)) (6.35) X (1P1, x513 + ij, X514 ‘1’ Nyj, x515](“z, r) + [P131313 '1' ~11}, r5114 '1’ Nyj, 1'55““; 11)) + [P13 y<1>é’) +Pj, X91110] X Ma([Pi, x5113 + in, x5114 + Nyi, 21515105, 1') + [P1, Y8k3 + N111, r514 + Nyi, y5115](“z,x))] dV all other [ka 1111.... = 0 (G36) APPENDIX H 201 APPENDIXH 11 |S|°m Mlllv K, E 11] The tangent stiffness components are obtained via Eq. (E4). They are as follows: T kl kl J kiijJH = [kimjbln l + l[‘2'Pi. X¢t(1k)[Pj,X513 +ij,X814+Nyj, x5151M13(uz,x) v 0 + (MaMB[P1, X513 ‘1' ”xi, x514 '1' Nyi, x5115” X [P1,x513 +ij,X514+Nyj,X615][(uz,X)2+ (“x,x)]] dV (H-l) i,j =1,2,...,4 k,l = l,2,...,7 01, = 1,2 ion = 4(01—1)+i = l,2,...,8 113 = 4(B—1)+j = l,2,...,8 T kl kl 1 kia,jB]12 = [1908113112 1 "’ [[5]: X¢hk)[Pj, 1513 +ij, y514+Nyj, r5151M13(“z.y) v, ’ + (MaMBIP1,x51¢3 + in,X8k4+Nyi,X8k5]) (H2) 1 X (5(1)); 1513 ‘1' N11}, 1'514 'J’ Nyj, r515)(“z,x)(“z. r) 1 + (P 1., ,5,3 + ij.’ X51, + Nyj’ X515)|:(uy, ,) + 5(114 Y)le dV r kl kl [ kia,jl3]l3 = [kia,jB]13 "‘ [[MaMBIP1,x51¢3 +Nx1,x5k4+Ny1, x515] (H-3) v 0 x [P13 111513 + ij’ x514 + Nyj’ x515](“z, 2)] (1V 202 r kl kl [ kia,j13]16 = [kidjblm 1 + I [5(P1,xg‘)+ Ma[Pi’ x513 + ”111,115“ + N111, 115/15““; x)) V0 xMfl([PJ-’X5,3 +ij,x514+Nyj,x515](“z, ,1) (HA) + “DJ. y613 + ij’ y6,4 + N”, y615](uz, x)) + MaM31P1,x513 + N 5,,4 + M44515] X [Pj,x513 +ij, 11814 + Nyj,x515][ux, 1’ + uy,X + (uz,X)(uz, 1’)“ dV xi, 1: r kl kl [ km, 10121 = 1km, 111121 1 ‘1' I [2P1- ,‘1’1JJIP1,x513 + ”111', x514 + Nyj, 1151515411“; x) v, ’ + M 61M (31” i, ran ‘J’ N xi, y5k4 + Nyi, r5k5] (H5) 1 x (5(Pj, X513 + ij, X514 + Nyj, x515)(“z, x)(uz, y) 1 2 + (P1, 1’813 + Nun/514 '1' Nyj, y515)[(ux, X) + i012, X) ])] dV T kl kl [ kia,jB]22 = [kidjfllzz 1 V0 ’ (H.6) '1' (MaMB[P1, r5113 ‘1' ”xi, 1151.4 ‘1' Nyi, r5115” x [P13Y813 + ij, ,6“ + N”; y515][(uz, y)2 + (uy, Y)]] (W T kl kl [ kia, 113123 = [kid 113123 + I[MaMB[Pi, Y5k3 + Nx1,r5k4 + Nyi, Y5k5] (H-7) v 0 X [P1, 1'513 + N11}, 1514 + Nyj, 1515““; 2)] (IV 203 T kl kl [ kim113126 = [hm-11121, 1 + 1 [5031, 1"?th "’ Malpi, Y5k3 + ”111, 1’514 + Nyi, y5k5](uz, y)) V0 x MB([PJ-’ ,6,3 + ”11111514 + N”, X65101, y) (H8) '1’ [P1, 7513 + ij, y514 “J” Nyj, r515](“z. 11)) + MaMfllpi, 1'513 + Na, r514 + 1Vyi,1’5k5] X 1P); 1'513 "’ ij, r514 + Nyj, 1251511141, r + “y, x + (“8 XX“; 1')“ CW T kl kl [ kia,jB]3l = [kidj13131 1 — _. _ + 1 [5M ,1 4MBIP1513 + IVA-51.4 + N,.5.sl (H9) V0 X [P1, x513 4' ij, 11514 + Nyj, x515](“z, x” dV where NJ", N . are defined in Eq. (G.5). yr T 81 81 [ k1,}13131 = “1.113131 1 (H.10) v, ’ T kl kl [ kra, 18132 = “‘11:, 16132 1 _ _ _ + 1 [5M 1, 21411191513 + ”..,-514 + ”,..-5.51 (11.11) v, ' x [P12 1513 + ij’ r514 + Nyj, 1'55““; Y)] dV T 81 81 [ k1,;13132 = ”1.113132 (H.12) l V, ’ 204 T kl kl [ kia,jB]36 = ”101.1356 1 _ _ _ + J [55401. ZMB[Pi5k3 "’ Nx15k4 + Nyi6k5] (H. 13) X ([Pj, x513 "' ij, x514 + Nyj, x515](“z, y) + [P]; y513 +ij, y514 +Nyj, y515](“z,x))] (W T 81 81 [ kl,jB]36 = [k1,jbl36 l + I [5M3 ZMWPJ: x513 + ij. x514 + Nyj.x5151(uz, y) (11.14) v, ’ + [P111513 + ij’ YE),4 + Nyj’ ,615](uz, x) )] dV T kl kl [ kia,jB]6l = [kimjfllm l + I [EMOIMBIPJL X513 + ij, X514 + Nyj, x515] Va X ([P 1, 115113 + ”111,115“ + Nyi, 7515]“; x)2 + [P8115113 + ”n, x5114 + Nyi, x5k5](“z, x)(“z, y)) (H.15) + MaMfisz, x)2 + (ux’ 11)] X ([Pi, x5113 + ”n, x5k4 "' Nyi, X8k5][Pj, Y5k3 + ij, y5k4 '*' Nyj, 115115] + ”’1, r5113 + ”111,115“ "‘ Nyi, YakSHPj, x513 + ij, X514 + Nyj, x5151) l + 5M 3“” j, x5113 + N 5114 + N 51:51“D 1, chg‘) + Pi, #9)]("5 X)] “W xj,X yj,X 205 T kl kl [ km. 113162 = [kimjfiltsz l + I [EMGMBIPJL Y513 + ij, Y514 + Nyj, y515] V0 X (IP13 x5113 4' in, x5114 + Nyi, x5k5](“z, y)2 "' [P1, y5k3 "’ in, y5k4 + Nyi, 751151“; x)(“z, y)) (H.16) 1 + MGMB[§(uz’ ,)2 + (11% g] X (“’13 x5113 + ”111,115“ + Nyi, x5k5][Pj, 115113 + N + [P1, Y8k3 + in, y5k4 + N 1 x j, Y5k4 + ”N. Y5k5] yi, Y8k5][Pj, x513 + ij, x514 + Nyj, x5151) 51:4 + N 5115]“); Y‘I’g‘) +Pi,X‘P((xk)](uz, fl] dV xj,Y yj,Y T kl kl [ kia,jB]63 = [kia,j[3]63 + [MaM (uz ) Jo B ’2 (H.17) X (“D1, x5113 4' in, x5114 "‘ Nyi, X6k5][Pj. Y5k3 + ij. Y5k4 + Nyj. 1’5“] + [P1, y5k3 4' in, 115184 + Nyi, YakSHPj, x513 + ij, x514 + Nyj, 115151)] 4" T kl kl [ ioz,jB]66 = [kia,jB]66 l + j [5541391. ,cbgk) + Pi, xwgn V 0 X ([Pj, x513 + ij, 111514 + Nyj, x515](“z, y) + [P]; 75113 + ij, 115114 "' Nyj, 115/15““; 11)) l + iMaMBHPj. x513 + ij. X614 + NH" X5’5](uz' Y) + [P]; Y5k3 + ij, 215114 "’ Nyj, Y5k5](“z, X)) X ([P 1, X8k3 4' in, x5114 "’ Nyi, x5k5](“z, y) 4' [P13115113 + ”11115114 + Nyi, Y5k5](“z, x)) + MaMpKux, y) + (My, x) + (uz, x)(uz, y)] X ([P 1, X513 + in, x5114 + Nyi, 115/151”}; Y5k3 + ij, y5k4 "’ N y j, Y5k5] "' [P1, Y5k3 + ”111,115“ + Nyi, 115/15]”); x513 + ij, x514 + N 5151)] (W (H.18) yj. X T ' ' kl all Other [ kt}, 181nm = [kim jBlmn (H.19) ICHIGAN STRTE UNIV. LIBRARIE W”WI1WWW”WWI“”IHIWIWINWWII 31293015671856