vw v”- ' ‘4 3).: — . . - 2'; v .A ..~A ‘ -A “.15! ’f‘ .- -,:. -... 44." .-':"':?~'::. 3;: I" "v - 22w. 3-; ,— ‘V‘. iv. -,_‘_, ~n.‘. my; 1' Mas " s ‘ I ‘ R L 2': 31' r'! N‘ "L‘t‘g’ hi” Ilillili'liiililiiiiilii'iliili'fliiifliiflfililil 31293 01712 6578 This is to certify that the dissertation entitled NONLINEAR RANDOM VIBRATION OF LAMINATED FRP PLATES USING HIGH-ORDER SHEAR THEORY presented by Joowon Kang has been accepted towards fulfillment ofthe requirements for Ph.D degree in Mneering Major professor Date 4/0791/96’ MSU is an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE we, WW 1m CICIWM.“ _4—-—— NONLINEAR RANDOM VIBRATION OF LAMINATED FRP PLATES USING HIGH-ORDER SHEAR THEORY By Joowon Kang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirement for the degree of DOCTOR OF PHILOSOPHY Department of Civil and Environmental Engineering 1998 ABSTRACT NONLINEAR RANDOM VIBRATION OF LAMINATED FRP PLATES USING HIGH-ORDER SHEAR THEORY By Joowon Kang The nonlinear response of laminated fiber reinforced plastic (FRP) plates modeled with finite elements and excited by stochastic loading is studied. The use of FRPs has been extended rapidly for structural applications in recent years due to their outstanding mechanical properties. Most FRP materials have strong anisotropic properties and exhibit significant nonlinearity in the shear stress—strain law. A high-order shear theory is used to account for the variation of strains through the thickness, since Kirchhoff and Mindlin plate theories are usually inadequate for modeling laminated FRP plates of reasonable thickness. Approximate nonlinear random vibration analysis is performed using the method of equivalent linearization to account for material nonlinearity. The solutions are obtained using an iterative approach, where linear random vibration analysis is performed in each iteration. A computer program is developed and implemented for the nonlinear random vibration analysis. Cantilevered FRP plates consisting of single-ply and three-ply plates with rectangular geometry and Boron/Epoxy Narmco 5505 FRP material are con- sidered in numerical examples. Various types of fiber orientations and stacking sequences are considered. The results show that the effect of nonlinearity on the responses becomes more pronounced as the excitation level is raised, and hence the number of iterations required for convergence also increases. The material nonlinearity significantly affects some displacement, strain, and stress responses, depending on the fiber orientations, stack- ing sequences, material properties and plate geometry. Use of the high-order shear defor- mation theory also significantly affects the responses in comparison to simpler theories such as the classical and first-order shear deformation theories. A formulation for deter- ministic dynamic analysis is also developed and some deterministic simulations are per- formed to verify the accuracy of the approximate nonlinear random vibration method. The random vibration analysis is found to be sufficiently accurate and considerably more cost effective than the use of deterministic simulations. Dedicated to my father, Sang-Tae Kang iv ACKNOWLEDGEMENTS First and foremost, I would like to express my sincere gratitude to Dr. Ronald S. Harichandran for his guidance, encouragement and support. I would like to thank you my committee members, Dr. F. Hatfield, Dr. D. Liu and Dr. D. Feldman. Their assistance in answering questions and providing materials in areas associated with my research was invaluable. Most importantly, I could not have completed this work without the support and understanding of my wife, Mee-Kyung and my children Sung-Hyun and Sung-Woo. They are the ones who have sacrificed the most during the years of my studies. Finally I would like to express sincere gratitude to my father, mother, brothers and parents-in-law. They always encouraged me to go one more step until I reached my goal. Table of Contents List of Tables 1x List of Figures.. ............................... - - - . ...... ........ . .............. x CHAPTER 1 Introduction. 1 1.1 Dynamic Analysis of FRP Plates .................................................................... 1 1.2 Objective ......................................................................................................... 2 1.3 Overview of Random Vibration ...................................................................... 3 CHAPTER 2 Literature Review...--.- - ....... . ..... - - ........ 7 2.1 Material Nonlinearity of FRPs ........................................................................ 7 2.2 Laminated Plate Theories ............................................................................. 14 2.3 Equivalent Linearization Techniques ........................................................... 19 2.4 Random Vibration Analysis of FRP Laminates ............................................ 21 CHAPTER 3 Finite Element Formulation............................. ..... . ...... 24 3.1 Introduction ................................................................................................... 24 3.2 Nonlinear Constitutive Equations ................................................................. 25 3.3 Classical Laminate Plate Theory .................................................................. 29 3.4 First—Order Shear Deformation Theory ........................................................ 31 3.5 Higher-Order Shear Deformation Theories .................................................. 33 3.5.1 Displacement Field Equations .................................................... 33 3.5.2 Strain-Displacement Relations .................................................... 36 3.5.3 Stress-Strain Relation .................................................................. 37 3.5.4 Finite Element Formulation ........................................................ 40 3.5.4.1 Extensional Deformation .......................................... 40 3.5.4.2 Shear Deformation .................................................... 42 3.5.4.3 Bending Deformation ................................................ 43 3.5.4.4 Strain-Nodal DOF Relations ..................................... 44 3.5.4.5 Linear Element Stiffness Matrix ............................... 49 3.5.5 Consistent Mass Matrix .............................................................. 52 3.6 Numerical Integration ................................................................................... 54 CHAPTER 4 Equivalent Linearization Method 58 4.1 Introduction ................................................................................................... 58 vi 4.2 Equivalent Linear Equation of Motion ......................................................... 59 4.2.1 Gaussian Approximation ............................................................ 62 4.3 Nonlinear Element Stiffness Matrices .......................................................... 64 4.4 Random Vibration Analysis .......................................................................... 71 4.5 RMS Strains and Stresses in a Laminate ...................................................... 73 4.5.1 Computation of RMS Strains ...................................................... 73 4.5.2 Computation of RMS Stresses .................................................... 75 CHAPTER 5 Numerical Results 82 5.1 Introduction ................................................................................................... 82 5.2 Model Description ........................................................................................ 82 5.2.] Selection and Modification of Nonlinear Model ........................ 82 5.3 Numerical Results for Single-Ply Plate ........................................................ 83 5.4 Numerical Results for Three-Ply Laminated Plate .................... - ................... 87 5.4.1 Symmetric stacking sequence ..................................................... 87 5.4.2 Antisymmettic Stacking Sequence ............................................. 96 5.5 Comparisons with Lower Order Theories ................................................... 102 5.5.1 Symmetric three-ply laminated plates ...................................... 103 5.5.2 Antisymmettic three-ply laminated plates ................................ 111 5.6 Effects on 2-3 Plane Isotropy Assumption ................................................. 112 5.6.1 Single-Ply Plate with 2-3 Plane Isotropy .................................. 114 CHAPTER 6 Deterministic Nonlinear Dynamic Analysis ._ man—3119 6.1 Introduction ................................................................................................. 119 6.2 Incremental Formulation of Equation of Motion ........................................ 120 6.2.1 Implicit Time Integration .......................................................... 120 6.2.2 Formation of Tangent Stiffness Matrix ..................................... 121 6.2.3 Formation of Element Load Vector .......................................... 126 6.3 Convergence Criteria .................................................................................. 127 6.3.1 Convergence Criterion Based on Displacements ...................... 128 6.3.2 Convergence Criterion Based on Out-of-balance Load ............ 129 6.3.3 Convergence Criterion Based on Internal Energy .................... 129 6.4 Procedure for Implicit Integration of Systems with Material Nonlinearity using the Newmark [3 Method ..................................................................... 130 6.5 Numerical Results ....................................................................................... 133 CHAPTER 7 Conclusions and Recommendations 137 7.1 Summary and Conclusions ......................................................................... 137 7.2 Recommendations for Future Research ...................................................... 139 vii LIST OF REFERENCES . ....................................................................... 141 APPENDIX A .......................................................................................... 147 APPENDIX B ........................................................................................... 150 viii TABLE 5.1 TABLE 5.2 TABLE 5.3 TABLE 5.4 TABLE 5.5 TABLE 5.6 TABLE 5.7 TABLE 5.8 List of Tables First five natural frequencies from linear (first iteration) and nonlinear (last iteration) analysis for the excitation intensity SD = 100,000 N 23 ............... 84 Number of iteration required for convergence ........................................... 85 First five natural frequencies from linear (first iteration) and nonlinear (last iteration) analysis for the excitation intensity SO = 100,000 st ............... 91 Number of iteration required for convergence ........................................... 91 First five natural frequencies from linear (first iteration) and nonlinear (last iteration) analysis for the excitation intensity SO = 100,000 N25 ............... 97 Number of iteration required for convergence ........................................... 98 Natural frequencies predicted by the three theories-symmetric three-ply ................................................................................................... 103 Natural frequencies predicted by the three theories-antisymmetric three-ply ................................................................................................... 111 ix Figure 1.1 Figure 1.2 Figure 1.3 Figure 3.1 Figure 3.2 Figure 3.3 Figure 3.4 Figure 3.5 Figure 3.6 Figure 3.7 Figure 3.8 Figure 3.9 Figure 3.10 Figure 3.11 Figure 3.12 Figure 3.13 Figure 3.14 Figure 3.15 List of Figures Spectral density function of white noise ...................................................... 3 Gaussian response process with three-RMS ss design ................................ 4 Overview of random vibration analysis ....................................................... 6 Global and material coordinate systems for a lamina ................................ 25 Fit of approximate shear stress-strain law ................................................. 28 Kirchhoff’s assumption .............................................................................. 30 Assumption for first-order shear theory ..................................................... 31 Assumption for higher-order shear theory ................................................. 34 Typical normal strain variation through the thickness ............................... 38 Typical transverse shear strain 713 variation through the thickness ................................................................................. 39 Typical transverse shear strain 723 variation through the thickness ................................................................................. 39 Typical 4-noded isoparametric quadrilateral element ................................ 40 Nodal displacements for extensional deformation of a plate element ........................................................................................... 41 Nodal displacements for shear deformation of a plate element ................. 42 Nodal displacements for flexural deformation of a plate element ............. 43 Geometry of laminated plate and layer numbering system ....................... 50 Cantilevered plate loaded in flexure .......................................................... 57 Variation of normalized tip displacements with a length-to-thickness ratio of the plate .................................................................................................. 57 Figure 5.1 Figure 5.2 Figure 5.3 Figure 5.4 Figure 5.5 Figure 5.6 Figure 5.7 Figure 5.8 Figure 5.9 Figure 5.10 Figure 5.11 Figure 5.12 Figure 5.13 Figure 5.14 Figure 5.15 Single-ply plate loaded in flexure .............................................................. 83 Variation of absolute r.m.s. shear strains with excitation level for single-ply .............................................................................................. 86 Variation of normalized r.m.s. displacement with excitation level for single-ply .............................................................................................. 87 Variation of normalized r.m.s. normal stresses with excitation level for single-ply .............................................................................................. 88 Variation of normalized r.m.s. shear stresses with excitation level for single-ply .............................................................................................. 88 Variation of normalized r.m.s. normal strains with excitation level for single-ply .............................................................................................. 89 Variation of normalized r.m.s. shear strains with excitation level for single-ply .............................................................................................. 89 Three-ply symmetric laminated plate loaded in flexure ............................ 90 Variation of absolute r.m.s shear strains with excitation level for symmetric three-ply ..................................................................................................... 92 Variation of normalized r.m.s displacement with excitation level for symmetric three—ply ................................................................................... 93 Variation of normalized r.m.s. normal stresses with excitation level for symmetric three-ply ................................................................................... 94 Variation of normalized r.m.s. shear stresses with excitation level for symmetric three-ply ................................................................................... 94 Variation of normalized r.m.s normal strains with excitation level for symmetric three-ply ................................................................................... 95 Variation of normalized r.m.s shear strains with excitation level for symmetric three-ply ................................................................................... 95 Three-ply antisymmetric laminated plate loaded in flexure ...................... 96 xi Figure 5.16 Figure 5.17 Figure 5.18 Figure 5.19 Figure 5.20 Figure 5.21 Figure 5.22 Figure 5.23 Figure 5.24 Figure 5.25 Figure 5.26 Figure 5.27 Figure 5.28 Figure 5.29 Variation of absolute r.m.s shear strains with excitation level for antisymmetric three-ply ............................................................................. 99 Variation of normalized r.m.s displacement with excitation level for antisymmetric three-ply ............................................................................. 99 Variation of normalized r.m.s normal stresses with excitation level for antisymmetric three-ply ........................................................................... 100 Variation of normalized r.m.s shear stresses with excitation level for antisymmetric three-ply ........................................................................... 101 Variation of normalized r.m.s normal strains with excitation level for antisymmetric three-ply ........................................................................... 101 Variation of normalized r.m.s shear strains with excitation level for antisymmetric three-ply ........................................................................... 102 Variation of absolute r.m.s. 712 shear strain with excitation level for symmetric three—ply ................................................................................. 104 Variation of absolute r.m.s. 713 shear strain with excitation level for symmetric three-ply ................................................................................. 105 Variation of absolute r.m.s. 723 shear strain with excitation level for symmetric three-ply ................................................................................. 105 Variation of normalized r.m.s. displacement with excitation level for symmetric three-ply ................................................................................. 106 Variation of normalized r.m.s. 611 with excitation level for symmetric three-ply ................................................................................ 107 Variation of normalized r.m.s. 0'22 with excitation level for symmetric three-ply ................................................................................ 107 Variation of normalized r.m.s. 1:12 with excitation level for symmetric three-ply ................................................................................ 108 Variation of normalized r.m.s. 113 with excitation level for symmetric three-ply ................................................................................ 109 xii Figure 5.30 Figure 5.31 Figure 5.32 Figure 5.33 Figure 5.34 Figure 5.35 Figure 5.36 Figure 5.37 Figure 5.38 Figure 5.39 Figure 5.40 Figure 6.1 Figure 6.2 Figure 6.3 Variation of normalized r.m.s. 1:23 with excitation level for symmetric three-ply ................................................................................ 109 Variation of RMS 1:13 through the thickness for symmetric three-ply ................................................................................ 110 Variation of RMS 123 through the thickness for symmetric three-ply ................................................................................ 110 Material direction with 2-3 plane of isotropy .......................................... 113 Variation of absolute r.m.s. shear strains with excitation level for sin gle-ply .................................................................................................. 1 15 Variation of absolute r.m.s. shear strains with excitation level for single-ply .................................................................................................. 1 15 Variation of normalized r.m.s. displacement with excitation level for single-ply .................................................................................................. 1 16 Variation of normalized r.m.s. on with excitation level for single-ply .................................................................................................. 1 17 Variation of normalized r.m.s. 022 with excitation level for single-ply .................................................................................................. 1 17 Variation of normalized r.m.s. shear stresses with excitation level for single-ply .................................................................................................. 1 18 Variation of normalized r.m.s. shear stresses with excitation level for single-ply .................................................................................................. 1 18 Convergence criteria based on displacement increments and out-of-balance loads ................................................................................ 128 Typical simulations of linear and nonlinear responses ............................ 134 Variation of nonlinear r.m.s. displacement with truncated length of initial transient part ................................................................................. 136 xiii Figure 6.4 Variation of linear and nonlinear r.m.s. displacements with number of simulations .............................................................................................. 136 xiv CHAPTER 1 Introduction 1.1 Dynamic Analysis of FRP Plates Fiber reinforced plastics (FRPs), a group of advanced composite materials consist- ing of strong and continuous fibers of different chemical origin embedded in a plastic matrix, have been developed and used for a great variety of engineering applications in aeronautical, mechanical, chemical and other industries over the past 20 years. FRPs have outstanding mechanical properties, such as high strength to weight ratio, excellent corro- sion resistance, and very good fatigue characteristics, etc, and are being increasingly used in civil engineering. Considerable research has been done to characterize the static and dynamic response of structures made of FRPs. Much of the work on dynamic responses is based on deterministic analysis. Relatively little work has been done on the response of FRP structures exposed to stochastic dynamic loads. Laminated FRP plates and shells used in aeronautical and aerospace vehicles are often subjected to random loads. The behavior of FRP structures is considerably different from that of structures made of conventional materials such as steel and aluminum, because of the many different types of constituent materials and numerous options in the fiber orientation and ply stack- ing sequence. Another noticeable feature is elastic nonlinear material response when shear loading is involved. The research described in this thesis focuses on laminated FRP plates. Early work on laminated FRP plates used classical plate theory with Kirchhoffs assumption. This the- ory yields sufficiently accurate displacements and in-plane stresses for thin plates. How- ever, since the ratio of the elastic in-plane modulus to the transverse shear modulus is very large in laminated FRP plates, the classical theory which neglects transverse shear defor- mation is inadequate even for moderately thick laminated FRP plates. Many plate theories have been proposed to include the effect of shear deformation, of which, the laminated version of the first-order shear deformation theory developed by Reissner (1945) and Mindlin (1951) is the simplest. This theory assumes a linear distribution of the in-plane normal and shear stresses over the thickness which results in nonzero transverse shear stresses, but does not produce the nonlinear transverse shear stress distribution through the plate's thickness. High-order shear deformation theories can overcome the limitations of the first-order theory by introducing additional degrees-of—freedom (DOF). The third order theory proposed by Reddy (1987) accounts not only transverse shear effects but also pro- duces a parabolic variation of the transverse shear stress through the thickness of the plate. Random vibration analysis of structural and mechanical systems has gained increasing acceptance in many engineering applications. Many structures subjected to ran- domly varying dynamic loads such as earthquakes, wind loads, ocean waves and wind tur- bulences are commonly encountered in buildings, dams, nuclear facilities, offshore structures, aircraft, etc. Unfortunately, some of these structures exhibit nonlinear behavior due to material properties or geometric configurations. The analysis of nonlinear systems under random loadings can seldom be accomplished exactly and approximate methods must be used. The method of equivalent linearization is widely used due to its robustness with even strongly nonlinear systems. In this thesis, a random vibration analysis technique for laminated FRP plates modeled with finite elements and including material nonlinearity is developed. The tech- nique is based on the method of equivalent linearization. 1.2 Objective The objective of the research is to perform random vibration analysis of thick lam- inated FRP plates exhibiting nonlinearity in the shear stress-strain relations. Shear defor- mations are to be included, and in order to include complex boundary conditions a finite element-based approach should be used. 1.3 Overview of Random Vibration Dynamic analysis of structural systems is used to predict responses under specified dynamic loads. When structural systems are exposed to stochastic dynamic loads, random vibration analysis can be used to obtain statistics of the random responses. One of basic assumptions in random vibration analysis is that the structural system is known and deter- ministic, while the random excitations have well-defined statistical properties. In many applications, the excitations are assumed to be stationary Gaussian random processes. This assumption simplifies the analysis and is yet quite reasonable in describing physical loads. For stationary random vibration analysis, the frequency domain approach is simpler and often preferred over the time domain approach. In the frequency domain approach, a stationary random process is characterized through a spectral density function (also known as the power spectrum). Many physical random loads may be approximated by white noise, which leads itself to simple solutions. White noise is defined by constant power over all frequencies (i.e., a constant spectral density function). This spectrum is shown in Fig. 1.1. More realistic spectra of physical loads can often be obtained by filtering a white noise through a linear dynamic system. S,(w) 0 0) Figure 1.1 Spectral density function of white noise Once the random excitation is defined, the statistics of the dynamic displacement, strain and stress responses of the structural system can be determined. The intensity of the dynamic responses about their mean values is usually expressed in terms of the standard deviations, which are identical to the root-mean-square (RMS) values when the means are zero. A variety of criteria can be used for the safe design of components depending on the failure modes. In many cases, response values exceeding some deterministic strength value can jeopardize the design and the probability distribution of the extreme value of the response plays an important role in the safe design of structural system. A typical criterion for this type of situation is the three-sigma criterion. Consider S (t) to be a zero-mean sta- tionary response process with RMS 0's and the deterministic strength to be R. A simple criterion for design criterion is R2365 (1.1) ; 3 5' fs“) Figure 1.2 Gaussian response process with three-RMS as design If S is a Gaussian process, it will exceed 305 only for about 0.3% of the time. Fig. 1.2 shows a typical response process and the probability density function of the response f 3(5). A more accurate approach is to assure a desired reliability, i.e. P[IS(t)l < R, O < t < T] = P R. It can be shown that this probability is closely related to Us. Thus Us is a primary response quantity to be computed in random vibration analysis. The step involved in stationary random vibration analysis is illustrated Fig. 1.3. EXCITATION O Random O Stationary Process STRUCTURE O Known and Deterministic O Linear or Equivalent Linear System RESPONSE O Stationary Process O RMS Value l DESIGN CRITERION O Three-Sigma Design O Reliability Figure 1.3 Overview of random vibration analysis CHAPTER 2 Literature Review 2.1 Material Nonlinearity of FRPs Most FRPs generally exhibit physically nonlinear stress-strain behavior. Experi— ments show that stress-strain response in the longitudinal direction (fiber direction) of a unidirectional lamina is quite linear. However, the stress-strain response in the transverse to the fiber direction is often mildly nonlinear and the shear stress-strain response is quite nonlinear. This nonlinearity is mostly caused by the nonlinear matrix (resin) materials. This means that for laminate orientations and loadings which induce shear deformation, the resulting response will be significantly nonlinear. Examples of FRPs with nonlinear stress-strain behavior are: boron/epoxy and graphite/epoxy which have highly nonlinear shear behavior; boron/aluminum, which displays nonlinear behavior for extension trans- verse to the fibers as well as shear nonlinearity; and carbon/carbon which displays nonlin- earities in all principal material directions. Therefore, the analysis of FRPs should consider material nonlinear effects that should be based on properly formulated and suffi- ciently accurate constitutive models. Consequently, in recent years considerable research has focused on the development of nonlinear constitutive models for FRPs. The many models that have been proposed can usually be classified into two groups: nonlinear elas- ticity theories and classical incremental plasticity theories. Petit and Waddoups (1969) simulated the nonlinear stress-strain response of a lam- ina by using a piecewise linear method. The loading is proportional and is applied in small increments to compute the current stress-strain states. In this method, the nonlinear lamina material behavior is accounted for by incrementally applying the average laminate stresses and using a simple integration procedure. The first increment in the laminate strains is cal- culated with the assumption that the laminate behaves linearly over the applied stress increment. This statement becomes, [As] = [A'lnlel (2.1) n+1 n+1 where [A'] is the initial laminate compliance matrix. The increment in the laminate strains, A8 , is added to any previous strains to determine the current total laminate strain [51,, H = [8],. +[A81n +1 (22) Adams (1970) and Foye (1973) used finite element techniques to study the nonlin- ear stress-strain behavior of a composite lamina by considering the fibers and matrix prop— erties at the micromechanical level. Adams proposed the finite-element micro-mechanical model based on a unit-cell structure. Two particular periodic arrays of fibers were exam- ined: a rectangular array which includes the square array as a special case, and diamond array which includes the hexagonal array as a special case. The transverse response of uni- directional composite materials were investigated by using a plane strain finite element program. Hahn and Tsai (1973) presented a constitutive equation for in-plane shear stress and strain of a composite lamina by studying boron/epoxy and graphite/epoxy FRPs. Under axial loadings in the longitudinal and transverse directions of unidirectional FRPs, they observed a linear stress-strain relations, while under shear loading, they observed a nonlinear strain-stress relation. They modeled the nonlinear behavior by using a comple- mentary elastic energy density which contained a fourth order term for in-plane shear stress. This model predicted the strain response under uniaxial offaxis loading (i.e., load- ing in other than principal material directions) fairly well compared to experimental mea- surements. The strain-stress relation they developed is where 66 , 426 are shear stress and shear strain and S 66 , S 6666 are material constants. In this approach no distinction is made between tensile and compressive behavior of the lamina, and problems of failure are not considered. Hashin, Bagchi, and Rosen (1974) presented constitutive equations by curve fitting of the experimental stress-strain data for a composite lamina using the Ramberg-Osgood model which is also capable of representing the nonlinear behavior. Sun, Feng, and Koh (1974) modeled composites by employing nonlinear matrix layers alternating with effective linearly elastic fibrous layers. Since the nonlinear behav- ior of composites are essentially due to the matrix materials, they separated the effect of the matrix from that of the fibers. The matrix material may be taken to be isotropic. They replaced a layer of unit thickness of the unidirectionally reinforced composite medium by a pair of matrix layer and an effective reinforcing layer. The thickness of the matrix layer and the effective reinforcing layer are denoted by dm and d f , respectively. The values of the thickness are assumed to satisfy d! df+dm = n (2.4) where n is the volume fraction of the fibers. Since in most FRPs only the matrix material behaves nonlinearly, they assume that the effective reinforcing layer is orthotr0pic and lin- early elastic with the stress—strain relations given by ' f ofx Cf 11 Cf 12 0 5 x 01y = Cftz Cf22 0 8f? (25) Ony _ 0 0 €1.66 8fry It should be noted that the elastic constants Cfij are not the same as the fiber elastic con- stants. The constitutive relations for the nonlinear matrix are expressed in the form 10 ‘5 x Cm11+Cm11 Cm12+Cm12 0 3 I “my = Cm12+5m 12 szz-t-szz O 5my (2'6) Gmxy _ O O Cm66+5meg mey where Cmij is a nonlinear stress-strain relationship proposed by Kauderer (1958). The gross stress-strain relations for the composite can be derived by the procedure used by Post- ma (1955). Jones and Nelson (1975) modeled nonlinear material properties as functions of the strain energy density, U. If the material properties change monotonically for increasing values of the strain energy, a unique relationship exists for the properties in terms of U. The curves between material properties and strain energy can be obtained from uniaxial stress-strain data for the material and approximated by C.- Material Propertyi = Ai[1 —Bi(-Ug-) ] (2.7) 01' where i is the designation of the specific material property. The energy function, U is given by The material constants A ,- , B i , and C i are the initial slope, initial curvature, and change of curvature of the ith stress-strain curve. The quantity U 0 is used to non-dimensionalize the strain energy portion of the material property equation. The J ones-Nelson model is used in an iteration procedure which converges to the state of stress and strain with a secant mod- ulus corresponding to an equivalent linear elastic body. The model accurately predicts the behavior of ATJ-S graphite under uniaxial loading. Jones and Nelson (1976a) also predict- ed biaxial softening behavior in tension and compression of ATJ-S graphite. They present- ed further theoretical-experimental correlation of a nonlinear material model for ATJ-S graphite (Jones and Nelson 1976b). 11 Sandhu (1976) represented experimental stress-strain data by using cubic spline interpolation functions. The splines yield a smooth stress-strain curve from which accurate moduli of elasticity over the entire range of the curves can be determined. Using the stress-strain data, the behavior of laminates under incremental loading can be predicted. To define the general incremental constitutive relationship for anisotropic materials, it is assumed that l) the increment of strain depends upon the strain state and the increment of stress and 2) the increment of strain is proportional to the increment of stress. Using these assumptions, the incremental constitutive law can be written as where dei and do]. are strain and stress increments and the incremental compliances Si}- are functions of strain 81-. The Sandhu (1976) and Petit and Waddoups (1969) techniques are similar in that both methods allow all of the lamina stress-strain curves to be nonlinear. However, Petit and Waddoups used a piecewise linear function, while Sandhu used a piecewise cubic spline function. In the Petit and Waddoups' technique it is assumed that the tangent modu- lus relative to one material axis is not affected by the presence of stresses in the other directions, whereas the effect of transverse stresses is allowed for in the determination of the tangent moduli in Sandhu's technique. Lamination theory is used to compute laminate strain increments under the applied stress increment. These laminate strain increments are added to those obtained previously to determine current strains in both methods. But in the Petit and Waddoups approach these strains are used to compute the laminate compliance for the next load increment, while in Sandhu's analysis the strains are used to determine average laminate compliances for the same load increment. Further, Sandhu used the revised compliances to compute a new set of laminate strains and iterated this procedure until the difference between two consecutive sets of laminate strains was less than a pre- scribed tolerance. 12 Jones and Morgan (1977) extended the Jones and Nelson (1976) model which was only applicable up to a specific strain energy value. They removed the limitation by using extrapolations of the available stress-strain curve and mechanical property-energy curve for strain energies above available stress-strain data. An interesting feature of their model is the ability to simulate the nonlinear behavior of material in all principal material direc- tions. Amijima and Adachi (1979) presented a simplified method of predicting the non- linear stress-strain curves up to the tensile and compressive failures for an unidirectionally orthotr0pic lamina, symmetric biaxial and triaxial laminates. The analytical procedure is based on classical laminated plate theory. By applying this theory to the small stress or strain increments of the stress-strain curve, the complete nonlinear stress-strain curve is predicted for various laminations. Based on experiments, they assumed the longitudinal and transverse stress-strain responses to be linear and the inplane shear stress—strain curves in the principal axes to be strongly nonlinear. The nonlinear inplane shear stress-strain curve was divided into small sections, and each section was considered to be linear. Dvorak and Bahei-El-Din (1979) used a micromechanical model consisting of elastic fiber and an elastic-plastic matrix. They explored the elastic-plastic behavior of FRPs with self-consistent models. The overall plastic strains and stresses were calculated for selected material systems, using Hill's (1965) model and a modified scheme. Kuppusamy, N anda and Reddy (1984) presented a materially nonlinear analysis of laminated FRP plates. The modified Ramberg-Osgood relationship was used to compute the principal elastic moduli for each lamina. The modified Ramberg-Osgood law is (Ei—Ep)e m l/m O1P o = (2.10) P and based on this the tangent modulus defined as 13 E-—E E, = d—G = ‘ P +13 (2.11) de m (l+m)/m p [Mari—El,» ] 0'P The tangent modulus is computed at every strain level and the stiffness matrix in the finite element formulation is updated accordingly. Sun and Chen (1989) developed a one-parameter fiow rule for orthotr0pic plastic- ity to describe nonlinear behavior of FRPs. Since most FRPs exhibit insignificant plastic- ity in the longitudinal direction, a plastic potential containing a single parameter was found adequate. Effective stress and plastic strain increment were introduced to yield a universal nonlinear stress-strain relationship for the orthotr0pic composite. They further developed the theory presented by Kenaga et a1. (1987) by incorporating the fact that for most unidirectional FRPs the stress-strain relation in the fiber direction is essentially lin- early elastic. The plastic behavior of composites is then characterized by a single parame- ter. Nanda and Kuppusamy (1991) developed an elastic-plastic model for anisotropic materials including both perfect and work-hardening plasticity. The incremental formula- tion with isotropic work-hardening was used. In the incremental formulation it is assumed de = de‘ + dep (2.12) where d8 is the incremental strain, dee is the elastic incremental strain, and de” is the plastic incremental strain. The elastic incremental strain is related to the incremental stress through Hooke's law for anisotropic materials: de" = D”'do (2.13) where D is the elastic constitutive matrix. The plastic incremental strain is obtained from the normality condition: l4 8F(o, k) I). de —}t do (2.14) where F (0', k) is the yield function, k is the work-hardening parameter, and it is a propor- tionality constant to be determined from the consistency condition. Vaziri (1991) developed a plasticity—based constitutive model for the nonlinear behavior of laminated FRPs up to and including ultimate failure. The proposed model for single layers of FRPs was derived within the framework of rate-independent theory of orthotr0pic plasticity. The elastic constitutive equation for orthotr0pic material in plane stress is given by do, = @3de = ij(dej—def) (2.15) where doi and dei are total stress and strain increments respectively. Superscripts e and p denote elastic and plastic portions. During plastic loading, if the consistency condition df = 0 ,where df is the quadratic criterion for yield of an orthotr0pic material is satisfied, the final constitutive equation leads to ‘3 a a e. do.= 9— "" "’ " "J de. (i,j=l,2,6) (2.16) ' i " uH'+a..Q;..a.i ’ The individual layer constitutive equations are superimposed using classical lamination theory to yield the global laminate response. 2.2 Laminated Plate Theories Laminated plate theories, originating from the classical theory based on Kitch- hoff’s assumption, have been significantly refined over the last two decades for application to laminated FRP plates. All theories begin by means of a hypothesis about how the dis- placements vary in the thickness direction of the laminates. The assumed displacements are usually expanded as polynomials with respect to the thickness coordinate. Laminated plate theories can be divided into three categories depending on the type of assumed dis- 15 placement field: global theories (also called single-layer laminate theories), layerwise the- ories and zig-zag theories. In early studies, classical plate theory was used to analyze the laminated FRP plates. It is well known fact that the classical laminate plate theory (CLPT) based on the Kirchhoff’s assumption stipulating that normals to the midplane before deformation remain straight and normal to the plane after deformation underpredicts deflections and overpredicts natural frequencies and buckling loads. These results are due to the neglect of transverse shear strains in the classical theory. The displacement equations of CLPT are u = uo(x.y)—z(%¥) 2.17 v(x, y, Z) 8w) ( ) v0(x, y) _ Z('a_y W(x. y. z) = wo(x. y) where u and v are the in-plane displacements with respect to x and y coordinates, respec- tively, and w is the out-of-plane displacement. dw/ax and dw/ay are flexural angles about the y and x axes. Because transverse shear deformation is neglected, the CLPT can yield significant error for thick plates. For FRP laminates, the ratio of the effective elastic in-plane modulus to the transverse shear modulus is much higher than that for isotropic ma- terials. These high ratios render the CLPT inadequate for the analysis of most FRP plates. An adequate theory must account for transverse shear deformations. A more refined theory is the first-order laminate plate theory. The first-order shear deformation theory (FSDT) for isotropic and homogeneous plates was developed by Reissner (1945) and Mindlin (1951). They assumed the same level of approximation resulting in displacements of the form “(x9 y: Z) = “0(x9 y) + ZWx(x9 y) We y. z) = vo(x. y) + zwy(x. y) (218) W(x, y: Z) = W0(x’ y) 16 where uO , v0 , w0 denote the displacement components of a point (x, y) on the midplane, and nix , my are the rotations of the normals to the midplane about the y and x axes, respec- tively. The laminated version of the FSDT was developed by Yang, Norris, and Stavsky (1966) with the same assumed displacement fields. The FSDT is based on a linear distribu- tion of the in-plane displacements and a constant of the transverse displacement through the thickness, which results in nonzero transverse shear stresses at the plate's bounding planes and hence necessitates a shear correction factor. Yang, Norris, and Stavsky (1966) only considered the frequency equations for the propagation of harmonic waves in a two-layer plate. Whitney and Pagano (1970) improved this theory by introducing reduced stiffnesses and a shear correction factor. The FSDT yields a constant value of transverse shear strain over the thickness, although it is known that in reality the distribution of trans- verse shear strains is nonlinear. Higher-order shear deformation theories that include high- er-order terms in the displacement fields are capable of capturing this nonlinearity, and the third-order theory is popular. Moreover, when a parabolic distribution of transverse shear stresses across the plate thickness is adopted, the boundary condition on the bounding planes of the plate could be satisfied and there is no need for a shear correction factor. Lo, Christensen and Wu (1977) presented a third-order displacement field which includes both in-plane and out-of-plane modes of deformation. Their theory gives a much better approximation to the exact elasticity solution for laminated plates than CLPT. They assumed the following displacement field is assumed which is consistent in the sense that the transverse shear strains due to in-plane displacements u and v are of the same order in 2 as those determined by the transverse displacement w: u(x. y. z) = u0(x. y) + zv.(x. y) + 2290:. y) + z3¢.(x. y) v(x. y. z) = vo(x. y) + zv,(x. y) + zzc,(x. y) + 23¢,(x. y) (2.19) m. y, z) woo, y) + aux. y) + zzczcc. y) 17 where “0’ v0, W0 and tux, w), are physical DOF as for the FSDT, while 3; and (i) may be viewed as higher-order rotations. Reddy (1984) developed the third-order shear deformation theory (TSDT) of lami- nated FRP plates with the same dependent unknowns as the FSDT. Reddy’s (1984) theory is not only capable of including transverse shear strains, but also of representing the para- bolic distribution of the transverse shear strains through the thickness. Consequently, there is no need to use shear correction coefficients when computing the shear stresses. It yields deflections and stresses that are more accurate than the FSDT and satisfies the zero trac- tion boundary conditions on the surfaces of the plate. This theory is adopted for the non- linear random vibration analysis developed in this thesis. Details are given in Chapter 3. Reddy’s theory gives good results for overall responses such as global deflections of FRP laminates, and is computationally efficient. However, it can give rise to errors in local effects such as interlaminar stresses and strains. Layerwise theories must be used if accurate interlaminar stress are required. Lay- erwise theories consider each FRP layer individually and enforce the condition of continu- ous interlaminar stresses at all layer interfaces. These theories yield excellent results for global and local responses, but because each layer is considered separately, the number of DOFs increase rapidly when the number of layers is large and result in high computational costs. In 1987, Reddy proposed a layerwise displacement plate theory. In his layerwise theory, the displacement field is expanded as a linear combination of the thickness coordi- nate and undetermined functions of position within each layer. “(1,y,Z) = “0(X, y) + U(x: y! Z) We y. z) = vo(x. y) + V(x. y. z) (2.20) W(x, y: Z) = W0(x, y) + W(x9 Y: Z) where, 18 Q 11 2 U,(x.y)<1>,-(z) j 1 Vj(x, y)J-(z) (2.21) l < H 11' M a 11 J W ,(x, y1 ,(z) 1 E II “M” J and u0 , vO , and wO are displacement components of the midplane, and n is the number of layers of the laminate. The (1)}- denote the interpolation functions of the thickness coordi- nate and U j , V j. and W j are nodal values of U, V and W at the nodes in the jth layer. The functions (1)]. are piecewise continuous functions, defined only on two adjacent layers, and can be viewed as the interpolation functions associated with the jth interface of the layers through the laminate thickness. Because of this local nature of (Dj, all transverse strains will be discontinuous at layer interfaces, allowing the possibility that the interlaminar trans- verse stresses computed from constitutive equations can be continuous. Lu and Liu (1992) presented a layerwise laminate theory which accounts for both the interlaminar shear stress continuity and transverse shear deformation. The theory gives excellent results for both stresses and displacements and more importantly, the interlami- nar shear stresses are computed directly from the constitutive equations instead of being recovered from the equilibrium equations. They developed an interlaminar shear stress continuity theory with a refined displacement field considering Hermite cubic shape func- tions through the thickness. The nodal displacements and rotations are used as indepen- dent variables. Another class of accurate theories are layerwise zig-zag theories. These theories use displacement fields similar to those discussed for single-layers but apply the condition of continuous transverse shear stresses at layer interfaces to the in-plane displacement fields. They assume piecewise polynomial variations of the in-plane displacements through the thickness and can satisfy the conditions of continuous transverse shear l9 stresses. DiSciuva (1985) proposed a first-order zig-zag theory where a piecewise linear variation through the thickness was imposed on the in-plane displacements. His assumed displacement fields are 8w “(xt y: Z) = “0(x’ y) + Z(Yx - —ax_0) + 2¢k(z — Zk) k a 2.22 v(x, y, z) = v0(x, y) + {‘5 - 7:9) + 29,42 — zk) ( ) k W(x. y. Z) = wo(x. y) where u0 , v0, and wO are the values of the displacements on the reference surface and 7x and yy are shear rotations. The (pk and 0k are functions to be determined by satisfying the continuity of interlaminar shear stresses at interface k. Unlike layerwise theories, this the- ory has a finite number of DOFs regardless of the number of layers in the laminate, and can therefore result in computational savings. However it does not satisfy the traction free boundary conditions at the surfaces and yields constant transverse shear stresses through the thickness. Many higher-order zig-zag theories such those by Murakarni (1986) and Di- Sciuva (1987) have been developed to overcome those drawbacks. More details can be ob- tained from the excellent review papers by Noor and Burton (1989), Kapania and Raciti (1989) and Reddy (1990). 2.3 Equivalent Linearization Techniques Random vibration analysis of structural systems has been deve10ped considerably over last three decades. However, in contrast to deterministic dynamic analysis, random vibration analysis is more complicated and for nonlinear problems exact solutions cannot often be obtained. Therefore various approximate methods have been introduced. Among them, the equivalent linearization technique is a versatile method capable of handling even strong nonlinearities. In this method, equivalent linear equations of motion are obtained 20 from the nonlinear equations by minimizing the difference between the two equations in a least square sense. Equivalent linearization was first proposed by Caughey (1959), who considered the replacement of a nonlinear oscillator by a linear one with properties computed through a mean-square minimization criterion. Foster (1968) first applied the equivalent lineariza- tion technique to multi-DOF systems with random excitation. He derived the equation of motion which defines optimum mathematical stiffness and damping matrices for a linear- ized system. However, his approach is not applicable to general problems. Later, Iwan and Yang (1972) formulated an approach to solve general multi-DOF systems with random excitation. Their approach is based on the concept of decomposing the nonlinear function into a sum of simple approximate nonlinear terms, each of which depends on the displace- ment and velocity vectors. However, the choice of approximate terms is an important fac- ' tor controlling the accuracy of the solution. Also, when mass is introduced, the identification of these terms can become complicated. The most general method for equivalent linearization was introduced by Atalik and Utku (1976) for structural dynamic problems. This approach had previously been applied to electrical engineering problems in Russia by Kazakov (1965). Atalik and Utku consid- ered a general multi-DOF nonlinear system in which the nonlinearities are represented by a vector, {4)}. They showed that the coefficients m c,” and kij which minimize the i1" 1 magnitude of the mean square error {e} = {¢}-[M]{f}-[Cl{i}-[Kl{x} (2-23) have the forms 21 29¢,- mi]- : Elf-xi] -34” C.) = ELIE? (2.24) -34); k” = Eta—x}- Much attention also has been given to modeling nonlinear hysteretic systems. Wen (1976) developed a generalized restoring force model and derived closed form solutions for the linearization coefficients by assuming a Gaussian response. Due to its versatility and stability in obtaining convergence, this model has been widely used in random vibra- tion analysis. It was further refined by Barber and Noori (1984) to model damage accumu- lation, degradation and additional hysteretic characteristics. Many modified versions of equivalent linearization also have been introduced. Izumi, Zaiming and Kimura (1989) presented a version of linearization where the expecta- tion is taken with suitable weighting functions. Iyengar (1988) introduced a method which approximates a nonlinear system by a higher-order linear system. Equivalent linearization was applied to nonlinear systems with non-white noise excitation, when the excitation is treated as the output of a linear filter excited by white noise (Socha 1990). Also, an equiv- alent linearization technique for systems subjected to nonstationary random excitations was introduced by Iwan and Mason (1980). Kimura and Sakata (1981) considered non- symmetric nonlinear systems subjected to nonstationary random excitations. Many gener- alizations and applications are summarized by Roberts and Spanos (1990) and Socha and Soong (1991). 2.4 Random Vibration Analysis of FRP Laminates The majority of investigations on FRP laminates were accomplished within the framework of static and deterministic dynamic theories. FRP structures in plate and shell forms are used for many engineering applications because of their high performance char- 22 acteristics. Many of these FRP structural components are exposed to randomly varying dynamic loading, e.g., aircraft components, marine structures, etc. However, relatively few investigators have explored random vibration analysis of FRP laminates. Witt and Sobczyk (1980) analyzed the dynamic response of a symmetric cross-ply laminate in cylindrical bending to an exponentially decaying random load. The effect of the fiber orientation and correlation parameter of the external load on the mean-square ver- tical displacement was investigated. Witt (1986) later extended the formulation using a first-order shear deformation theory to a square antisymmetric angle-ply laminate. An approximate solution was obtained by applying the Galerkin method and displacements were developed using a series of independent orthogonal linear functions satisfying the boundary conditions of the plate. Cederbaum, Librescu and Elishakoff (1989) used the first-order shear deformation theory to analytically study the linear random response of FRP laminated plates. They determined the mean square response of a cross-ply laminate subjected to a random noise with exponential decaying correlation function, and of an angle-ply laminate subjected to white noise. They extended the formulation to include the high-order shear deformation theory developed by Librescu (1986). They dealt with two random load fields: In the first case, the random excitation was modeled as a point load, random in time with constant spectral density; in the second case, it was modeled as a turbulent boundary layer pressure fluctuation (1988). Mei and Wentz (1982) presented an analytical methods for determining the large amplitude response of antisymmetric angle-ply laminated rectangular plates subjected to wide band random excitation. Karmen-type geometric nonlinearity was considered in the formulation. They considered only a single mode Galerkin approach to obtain a nonlinear modal equation and employed the equivalent linearization method on this equation. How- ever, transverse shear effects were neglected and only in-plane responses were considered. Later Mei and Prasad (1989) derived nonlinear equations of motion for symmetric lami- 23 nated plates considering transverse shear deformation using a first-order shear deforma- tion theory. It is apparent that the first-order shear deformation theory is widely used in random vibration analysis. Singh et al. (1989) studied the random response of an antisym- metric angel-ply FRP plate subjected to random lateral load applied on its surface. Harichandran and Hawwari (1992) performed nonlinear random vibration analysis of a single layer FRP plate loaded in extension. This was the first application of the finite element method and equivalent linearization for the nonlinear random vibration analysis of FRP plates. They considered nonlinearity in the in-plane shear stress-strain law. Later, Harichandran and Naja (1997) extended the formulation to the bending deformation of FRP laminated plates. They showed that the nonlinearity in the response is significant for large excitations. However, their formulation was limited to the classical laminate theory and transverse shear deformation was not considered. CHAPTER 3 Finite Element Formulation 3.1 Introduction As the use of FRPs in structural elements has increased, many engineers have devoted considerable attention to the behavior of FRP structures. Typical structural ele- ments made of FRPs include laminated plates and shells. These structures are composed of plies or layers of fiber reinforced material that are bonded together. The FRP structures studied in this work are FRP laminated plates. Numerous FRP laminated plate theories with varying degrees of complexity and accuracy have been developed. However, analyti- cal solutions of these theories are only available for a very limited set of geometries, load- ings and boundary conditions. In general, it is necessary to use approximate methods which can be applied to complicated structures with arbitrary shapes, loads and support conditions. The finite element method, in which the plate is discretized into numerous ele- ments, is one of the most versatile approximate analysis techniques. Within various types of finite element formulations, the simplest and the most widely used one is the displace- ment-based finite element method. In this method, the primary unknowns are taken to be the displacements at the nodes of the discretized plate. Displacements within a discrete element are interpolated from the displacements at the element nodes by means of interpo- lation (or shape) functions. Applied loads are replaced by equivalent nodal loads. Mass, damping and stiffness matrices are computed for each element and then assembled to form global matrices which describe the behavior of the entire plate. After specifying the restraints, linear algebraic equations or a set of ordinary differential equations are solved for static or dynamic problems, respectively, to determine the nodal displacements of the plate. Derived quantities, such as strains and stresses are computed from the nodal dis- placements. 24 25 3.2 Nonlinear Constitutive Equations Based on experimental results, it is well known that a unidirectional composite lamina displays essentially linear behavior when the direction of loading is parallel to or perpendicular to the fiber direction. However, when loaded in shear, the behavior is signif- icantly nonlinear behavior. Fig. 3.1 shows the global coordinate system x-y and the local or material coordinate system 1-2 for a typical unidirectional composite element. Several nonlinear constitutive models have been proposed, of which the Hahn and Tsai (1973) model is simple and easily implemented into a finite element code. Hahn and Tsai used the complementary energy density theory to derive a strain-stress relation. Their strain energy density function, which includes a fourth order term, is expressed as — 1 2 l 2 1 1 2 l 4 where W“ is the reduced complementary energy density function. The strains are related BW“ a—oi'r are stresses and strains respectively. Therefore the strain—stress relation is given by to the complementary energy density through, 8,- = i = l, 2, 6 , where o and e k\\\\\ ////,. W . Figure 3.1 Global and material coordinate systems for a lamina ..\\\\ X 26 81 ‘511512 0‘ 01 0 2 82 = 512522 0 62 +S*6616 0 (32) in which 81 and 01 are the normal strain and stress in the fiber direction, 82 and 0'2 are the normal strain and stress normal to the fiber direction, and 76 and 16 are shear strain and stress in the in-plane material coordinates. The 3x3 square matrix is the linear compli- ance matrix [S]. The constant, S * 66 describes the nonlinear behavior of the composite lam- ina loaded in shear. The Hahn and Tsai model is limited to two-dimensional cases with only in-plane stresses and strains. In this work, the model is extended to three-dimensional cases by assuming ”e ‘ ' ' ro ‘ O 1 811512 0 0 O 1 O 32 512522 0 0 O 02 3 * (76>: 0 05660 o<16>+156616> (3.3) 3 74 o o 0 S44 0 14 55414 where subscripts 4, 5, and 6 are equivalent to shear along the material directions 2-3, 1-3, and 1-2, respectively. The shear stress-strain relations may be written as tq = gq(yq), q = 4,5,6, where gq(yq) is the real root for 'tq of the function 5* “1:: + S quq - 'yq = 0. Separating out the linear term, g q('yq) may be expressed as ngq) = (fqoqnyl-qqhq. q = 4,5,6 (3.4) Consequently, the complete stress-strain relation, which is the inverse of the strain—stress relation may be written as 27 fo“ —Q11 Q12 0 0 0 f8“ ”00 O O 0 - r8“ 52 Q12 Q22 0 O O 82 00 0 O O 82 <16): 0 0 Q66 0 0 1 1 zV Figure 3.5 Assumption for higher-order shear theory tions. The requirement that the transverse shear stresses vanish at the top and bottom sur- faces of the plate, i.e. h Tyz(x, y, :2) = 0 h (3.19) sz(x, y, ii) = 0 where h is the thickness of the plate, is equivalent to the following conditions on the corre- sponding strains: (3.20) 35 Using the strain-displacement relations 3v 8w 2 3w Y)'z=a—Z+§;=Wy+22§),+3z gy'i'a—y: 32 ( - 1) a a 2 a M = 8—:4-3—‘2? = wx+22§x+3z §x+§—:-’ and settin - 0 at - +13 ields g sz — Z — *2 ’ y - +211§+3é2§ +a—i-v-O ~sz‘wy 2 y 2 y By- h 2 a (3.22) h w ,2 . ,y.2(-§)g,.3(-§) 9.5; = 0 Solving the above two equations for Q and C), yields t, = o 4 8w (3.23) Sy = “372(33- + lily) Similarly, using the condition 'sz = 0 at z = i; , yields 6, = 0 4 3w (3.24) c. - 37,305: +le Substituting Eqs. 3.23 and 3.24 into Eqs. 3.18 and simplifying them, the displacement field in Eqs. 3.18 becomes _ '- 4 z 2 3W 1 u — “0(x: y)+z Wx_§(}-1)(Wx+0;)- (3.25) 2 _ V = ”0(X’Y)+ZLWY-3(h) (“Yuk-331 W0(x7 y) E 11 Eq. 3.25 indicates that in-plane displacements have a cubic variation through the thickness direction. Therefore, the transverse shear strains have a parabolic distribution through the 36 thickness of the plate and satisfy the stress free condition at the top and bottom surface of the plate. Consequently shear correction factors are not needed because the parabolic shear stress/strain distributions are properly represented. 3.5.2 Strain-Displacement Relations The linear strain-displacement relationships for small deformations are: _Q_u___a_“_0+3‘fi_ 3_4_ hZwa+al20 8x 8x Z3x Z_ 3x 8x2 0v __ 3120+ aWy -3—Z3 4 _:2(a‘|,y+ a _sz] 8’ = 3)” _ 3): Za—y 3y 3)? _Bw 82 :5; = O _au+ 8v 'ny —-a-y +0; (3.26) +— By 8x 2 3y +$+ axay auo 9% “(al. LWy)_ 3 4 Jay. av, 32%] = + z — -- 2 d—y B—x 8v aw _ aW0 242y+( +3_w0) sz=az+ a—y W)’ y+-87—Zh2 3y _8_u +11_w_ 8W0 24 +8w0 xz‘ a? ax “5+? Zfzi '07) Note that the above equations contain higher-order derivatives of WO than the first-order equations. This is an important feature when the displacement-based finite element formu~ lation is performed. The strain-displacement relation can be written in the following con- cise matrix form 37 EyO Kyl 2 0 3 xy3 {8} =1Yx20i+zixxy1i+u 0 Hz < “mi (3.27) szO 0 Kyzz 0 L‘YXZOJ L O 2 kazzJ \ O J where auo _ 301x _ 4 801x dzwo ‘0 - 0)? x1 - 8x ’ sz' - 3h2 8x 8x2 2 avg 30!), 4 3w), 8 W0 8‘0 _’ '1 = _’ K 3 + ) a); 2 a); 2’ 3’12 6); ay2 \l’ +8w0 4 \1’ +8w0 73'20 = y —’ 1'22 = ”—2( y —) 3y h 3y (3.28) dwo 4 dwo 7x20 = “Ix-+7, szZ " —;l—2(Wx+—87) duo avo 301x aw), 190 - 8—y+0—x_’ m - By ‘33: 2 1c 4 801x + 01.1, 3 W0 ”’3 3,12 3y 8x dxdy 3.5.3 Stress-Strain Relation The stresses in a typical kth lamina in global coordinates can be written as ro "‘ 2' - ' qk'e i x Q11 Q12 Q16 0 0 x or Q12 Q22 Q26 0 0 8? ”'9’ = Q16 Q26 Q66 0 0 ‘nyi (329) TYZ 0 0 0 Q44 Q45 YR thz. _ O O O Q45 §55_ thzl 38 where the Qij (i, j = 1,2,4,5,6) are the transformed material stiffness coefficients for the kth lamina. Since the normal strain ez is neglected in Reddy’s (1984) theory, the corresponding normal stress oz can be omitted. Only linear part of the constitutive equations is considered here. The nonlinear part of the constitutive equations is considered in Chapter 4. Typical normal and transverse shear strain variations through the thickness of sym- metric laminates are illustrated in Figs. 3.6 to 3.8. Reddy’s theory gives the parabolic dis- tribution of transverse shear strains while FSDT gives a constant distribution of transverse shear strains within each layer. However, it should be recognized that a theory with a glo- bal assumed displacement field yields discontinuities across layer interfaces for some shear strains that should be continuous. In Figs. 3.7 and 3.8, exact 3D analysis would yield continuous distributions for 713 and 723 in the z—direction. Thickness Normal Strain Figure 3.6 Typical normal strain variation through the thickness 39 1 Thickness L A —- Reddy - - FSDT Shear Strain 713 for Three-Ply Figure 3.7 Typical transverse shear strain 713 variation through the thickness I 1 . 2 - I l I l l I m i I 8 . l i . ' .9. ' l5 1 . l v I l l I l ' — Reddy ’ : — - FSDT Shear Strain 723 for Three-Ply Figure 3.8 Typical transverse shear strain 723 variation through the thickness 40 3.5.4 Finite Element Formulation Based on the displacement field proposed by Reddy (1984), a finite element is developed using a 4-noded rectangular element with 7 degrees of freedom at each node. The isoparametric formulation is considered to map from a rectangular to arbitrary quadri- lateral shape. A typical 4-node plane isoparametric element is shown in Fig. 3.9. In iso- parametric elements, the natural E_,— 11 coordinate system is used. For a rectangular element with dimensions 2a and 2b with x = O and y = 0 at the element center, the dimen- sionless isoparametric coordinates E and n are g = x/ a and t] = y/ b. The consistent mass matrix can be used to represent a continuous distribution of mass. The same shape functions used to compute the element stiffness matrix are used to compute mass matrix. The consistent mass matrix is more accurate than a lumped mass matrix for flexural prob- lems, involving beams and/or plates (Belytschko 1976). 3.5.4.1 Extensional Deformation Consider the plane rectangular element with the nodal displacements shown in Fig. 3.10. Figure 3.9 Typical 4-noded isoparametric quadrilateral element 41 T {“16} = {11], 121,142, V2, “39 V3, “4, V4} (330) This element is based on the bilinear displacement field: “(KS1 T1) = a1 '1' “2S + “311 '1' “4%" (3.31) V0“, 11) = b] + 1225+ b3" ‘1' [24:11 where “0 and v0 are in-plane displacements, and l; and n vary from -1 to +1. The dis- placements within the element are can be interpolated from the nodal degrees-of- freedom (DOF). The in-plane displacements are interpolated from the nodal DOF through 0 {ule}= [N'l{u1e} (3.32) {“0{9"1}}_ NI 0 N2 0 N3 0 N4 0111,01v201v301v4 V01; 11} — where N i , i = l, 2, 3, 4 , are interpolation functions given by V‘4 ZA 1'] = y/b v3 u4 / / / u3 4 3 + > § = x/a “1 1 2 112 7 / ’ v1 v2 Figure 3.10 Nodal displacements for extensional deformation of a plate element 42 N860) = fid—th—m N2(§.n) =5<1+t)(1-n) (3.33) N3(§.n)=}1(1+§)(1+n) Nacn) = j(l-§)(1+n) 3.5.4.2 Shear Deformation The out-of—plane shear rotations are interpolated from the nodal shear rotations shown in Fig. 3.11 with the same shape functions used for the in-plane displacements: 0 {uze}= [N'l{u2,.} (3-34) 0N10N20N30N4 {v.{§.n}} _ NI 0 N2 0 N3 0 N4 11515.11} where {“212} is a vector of the 8 element nodal shear angles shown in Fig. 3.11, T {“28} = {Wyl’ —\1”xl! Wyzt "1&2, ‘11))39 _‘Vx3, ‘1’).4, —\I,x4} (3.35) - “‘1’ 3 z n - y/b x A / / +1,” 3 p § = x/a 2 “M » z\l’x2 Figure 3.11 Nodal displacements for shear deformation of a plate element 43 and N i, i = 1, 2, 3, 4 , are the interpolation functions in Eq. 3.33. 3.5.4.3 Bending Deformation The simple rectangular plate-bending element originally developed by Melosh, Zienkiewicz and Cheung is adopted (Melosh 1963, Zienkiewicz and Cheung 1964). This element has three nodal displacements at each node, yielding 12 degrees-of-freedom per element as shown in Fig. 3.12. Therefore12 constants cl. are required to represent the dis- placement field for this element. The chosen polynomial approximation of w is 2 2 w = 61 + 0,: + c311 + c4§ + c590 + 660 (3 36) 3 2 2 3 3 3 + c7§ + c8§ n + 095,1] + c1011 + C119 11 + 612571 which is an incomplete fourth degree polynomial with a complete cubic of ten terms and two additional quartic terms. This polynomial does not preserve slope normal to the edges, i.e. aw/dy is not continuous across the edges. Elements violating any of the continuity conditions are called non-conforming elements, but despite this deficiency, the element -W Z 11 - W “’3 3" i / i/n W3’y 3 > g = x/a w2 2 wz’y Figure 3.12 Nodal displacements for flexural deformation of a plate element 44 gives good results. A comprehensive discussion on plate-bending element is presented by Zienkiewicz and Taylor (1991). In terms of the nodal DOF _ aw, aw, aw, aw, dw3 3w3 dw4 dw4 T 337 {“3e} - an—y'F—apwz, 57-73;:W3: ‘3‘),‘1—E—1W4: 737-87 (. ) the displacement is expressed as W0(§’n) = [N11N12 N13 N211V22 N23 N3111/32 N33 N41N42 N43]{“3e} (3.38) = [N"]{u3e} where N“ = é<1+§§))(1+nn,)(2 +§§1+nn.—§2—nz) N6 = -§bTI.-(1+§§,-)(l-nn,-)(1Him-)2. i= 1.2.3.4 (3.39) N13 = {lg-am 1-§§.~)(1 + 011,-)(1 + €622 3.5.4.4 Strain-Nodal DOF Relations Combining in-plane and out-of-plane displacements and all the nodal DOF together, the complete interpolation can be expressed as {u, v, w}T = [N]{ue} (3.40) where {ue} = {{“1e}{“2e}{“3e}} and [N] = [[N'][N'][N"]] . The strains given Eqs. 3.27 are expressed in four parts in terms of power of the z coordinate. It is also convenient to express the strain-nodal DOF relation in four corresponding parts. The first part of the strain field within an element may be expressed as 45 “Eto‘ 8’0 131 101 {01 ”not = [BI]{ue} = 1 3x8 {ue} (3.41) ”Y [0] [32]2x8 113312,,12 5x28 yzO t‘szOi 1n whtch the subscrrpts indicate the srzes of the matr1ces. The sub-matrices of [B ] are g1v- en by N1,x 0 N2,x O N3,x O N4,x O Nl,y Nl,x N2,y N2,x N3,y N3,x N4,y N4,x N o N 0 N o N 0 [32]: 1 2 3 4 (3.43) 0 N10 N20 N30 N4 Nll,y N12,y N13,) N21,y N22,)» N23,)» [B31= Nll,x N12,}: N13,x N21,x N22,x N23,): (344) N31,)! N32,y N33,y N4l,y N42,y N43.y:l N3l,x N32,x N33,): N41,): N42,): N43,x The notation , x and , y in the subscripts represent partial derivatives of shape functions with respect to x and y, i.e. N aNi laN‘ ' 1 2 3 4 LI - _a—x— - Za—g'! l - 9 9 9 (3 45) 3N, laN, ' N - = 1,2, 3,4 m = E - E'a?’ The elements of the [BI] matrix are therefore 46 l I NM = 7150 —n> NM, = in“ —§> 1 1 N23 = 470-11) NM = 730%) l l N3,x=;{('1(1+11) N3,y=-ZZ(1+§) 1 1 NW = —8—15(n—1)(3§2+n2+n—3) NM = -8ib(§—1)(§2+§+3n2-3) NW = {301+ 1)(n-1)2 N12,y = -§(3n+ 1)(n-1)(§— 1) NW = é<3§+ 1)(§-1)(n- 1) N11,, = 81,574+ Ina—1f N21,}: = §1(n-1)(3§2+n2+n-3) NM = $(&+1)<§2—§+3n2-3) N2“ = {7101+ nm- If N22,), = §(n—1)<3n+1)(§+1) (3.46) N23,}r == §<§+1)(3§- 1)(n— 1) NB,y -= fim— 1>(§+ If N31,x = -81—a(n+ 1)(3§2+n2—n—3) Nay = ~§1;(§+1)(§2-§+3n2-3) N32,}: = gm— 1)(n+ If N32,, = § = [B ]{u.} = 3‘8 {u.} (3.47) [01 [01 [01 O k. 0 J 0 Nl,x 0 N2,x 0 N3“ 0 N4,x [B4] = N1,y 0 NM 0 N3” 0 N4,y o (3.48) N1,x N1,y N2,xN2,y N3,x N3,yN4,x N4,y and each element is given explicitly in Eq. 3.46. The third part of the strain-nodal DOF is L: {Emu .} = u [01 [0] [0] [0] { e} (3.49) [3212x8 [3:512le in which [32] and [B3] are given by Eq. 3.43 and Eq. 3.44. The final part of the strain- nodal DOF relations is F 3}:2 ‘ (‘7)Kx3 (23),. 4 y3 3122 —T ny3 0 >= [B‘V1{u.} = [ O [0] [B4l3xg [3513x12]{u } (3.50) [0] [0] [0] where [B4] is given in Eq. 3.48 and [BS] is 48 [B5] = N11,xx N12,xx N13,xx N2],xx N22,.rx N23,xx N11,»? N12.yy N13.” N21,” N22.” N23.” 2Nll,xy 2NlZ,xy 2N13,ch 2NZ],er 2N22,xy 2N23,xy (3.51) N31,xx N32,xx N33,.rx N4l,xx N42,xx N43,xx N31,» N32.» N33.» N41,”, N42.» N43.» 2N3l,xy 2N32,er 2N33,xy 2N41,er 2N42,xy 2”43,ch The sign , xx , , yy and , xy represent second partial derivatives of the shape functions with respect to x and y. The elements of [BS] are: _ _3_ _ _ 3 Nll,xx " ‘4a2§(ll 1) N11,” - -47211(§- 1) = 1 ”‘2'“ 0 N12,” = -4—b(§— ”(311-1) 1 N21xx=—3§§(n—1) N21 yy= 23n(§+1) ’ 4a ’ 4b N22,xx = O l N22,yy = n<§+l)(3n-l) 1 N23.“ - Lam—”(3“ ” N23.” = 0 3 3 a N = O 1 32m N32,” = 375:4- 1)(3n+ 1) — _i N33,” ’ 40(Tl+1)(3§+1) N33,”, = 0 N41 xx=—3—2'E.9(n+1) N41 yy=_3_2n(§—1) ' 4a ’ 4b N 42’ xx = 0 (3.52) 1 N43,... = -4—a(n+ 1>(3§- 1) N43 W = o 49 N11,".y = —§;—b-(3E,2+3n2-4) N31,”. = —8;——5(3§2+3n2—4) N12,... = ~313m— 1)<3n+ 1) N32,.y = 511+ 1)(3n— 1) N”... = {315%- 1)<3§+ 1) N33,.y = --8—‘,;(§+ 1)<3§- 1) N21,” = g(171—)(3ngr3112—4) N41,”, = 53%;(3533112—4) 1 1 N22,... = 3201— 1)(3n+ 1) N42,.) = —§;(n+1)(3n— 1) l 1 N23,... = 8—l3(§+ 1)(3§— 1) N43,... = —37)(§-1)(3§+ 1) Combining the four parts, the complete strain-nodal DOF relation becomes 4 4 {a} = [B‘1{u.} + ztB"1{u.} + z2(—l-;5)13"‘1{u.}+ z3(—;;5)[Bwl{u,} [[31] + 21811“ zz(_}_:1_2)[Bm] + z3("3ihi)[3w]]{“e} (3.53) = [B]{ue} 3.5.4.5 Linear Element Stiffness Matrix FRP laminated plates consist of several laminae which are assumed to be perfectly bonded together. Each layer may be oriented at an arbitrary angle with respect to the glo- bal coordinates. Therefore the stiffness matrix of each layer must be computed indepen- dently and added to obtain the total stiffness matrix due to all laminae. The linear case is considered here while the nonlinear part is considered in the next chapter. The strain energy is {€}T[Q]{8}dV (3.54) NI!— Ue=j V: where {6} includes all strains including shear, [Q] is the constitutive matrix defined in Eq. 3.12, and V e is the volume of the element. Substituting Eq. 3.53 into Eq. 3.54 yields 50 U. = ;J{u.}TIB1T1Q1181{u.}dV= ${u.}T1K.1{u.} (3.55) where [K e] is the stiffness matrix of the element. For a laminate with N layers, the element stiffness matrix is N z. 1K.1= 2 j j [BthéltBldxdydz k=lA 112k _ N 1T HT 2 4 [HT 3 4 IV (3.56) _ 2 [Hus] +z[B] +2 (75)”? 1 +z(—;1—2)[B 17] k=1—l-lz,,_. — 1 n 2 4 1n 3 4 IV [Q] [81+le 1+2 (7)18 ]+z (372)“? ]]|J|d§dndz Zk—l where h is the thickness and zk is the height from the reference plane to the bottom of the kth lamina. IJ I = (jg-mu, y) is the determinant of the jacobian matrix for transforma- tion from the x-y coordinate system to the £41 system. A A 1 A 2 3 z0 z1 22 h .-- ...... 1 ....... @1396??? .......... Z ZR 1 zk zN 1‘ V V N V Figure 3.13 Geometry of laminated plate and layer numbering system 51 The integral with respect to z in Eq. 3.56 can be performed analytically and other integra- tion in terms of § and n is normally computed by numerical Gauss integration. Expanding Eq. 3.56 can be written as [Kc] = 2 [[klel+[k2el+m +[k15e]+[k16.,]] k=l where 1 l 1k,.1 = (z,—z,._,) J J [BllrléllBllllldédn 2 2 (zk—z [kze] = [k3e] = 3 -l—l 1 l -——2—"—i) J J 1B'1TIQ11B“JIJId§dn —1-1 (3_ 3 ) 1 1 M05.) J J [B‘]T[Q][B"I]|J|d§dn 2 h —l-l [k4e] = 4 2 (Kt-1:4) 4 ( 3h ) J J 1B‘1 1Q11B 1IJId§dn .1.1 [14.1 = -(—Z"—"—‘—) J J 1B 1’1 1Q11B1IJId§dn [keg] = (Zk [k7e] = 23 32-]: 4 -1—1 ‘) J J [B I1 1Q11B'1IJIdédn -1-1 1 l (iii-3(1) J J 1B“ 1 1Q1IB‘“1IJId§dn —l—l [kge] = 5 5 5 1 ‘ (Zk Zk-l)( 322)] J[BH]T[§][BW]IJ|d§dn —l-l (3.57) (3.58) 52 3 3 129,1 = (—fi’—3"—‘3(-—) J J 18" ‘1 1Q11B‘1IJldédn -l—l 4 4 12.0.1: W(-—) J J1B"1“1Q11B‘111122§dn -l-l (Z5 ZS ) 211 [kne]__£._"_‘_1_(_h iJJ‘J'w -l-l 44111-341 12e " 6 72— (z‘i-z4 ) “(132] = k 412-1( “(142] = (22—05(4) —1—1 6 [kISe] = (Zk -6Zk 1)( h2 7 7 “(1661] = (Zk 7Zk_l)( 3: —1—1 The Jacobian matrix [J] is .95 91. 1J1 = a: a: .22. 22. .311 311 3.5.5 Consistent Mass Matrix =13?) 111 T— 111 ] [QllB lllldédn 75)]?er 1 [QllB V]|J|d§am 4 l 1 7‘ 3,12)! JIBIV] [QllBllllldgdn -1-1 1 l T 5 ( 3:2” Ile] [QlanllJldgdn "3X! 3h2 2” l [B l 1Q11B“‘1|J|d§dn --l l l l 2) J J [B‘V1T1Q113‘V1llldgdn (3.59) Mass and damping matrices are required for dynamic problems. Like the stiffness matrix, the consistent mass matrix can also be obtained by assembling the element mass 53 matrices defined by 1M.1 = J p1N1‘1N1dV (3.60) V where p is the mass density of the material, and [N] denotes the same interpolation func- tions used to generate the element stiffness matrix. The damping matrix, [C] is usually specified indirectly through modal damping ratios. The consistent mass matrix for the plate element is 1M.1 = J p1N1‘1N1dV 112k = 2 121. J J J 1N1‘1N1llldédndz k=1 -1-.1221 (3.61) 1%.]... 101 101 ‘ = 101 11142.18x8 101 _ 101 101 1M3,1,,,,,_ Performing the integral with respect 2 analytically, the sub-matrices in Eq. 3.61 can be ex- pressedas l 1 1M1 = 2p,.(2,.-z,._1)J J1N'1‘1N'1IJIdédn k= --1— l l l 1M21= 2 121(21- 21- l) J J 1N'1‘1N'1IJIdtdn (3.62) = -l- l N l 1 [M3]: “pg (2,—2.4) J J 1N'1‘1N'1IJIdédn = 44 where [N'] ___ N1 0 N2 0 N3 0 N4 0 (3.63) 0N10N20N30N4 54 The elements of [N’] are given in Eq. 3.33. The element of [N"] are §(l—&)(1-n)(2-§—n-§‘—n‘) 911—611 +11)(1—n)2 8 -§(1+1:)(1—n)(1—§)2 é“ +§)(1-n)(2+§-n-§2—n2) §11+§)(1+n)(1—n)" §(I-§)(1-n)(1+§)2 §11+§)(1+n)(2+§+n—§‘-n‘) _gu +§)(1—11)(1 +11)2 §(1—§)(1+n)(1+§)‘ §(1-§)(1+n)(2—§+n-§‘-112> —‘§‘(1-§)(1-11)(1+11)2 —§(1 +5311 +11)(1—§)2 3.6 Numerical Integration (3.64) (3.65) With isoparametric elements, explicit integration of the stiffness matrix is either intractable or lengthy and error prone for implementation in a computer program. There— fore, numerical integration is preferred, and the method particularly suited for finite ele- ment analysis is Gauss quadrature. In order to use Gauss quadrature, the integrand must be expressed in fi-n coordinates with the origin at the center of the element. In Gauss quadra- 55 ture, the integral is approximated by the summation of the weighted values of the inte- grand at n specified points. The various matrices in Eq. 3.58 become 1 1 J 1B11‘1Q11B21lJldédn '1“ (3.66) N :1 N n n 2 22 WW )IJI1B,(§.-. n ,)1‘1Q11Bz(§.-. 11 ,-)1 =11 j k 1K..1 = 2 k where W1 and W 1- are weighting factors for points i and j, fii is the location of point i rel- ative to the center and n is the number of points at which the integrands are to be evaluated. The order n of Gauss quadrature to be used in the evaluation of stiffness and mass matrices depends on the desired degree of accuracy. Exact integration requires the highest order of the integrand to be 2n — 1 . Since the integrands in Eq. 3.58 contain up to sixth order terms, = 4 yields exact results. However, n = 3 yields acceptable results. It should be noted that for theories that include transverse shear deformation, Gauss integration can sometimes produce erroneous results when the length-to—thickness ratio of the plate is large, i.e. the plate is thin. For a thin plate the transverse shear deforma- tion can be ignored. However Gauss integration can sometimes produce very stiff element stiffnesses. This phenomenon is called shear locking (Cook, et a1. 1989). From Eq. 3.56, the element stiffness matrix can be decomposed into a bending stiffness, [kb] and a trans- verse shear stiffness, [ks] , where [kb] is associated with in-plane deformation, while [k3] is associated with transverse shear deformation. With this notation, in static analysis, ([ka + [k,]){d} = {p} (3.67) where {d} represents the DOFs and {p} represents nodal loads. The deflections {(1} should be governed by [kb] for a thin plate. However, as a plate becomes thin, Gauss inte- gration can result in a [ks] that dominates the stiffnesses yielding displacements that are very small. Shear locking can be avoided by careful and consistent formulation of the in- 56 terpolation functions, but this is often tedious. A more convenient method of avoiding shear locking is to use reduced or selective Gauss integration when evaluating the stiffness ma- trix. For example, one-point Gauss quadrature should be used to evaluate [k s] for a four- noded rectangular element. Fig. 3.15 shows the normalized tip displacements of cantilevered plates with length-to-thickness aspect ratio, UH. The plate has the configuration show in Fig. 3.14 and point loads are applied at the end in the z—direction. The fiber orientation is 0° so no twisting occurs. The normalized tip displacement is obtained the tip displacement for a plate by dividing by the tip displacement computing using beam theory which includes shear deformation. For the loading used, beam theory with shear deformation should pre- dict the plate displacements accurately. It can be seen that reduced integration gives rea- sonable results for all UH values for FSDT and Reddy’s theory while full integration produces poor results, especially for large UH values. For Reddy’s theory, full integration and reduced integration yield similar results. Classical theory, which neglects shear defor- mation, yields acceptable results for UH values larger than about 20, but yields poor results for thick plates. 57 Figure 3.14 Cantilevered plate loaded in flexure 1.6 '- l.4 - 1.2 - Normalized Tip Displacements — Beam Theory —- CLPT ' — ' FSDT (2x2 integration) ------ FSDT (2x1 integration) —- Reddy (3x3 integration) - - Reddy (3x2 integration) 00 P l . l 1 l 5 10 15 Ratio (UH) 20 25 30 Figure 3.15 Variation of normalized tip displacements with a lcngth-to-thickness ratio of the plate CHAPTER 4 Equivalent Linearization Method 4.1 Introduction In recent years random vibration analysis (RVA) of structural systems is being increasingly used for a variety of applications. Many structures are subjected to excitations such as earthquakes, wind loads, sea waves, acoustic and aerodynamic loads etc, which are random in natures. RVA is an efficient technique for considering such excitations in the design of buildings, aircraft, offshore structures, etc. A complication that arises often is that many of these engineering structures exhibit nonlinear behavior. As the deformation increases, the nonlinearities arising from geometric considerations or material properties tend to significantly influence system behavior. Therefore, for accurate analysis, the analy- sis of nonlinear system under random excitations must be considered. In random vibration problems, nonlinear dynamic systems are generally repre— sented by nonlinear differential equations. Generally three methods are used to solve these equations: (1) the Fokker-Planck approach, (2) the perturbation approach, and (3) the equivalent (statistical) linearization approach. 0 The Fokker-Planck method is the only one that is exact, but its use is quite limited because of severe restrictions that must be placed on the form of the nonlinearities and on the spectral density matrix of the excitations (Caughey 1971). 0 In the perturbation approach, the stochastically excited nonlinear system is treated in the same way as a deterrninistically excited ones. The assumed solution which is represented as an expansion in powers of some parameter which characterizes the nonlinearity is substituted into the original equations of motion to yield a set of lin- ear differential equations. A first-order approximation is obtained by solving the linear systems. The perturbation method can be used only if the dynamic system has weak nonlinearities (Crandall 1963). 58 59 ° The equivalent linearization technique is the most versatile and widely used ap- proximation method. In this method the solution of the system of nonlinear equa- tions is approximated by the solution of an equivalent system of linear equations. The parameters of the equivalent linear system are obtained by minimizing an error measure. It has been shown that the method of equivalent linearization gives good results even for strongly nonlinear systems (Roberts and Spanos 1990). This meth- od is employed here. 4.2 Equivalent Linear Equation of Motion Consider a general multi-degree of freedom nonlinear system which is described by the equation of motion [M]{ii}+[C]{u} +[K]{u} +{¢(i2',12, u)} = {PM} (41) where [M] is the mass matrix, [C] is the damping matrix, [K] is the linear stiffness ma- trix, and {¢(ii, u, u)} is a nonlinear internal force vector. The vectors {u}, {u} and {ii} are the generalized displacement, velocity and acceleration vectors, respectively, and {P(t)} is a stationary Gaussian random excitation vector with zero mean. When the func- tion {¢(1‘¢', 11, u)} is non-zero and nonlinear, it is difficult to obtain the exact solution of Eq. 4.1 and approximate analysis must usually be used. The equivalent linearization method is a versatile and readily mechanized method for obtaining an approximation solution of Eq. 4.1. In the process of obtaining an approximate solution, the nonlinear equation of motion in Eq. 4.1 is replaced by an equivalent linear equation of motion for which the exact solu- tion can be obtained. The replacement is optimized with respect to the difference between the original and the equivalent system of equations. The equivalent linear system that replaces Eq. 4.1 is ([M] + [M*]){ii} + ([C] + [C*]){11} + ([K] + [K*]){u} = {P(t)} (4-2) 60 where [M *] , [C*] and [K*] are additional mass, damping and stiffness matrices, respec- tively. The difference {e} between Eq. 4.1 and 4.2 is (Spanos and Iwan 1978) {e} = [M1{ii}+[C]{d}+[K]{u}+{¢(ii, 11, 14)} -([Ml+[M*1){fi}-([C]+[C*]){11}-([K]+[K*]){u} (4.3) = {$07.11,u)}-[M*]{ii}-[C*l{li}-[K*1{u} [M *1 , [C*] and [K*] are selected such that the magnitude of the error {e} is minimized. That is T 2 2 2 . . E[{e} {e}] = E[e1+e2+ ...+en] = mrnrmum (4.4) where E denotes the statistical expectation operator and e], e2,...en are the element of {e} . Since the operator E is linear and e,-2 is positive, Eq. 4.4 requires that D? = E{e,-2} = minimum (4.5) FromEq.4.3 2 2 Di = E{ei} 2 (4.6) = “{[i— 2 (MIMI-+6; u +k:-juj):l} i = l,2,...,n j=l where q) is the ith element of {‘1’} and m and k; are the (i, j) elements of the ma- ij’ 01'1" trices [M*], [C*] and [K*], respectively. Therefore D depends on m and k; and Ij’ cIj’ the necessary conditions for Eq. 4.6 to be true are a . ——7(D,2)=O j=1,2,...,n "‘7' *(D?) = 0 j = 1,2, ...,n (4.7) J __,_f.‘(D ) = o j=1,2,...,n 3k” '1 a aci- 8 61 The first equation in Eq. 4.7 may be expanded to 37(1)?) = 1.45141) = 424-381} m-- am.. am.. U U U = ZE{[¢1" 2 (mrsiis + C13“: + kfsus):|(-fij)} 3 =1 (4.8) n It: at: at: 3:1 n It It t = —E{¢,.uj}+ 2 [misE{iisfij}+cisE{zierj}+kisE{usiiJ-}] = 0 5:1 which yields E{¢,-iij} = Z [miszsafi “Larisa-j} + k:SE{usiij}] (4.9) 5:] Similarly, expanding the other two equations in Eq. 4.7 yield E{¢i11j} = 2 [asap-2,2,.) +c;E{usuj} +k:sE{usriJ-}] (4.10) 3:] and E{¢iuj} = z [m:5E{iiqu-} +c:sE{riqu-} +k;E{uqu-}] (4.11) s=l Eqs. 4.9 to 4.11 can be rewritten more compactly as {k;}‘ E[¢i{z‘4}] = E[{a}{a}‘] {c:}T 1: l,2,...,n (4.12) :1: T _{mi}- where 62 {21:112},121,1211‘ (4.13) and {m:}, {47} and {k2} are the ith rows ofthe matrices [M*], [C*] and [K*] , respec- tively. 4.2.1 Gaussian Approximation Assuming that the excitation is Gaussian allows the expectations in Eq. 4.12, and hence the elements of [M *] , [C*] and [K *] , to be evaluated in terms of the correlation matrix of the random vector {ii}. In order to accomplish this, the following property of Gaussian vectors is used: Consider a single valued function of n variables (10) = 6101.112, y") (4.14) where the y is a jointly Gaussian random vector with zero mean and q( y) is sufficiently smooth so that first partial derivatives with respect to y,-, i = 1, 2, ..., n , exist. If the re- II lation, Iq(y)l < Aexp[ Z yj], a < 2 for some arbitrary A and any y is satisfied, then 1' = 1 E1yq1y)1 = EIyy‘1B1Vq1y)1 (4.15) where V is the gradient operator defined by a a a T V = —,..,—— 4.16 [an 3);), 3),] ( ) Assuming that {4)} follows the condition given above and utilizing Eq. 4.15, then the left hand side of Eq. 4.12 becomes E[¢,{t¢}] = E[{&}{&}T]E[V¢i] i = 1,2, ...,n (4.17) Substituting Eq. 4.17 into Eq. 4.12 yields 63 ',.71 {k1} E1{12}{12}‘1 {6:}7 =E[{12}{12}‘]E[V¢,.] 1' l,2,...,n (4.18) :1: T d”"i}_ Since {i1} is Gaussian, its covariance matrix will be positive definite. When the covariance matrix is positive definite, the solution of Eq. 4.18 is _ _ ’31),“ {123‘ 8714—} * T = E V. =E< 3‘11“ '=1,2,..., 4.19 {Ci} 1 (11,1 301} i l n ( ) {mZ}‘ 3912;. ’ . 19W. Thus the element of the matrices [M *] , [C*] and [K *] are given by ‘ — E a“ 420 ”'0' - a—JJ— ( . ) k‘ - E 31),. 422 ,7 - {)7}. ( . ) Eq. 4.15 was originally developed by Kazakov (1965). Atalik and Utku (1976) used it for the first time in random vibration analysis and also gave a more detailed proof of above theorem. Subsequently Spanos (1980) pointed out that it is also applicable for non-stationary random vibration analysis. For the laminated composite plates dealt with here, the nonlinear force vector depends only on the displacements, thus the nonlinearity * * will affect only the stiffness matrix and consequently, mij = cij = 0. Therefore the equivalent linear system equation is [Ml{ii}+[Cl{d}+([Kl+[K*l){u} = {P(t)} (4-23) 4.3 Nonlinear Element Stiffness Matrices From the conservation of energy principle as applied to finite elements, the virtual work done by external nodal forces when the nodes undergo virtual displacements is equal to the virtual strain energy of internal stresses; i.e., SW = 8U (4.24) where 5W is the virtual work of external forces on the element and 8U is the virtual strain energy of internal stresses. To develop these quantities, small virtual displacements, Sue , may be assumed. Therefore, the external work of nodal loads becomes SW = {Sue}T{Re} (4.25) where {R e} are the nodal loads on the element. The internal virtual strain energy is 511 = J {ae}‘{o}dv (4.26) V Substituting of Eq. 4.25 and 4.26 into 4.24 and using the strain-displacement relation, {62:}‘ = {5a,}‘131‘ yields {Sue}T{Re} = J {Sue}T[B]T{o}dV (4.27) V Using the constitutive law in Eq. 3.10 for {0'} and the strain-displacement relation in Eq. 3.56 65 {62.}‘{R.1 = J162.)‘1BI‘11Q1121+111“1f11T1“{e}1dv V = {811,}T'J131713ne12v.» J1B1‘1T1“1f1IT1-T{e}dv] 'V V (4.28) = {82.1"J1B1‘1Q11B1{u.}dv 'V + J1B1‘1T1“1f11T1“1B1{u.}dVJ V rn whrch [f] = [drag(0, 0, f6(712), f4(723). f5('Y]3))] . Srnce {8ue} 1s an arbitrary vrr- tual displacement, the following condition must hold for Eq. 4.28 to be true: {R.} = [JIB1‘1Q11B1BV]{u.} V + J 1B1‘1T1"1f1111“1B1{u.}dV (429) V [Ke]{ue}+{¢e} where [K e] is the element stiffness matrix given by Eq. 3.58 and {the} is the elemental nonlinear force vector. More explicitly, We} J 1B1‘1T1“1f1111“1B1{u.}dv V 2 J IB1‘1T1"1f11T1“IBI{u.}dxdydz A k=l (4.30) Zk—l Using the functional forms for the elements of [F] given in the nonlinear force vector in Eq. 4.30 can be rewritten as n Zr {1).} = 2 J J 1B1‘m“If11T1"1B1{u.}dxdydz kzlAZk—r z (4.31) J 11¢..}+{42.}+{¢3.}1{u.1dxdydz = A Z11-1 k 1 66 where {41.1 = 2 26.113181‘111“ 1412210, 0, 1, 0, 0)11T1“1B1 14.32) i=1 {1..1=22.,-1§§1B1‘IT1‘Idiag10 0 o 1 0)11T1‘IBI 14.33) i=1 {4...} = 2 a..vi§1BI‘1T1“1diag10.o,0,0, 1)11T1“1B1 14.34) i=1 Since the nonlinear vector is only a function of the displacement vector, only an equivalent stiffness matrix results in the method of equivalent linearization, and is [BS] - E[a——{ 311). 1] 21in E[§{i—J({¢1e}{ue})]dxdydz N ‘* a (4.35) 4);?!“ lE[m({¢2e}{ue})]dxdydz N Zr +31 E1a{ ——}-1{¢3.1{u. })]dxdydz Consider the first term on the right-hand side (RHS) of Eq. 4.35. Using Eq. 4.32 it can be expressed as Zk N 2:14.21 E[fia;e_}({¢1e}{ue})]dxclydz bl (4.36) Zk N n =2 22164 J [B] ""[T 11311.4 ]dxdydz = A “=1 111 where 67 [7:] = [71“1diag(o,o, 1,0,0)1[T1“ (4.37) and B 21' 14.1 = E1517’1“‘2‘“e“1 e 311 (4.38) = E[figuzg + 211128173; 157131] with [128] being a unit matrix of size 28 by 28. Since 712 = [T3]{8} = [T3][B]{ue}. where [T3] is a row matrix consisting of the third row of [T]_T , Eq. 4.38 becomes 14,1 = E1117.11B1{u.1)“111+2i{u.11IT311B1{u.1)“"17311B11 = E1117311B1{u.1)“1111+E1221u.11IT.1IB1{u.1)“" = a.111+2i{B,117311B1 17.11311 (439) where a. = E1117.11B1{u.1)“1 {11.1 = B112.111T.11B1{u.1)““1 (4,40) _ 1 11 2 _4 III 3 __4_ 1v [B1 - [B]+z[B 1+2(h2)1B 1+z( 3h2)[B 1 and [8‘] through [BIV] are given by Eq. 3.41 to 3.52. ()11 may be expanded as 68 “=1 I 5 28 5 28 5 28 = 2 2 2 .. 2 E[“e,k,“e,k2-~“e,k2,] J'1=“‘1=1}'2=”‘2=l I'm-“‘21: (441) 21 21 ' X01731.) H 31.1.] = :1 5 8 5 28 5 28 21 21 = 2 2 2 2.- 2 22s in (1173111112....) j,=111,=1j,=1k2=1‘j 2,.=Ik.,.=1 =1 m=l =1 where “e,k, denotes the klth element of {ue}. {Bl} maybeexpanded as 21-1 {Bl}=EL{ue}[q2282T3131kfluek) 5] j-lkzl 28 =ud2222 iiwwmmm Jr=lk1=112=21k2=1 121-1=]k21—1 =1 H 731.]21—1 31.12”] (4'42) :1 =1 5 28 5 28 5 28 = 2 2 z 2 Z 2 E[{u€}(u¢2k1ue k2 'uevk21—1)] j1=lkr=1j2=1k2=l 1.21.1'411‘214:l 21 21 X[H 731...] II 31.14.] m=l :1 The pth element of {[3]} is 5 28 5 28 1312-2222 2 iElueMek “2121] jr=lk1=1j2=1k2=1 121-1=1k21—1=1 (443) 1.2; 1 1:122) 69 Using the same procedure for the second and third terms on the RHS of Eq. 4.35 yields Q N EJ 1 E[§T%({¢2e}{ue}):|dxdydz k=l ~ ‘*" Z (4.44) = 222214,.J J [B]T[T;][B][A2]dxdydz k: 11:] A, 75—1 1 and N Zr 6 kg]! j E[m({¢3e}{ue})]dxdydz ’ ‘*“ (4.45) N n z. = 2 2215, J J [B] T[T;][B][A3]dxdydz k=1i=l 42.- 1 where 17;] = [T]-l[diag(0,0,0, 1,0)11T1“ (4.46) 11;] = [71"[diag(o, 0,0, 0, mm“ (4.47) [A2] = 012[I]+2i{132}[T4][B] (4.48) [A3] = a3[1]+2i{B3}[T5][B] (4.49) 5 8 5 28 5 28 2i 2 2 2 2 - 22112....J l'IT41.JJTl J 1450) Jr=lkr=1j2=1k2=l 12:11‘2=1 l 70 5 28 5 5 28 {112} = J 2 2 2 2 2 E[{ue}(ue,klue’kz...ue,k2i_l)] 1 1.21-1= 1k21—1=1 21 21 (4'52) XLH 741.1 H 31.8..” =1 m=l 5 28 5 28 5 28 {[33} = { Z 2 2 Z 2 z E[{ue}(ue,klue,k2°”ue,kz,-_l)] jr=lkr=1j2=1k2=1 j21-1=1k21-1=1 (453) x 2:) 2%) and [T4] and [T5] are row matrices consisting of the fourth and fifth rows of [Tl-T re- spectively. The expectations in Eq. 4.41, 4.42, 4.50, 4.51, 4.52 and 4.53 are 2i-th order mo- ments. Since it is assumed that {u e} is a zero mean Gaussian vector, the even order moments may be expressed in terms of the second order moments by the following formu- 1a: E[ue’klue,k2...ue,k2i] = 2E[ue,klue,k2]E[ue,k3ue,k4]"‘E[ue, k21-1ue.k2.-] (4.54) where the summation is taken over all possible ways of dividing the 21' variables into 1' com- binations of pairs. The number of terms in the summation is l.3.5...(2i-3)(2i-1). When formulating the linear element stiffness matrix the integrand was expanded so that the integration in the z-direction could be performed analytically while the integra- tions in the 5 and 1] direction were performed numerically. However for the equivalent ele- ment stiffness matrix, expanding the integrand so that the integration in the z-direction can be performed analytically requires lengthy and cumbersome algebraic operations that are prone to error. Therefore, numerical 3-point Gauss quadrature is used in all three direc- tions and Eq. 4.35 is written in the form IK: 1 = E[Mi Jam] EMz: -l-l—l ‘—-.I-‘ ...— 1112 la4i EM: — I H H 11,2 1112 EM: -——-111 The Jacobian matrix [J] is —aO 71 l l l a... j J j [B]T[T:][B][A1]|J|d§dndC l j [BlrszllBllAzlllldédndC l l 1 as,- j J I[B]T[T;][81[A3]|J|d§dndC 0 [J]: Ob O zk—Zk-l _00 ————2 _ (4.55) (4.56) where a and b are the dimension of the rectangular plane element of size 2a x 2b in Fig. 3.9. 4.4 Random Vibration Analysis Random vibration analysis to solve Eq. 4.23 is performed using an iterative * approach with each iteration consisting of a linear analysis. The matrix [K ] in any given iteration is computed using the nodal displacement covariances from the previous iteration and the iterations are terminated when the covariances converge. The stationary response of multi-degree-of-freedom (MDOF) systems to random excitations is well-known and is only summarized here. Time domain or frequency domain techniques can be used for the analysis. The main steps in a frequency domain approach are as follows: 1. Using the mass matrix [M] and the stiffness matrix [K] + [K*] where [K *] = O in the very first iteration, determine undamped natural frequencies, (oj , and mode shape matrices, {(1) j} , for a chosen number of modes. Usually 0)}. and {¢j} are ob- 72 tained by solving the free vibration problem. The free vibration of a MDOF system is characterized by the generalized eigenvalue problem [[Kl—mleHWj} = {0} (4.57) whose solution yields (Di and {(1) j}. . Perform a linear random vibration analysis to determine the covariances of the nod- al displacement {u} . If the excitation P(t) is taken to be a stationary random pro- cess, the covariances of the stationary response may be computed through E[urus] = X X 241%; X j Hj(—w)Hk(m)S,m((o)dco (4.58) j=lk=1 J I=1m=1-.., where ¢rj are rjth elements of the mode shape matrix, M j = {oj}T[M]{¢J-} is the jth modal mass. H 1(0)) is the modal frequency response function for mode j given by H 1(a)) = 1 (4.59) 2 2 . in which co}. and C}- = C j/ (2 W) are the undamped natural frequency and damping ratio of mode j, respectively. H j(—0)) is the complex conjugated of H 1(a)) . S ,m((o) is the cross spectral density function of the nodal force excitations P, and Pm at the lth and mth DOF. Note that for synchronous loading only the auto spectra exist while the cross spectra are zero, and hence the double summation in Eq. 4.58 over I and m reduces to a single summation. For certain types of excitation spectra, closed form solutions can be used to rapidly compute the integral in Eq. 4.58 (Harichandran 1992). Use of closed-form solutions instead of numerical inte- gration yields considerable computational savings. However for more general exci- tations numerical integration will need to be used. 73 3. Compute the equivalent element stiffness matrices [K :] in Eq. 4.55 and assemble the global equivalent stiffness matrix [K it] . The three steps outlined above are repeated until convergence is obtained in the covariances of nodal displacements. It is convenient to check for convergence by using the nodal displacement variances and the mth iteration is assumed to have converged if 2 Jz(cui.m _ Gum - l) l 2 cum 1' < 3 (4.60) in which oui = ./E[u,.2]. 4.5 RMS Strains and Stresses in a Laminate The determination of the root-mean square (r.m.s) stresses and strains is essential for the analysis of laminated composites. The laminate is assumed to consist of several perfectly bonded laminae; i.e., the displacements are continuous through the laminate thickness so that no slip occurs between each lamina. The method of computing the r.m.s. laminae stresses and strains from the displacement covariances is presented in this section. 4.5.1 Computation of RMS Strains Consider the strain-nodal displacement relation given by Eq. 3.53: {a} = [IB‘] + lenl + zl-filwml + 23("':72)[Bw]]{"‘} (4.61) = [B]{ue} Since the nodal displacements are zero mean Gaussian processes, the mean value of the strains are E{8} = E[[B]{ue}] = [BlElue] = 0 (4.62) The covariance matrix of the {8} = [8x9 8)” ny’ sz’ sz [is] = 74 element strains in global T. ] 1S — , Elsi] 514.2,] Elm.) Ewan] £12,7le Elexeyl Elsi] EI2,~I,,I Elem.) £12,)le Elsa.) EIe,v,,I mi) E[v.,v,,1 E[nyvle 512.7,.) 512.7,.) E[v.,v,,1 1317;] Ema.) _Eiexvle £12,.)le Emmy.) Elvyzvle £173.] _ and may be computed through [2.1 E[({€}-E{€})({8}—E{8})T] E[{8}{6}T] E[[Bl{ue}{u,}T[BlT] [B]E[{u,}{ue}T][B]T = [BlliuellBlT where the covariance matrix of the element nodal displacements is [2...] = E[uf] E[ulvl] E[uluz] E[u1(—w4’x)] E[ulvl] E[vh E[vluzl E[v1(—w4’x)] E[uluZ] E[vluz] E[ui] E[u2(-w4’x)] E[ulw4’y] E[vlw4,y] E[‘u2w4,y] E[w4’y(—w4’x)] _Elu1(-w4,x)lE[v1(-w4,x)]E[u2(-w4,,)l E[(—w4,,)21 . q coordinates, (4.63) (4.64) (4.65) Using the coordinate transformation from global and material coordinates for strains, 75 {2'} = [Tl—TR} _T (4.66) [T] [B]{ue} . T . . . . . where {e } = ‘ {81, 82, 712, 723, 713} , the covariance matrix of element strains in material coordinates is 1 Elsi] £12,221 121817.21 £18,723] E[emn Elelezl Elsi] Elez‘m] Elezvzgl 51227131 lze'l= ElelvlzlElezvlzl Ely?) E[712723lEl'Y12'YI3] (467) E[8I723] E[82723]E[712723] E[‘Ygsl E[723713] _Elem3] Elszvlgl E[Ynm] 517237.31 Elvigl _ and may be obtained as [2...] = E[{E'}{€'}Tl _ —T T —T T — E[[TTl {8}{8}T([T] l) l (4.68) = [T]- E[{8}{8} ][T]- = [Tl'TlBlliucllBlTlTl'l 4.5.2 Computation of RMS Stresses In material coordinates, the stresses at any point are given by {0"} = [Q]{E'} + [diag(0, 0: f6(712), f4(723), f5(713))]{€'} (4-69) . T . . where {o } = {61, 62, 112, 1:23, 113} . It IS convenient to separate the stress components into normal and shear stresses: 76 {0'} = {O'N} +{T'} ol 0 02 0 (4.70) = 1 0 >+( T12> 0 I23 (0 , TI3 The covariance matrix of element stresses in material coordinate is computed through [2...] = E[{O'}{6'}Tl = E[[{6'N}+{T'}][{6'~}+{T'}]Tl = E[{o'N}{o'N}T+{t'}{T'}T+{O'N}{T'}T+{T'}{0'N}T] = E[{0'~}{0'N}Tl+E[{1:'}{T'}T] +E[{6'~}{I'}Tl+E[{I'}{6'N}T] = [2),] +Iz..I+I2,.N..1 + [20,217 (4.71) where I— — E[of] E[oloz] E[oltlz] E[oltz3] E[oltn] E[oloz] E[ofi] E[oztlz] 51621523] E[oer] [20'] = “011121 5152112] E[Tizl “7121231 E[T12113] (4'72) 2 E[01723] E[02123] “7121231 H1523] E[Tz3113] 2 _EWITB] “521131 E[712T13] E[Tz3113] E[Tn] . Consider the first term on the RHS of Eq. 4.71. Using the normal stress-strain relations, 77 lo‘ ' 1reW ' Q11Q12000 1 { } 0'2 Q12 QZZOOO 82 O" = >= 4 1 N )0 0 o 000 0 (4.73) O O O 000 0 101 _0 0 ooqLoJ = [Q11/HEW} the covariance matrix of normal stresses in material coordinate can be computed as [20le = E[{6'~}{6'~}Tl E[[QN1{€'N}{8'N}T[QN]T] T T (4.74) = [QN]E[{8'N}{8'~} llQNl = [QNHZe'NIIQNIT where _ 2 ., E[el] E[elez] 0 O O 2 E[e 8 ] E[e ] O O O [2,. l = 1 2 2 (4.75) N O O O 0 O 0 0 O O O 0 o o o o_ In expanded form, the second term on the RHS of Eq. 4.71 is .0 o o 0 o 1 O O O O O 2 [21.] = 0 0 Elm] E[T12123] E[T12713] (4.76) 2 0051le123] 51123] E[TZ3T13] 2 _0 o E[ranI E[tz3113] 511131 _ The shear stress-strain relation may be expressed as 78 ' 0 ‘ ‘0 0 0 0 o _ r 0 ‘ o o o o o 0 0 {1'} = 1712) = 0 0 Q66+f6(712) 0 O (7121 (4.77) 123 0 O 0 Q44+f 4(1'23) 0 723 3:131 .0 0 O 0 Q55+f5(713)_ (713, and each term in Eq. 4.76 can be evaluated individually. The first non-zero diagonal term is 2 E[(Q66712 "' 2 “61722+ 1] ] '=‘ (4.78) h 2 m 2' = Q66E[712] "' Q66 2 “6151712+ 1"” 2 061'2 ‘161'15‘[('Y12+1)('Y12+ 1)] i=1 i=1 j=1 1311le where 06,- are the parameters used to approximate the shear stress-strain law as described in Section 3.2. Similarly, the other two non-zero diagonal terms are E11331 = 2241511331)» Q44 2 (14.19113? 21+ 2 a..- 2 a,,EI(v§§* 511%“ ‘11 (4.79) i=1 i=1 j=l 511131 = 1235511131»: Q55 2 afiElvifl + 2 a5,- 2 4,5111%? 511?“ ‘11 (4.80) i=1 i=1 j=l Since the shear strains are zero-mean Gaussian random variables, expectations of their even powers that appear in Eqs. 4.78 to 4.80 can be expanded using the formula in Eq. 4.54 as E[yfi‘] = 1x3x5x ...x(2k—3)x(2k—1)xE[yf2]" (4.81) 21 2 k E[723] = 1x3x5x...x(2k—3)x(2k—1)><104 stec indicating that the nonlinearity 0.1 m B ll l Figure 5.8 Three-ply symmetric laminated plate loaded in flexure 91 TABLE 5.3 First five natural frequencies from linear (first iteration) and nonlinear (last iteration) analysis for the excitation intensity S0 = 100,000 N23 Mode [30°/0°/30°] [60°/0°/60°] First Iteration Last Iteration First Iteration Last Iteration (H2) (H2) (HZ) (Hz) 1 86.75 84.97 60.42 59.63 2 334.65 331.66 266.06 260.47 3 495.48 490.56 387.49 381.26 4 51 1.16 493.95 443.68 440.02 5 947.63 934.98 852.70 840.25 is more significant for CL = 60 . Fig. 5.9 shows the variation of the absolute r.m.s shear strains in material coordi- nates at the center of element 2 of the plate with the excitation level. In-plane shear strains are computed at the top surface of the plate while out-of-plane shear strains are measured TABLE 5.4 Number of iteration required for convergence Load level 30(103N2 - sec) No. of Iteration [30°/O°/30°] [60°/0°l60°] 10 20 30 40 50 60 70 80 90 100 ##-F-§-§UJUJWUJUJ hhehehbmmw A1- 9“ 11er 311 300 exp lat .231 i I ”IE I '3? Lt: Jim‘ 4 t I ‘11.. hi”. 92 0.0055 r I r 1 ' 1 — 712 fora=302 1 0-005 —- 7‘3 fora=30 Myra '—' 723 fora=300 _.-""‘- /-:l 0.0045 ...... .712 fora=600 , 0 / m —' 713 for a = 600 " .g OW "" 723 fora=60 - i.. 3 0.0035 L I I I J I ’ /_‘ ,2 0.003 - /// . U) U) 0.0025 ’4 5 0.002 5 0.0015 - ’ I .9 i / 4 0.001 V’;::/""' 1 omos , 1 l 1 l g l 1 i q 100 150 200 250 300 Excitation Level (S0)u2 (N seem) Figure 5.9 Variation of absolute r.m.s shear strains with excitation level for symmetric three-ply at mid-depth where they are maximum. At S0 = 105stec , the r.m.s in-plane strain is about 0.005 for the [60°/0°l60°] ply and 0.0043 for the [30°/0°/30°] ply and hence the expected maximum shear strains are 3 x 0.005 = 0.015 for the [60°/0°l60°] ply and 3 x 0.0043 = 0.013 for the [30°/0°l30°] ply. Since the approximate shear stress-strain law is valid for shear strains up to about 0.02, loads of up to S0 = 105 stec can be ana- lyzed. For any given excitation level, the shear strains for a = 300 are smaller than the corresponding shear strains for on = 60°. The variation of the normalized r.m.s. z-displacement at the free right comer nodes of the plate with the excitation load level is shown in Fig. 5.10. The displacement is nor- malized by dividing by the linear response which is obtained at the end of the first itera- tion. As with the single-ply plate, the figure illustrates the increase in the displacement due to material nonlinearity as the excitation level is increased. The results show that nonlin- 93 earity increases the displacement by 2.5% for the [60°/0°/60°] ply and 2.2% for the [30°/ 0°/30°] ply compared to the linear response at S0 = 105stec . The effect of nonlinearity is therefore less significant for the three-ply case than the single-ply case. The variation of the normalized r.m.s. normal and shear stresses in material direc- tions at the center of element 2 of the plate with the excitation level is shown in Figs. 5.11 and 5.12. The normalization is performed by dividing by the corresponding linear responses. Fig. 5.11 indicates that for in-plane stresses nonlinearity is more significant for the [30°/0°/30°] ply than for the [60°/0°/60°] ply. This behavior is similar to that of the single ply plate. Fig. 5.12 shows that the in-plane r.m.s. shear stress for the [60°/0°/60°] ply is most affected by nonlinearity, with a 15% decrease compared to the linear response at S0 = 105stec. 1 .025 l .02 1.015 1.01 Normalized RMS Displacement 1 ‘0 r l L l 1 l 1 l 100 150 200 250 300 Excitation Level (30)"2 (N. seem) Figure 5.10 Variation of normalized r.m.s displacement with excitation level for symmetric three- ply l .045 1.04 1.03 Normalized RMS Normal Stresses 94 — Unfora=300 —- Unfora=300 —° (Inforat=600 °°°°° Uzzfora=600 1.035 1 .025 1.02 ”—‘C‘ 111. 1 Excitation Level (80)“ 2 (N secm) Figure 5.11 Variation of normalized r.m.s. normal stresses with excitation level for symmetric three- ply 1.1 j , i I . I L w ———-—— a) ____....—- #### § 1.05 " I”/’ .q 5, ,,,/z a wf”::' _____________________________ ,2 1.0 F'—".:."'_ ____________ - m ~~ ______________________ g ........................ . '8 0.95 .............. - % T12 fora=300 .............. . E 7'13 fora=300 ......... c 0.9 -' 1'23 foraz=300 ......... _ Z """" 1'12 forat=600 ......... _' T13 fora=600 - - 1'23 fora=60o 0.85 I I 4 4 Excitation Level (S0)ll 2 (N seem) Figure 5.12 Variation of normalized r.m.s. shear stresses with excitation level for symmetric three-ply 95 1.06 l i I ——- e” for a = 300 CD _" €22 fora=303 ’______.. I: 1.05 -—- enfora=60 /// q E """ 622 for a = 600 // m // a _ / - E 1.04 /// O . Z /// g 1.03 L / ‘l “D g 1.02 E 0 Z 1.01 / 1.0 100 . . l/ 2 ll2 Exc1tat1on Level (SO) (N sec ) Figure 5.13 Variation of normalized r.m.s normal strains with excitation level for symmetric three-ply 1.12 0 r I r r l -— 712 for a = 30 —- 713 forat=300 ,’ E 1.1 --- 723f0f0=300 ,’ - .8 """ 712 fora: 600 / b —- r - 60° I ’ ' m 713 ora - 0 / /// 31.08 -- 7nfora=60 /// //7 _ .2 I //// /’ if, ,9/ /’ 1.06 - ,/’/ /' a E /’// /’ 'c ' /’ ‘ 0 . N _ / .1 a 1.04 // /. g /% /' """""""" o //‘ .............. '— " 1.02 - ' ---------- ,__..—-—--'""" Z /./‘:' / ”H."___.....:::'_:...—---—"""" .1 1.0 ”‘9‘““2' ' l A l 1 100 150 200 250 300 Excitation Level (80)" 2 (N sec 1,2) Figure 5.14 Variation of normalized r.m.s shear strains with excitation level for symmetric three- ply 96 0.1m 14y. ‘ _l Figure 5.15 Three-ply antisymmetric laminated plate loaded in flexure The variation of the normalized r.m.s. normal and shear strains in material direc- tions at the center of element 2 of the plate with the excitation level are shown in Figs. 5.13 and 5.14. The behavior of the normalized in-plane r.m.s. strains are similar to that of the normalized in-plane r.m.s. stresses. 8“ and £22 for the [30°/0°/30°] ply are the most sensitive to nonlinearity. Fig. 5.14 indicates that shear strains are significantly affected by nonlinearity. 5.4.2 Antisymmetric Stacking Sequence The three—ply cantilevered plate with fiber orientation 0t, 0 and —0t in the top, mid- dle and bottom layers, respectively, as shown in Fig. 5.15 is examined. TWO values of at, 30° and 60°,were used and the plate was loaded at the free end in the z-direction. The total thickness of the plate is the same as that of the plates considered in Section 5.3 and Sec- tion 5.4.1. The magnitude of the excitation, So, was varied from 10,000 stec to 100,000 N2 sec with 10,000 stec intervals. Table 5.5 shows the first five undamped nat- 97 ural frequencies of the plate from the first and last iterations of the analysis, for the excita— tion level S0 = 100, 000 stec. The values from the first iteration correspond to the case where nonlinearity is neglected, while the values from the last iteration include the effect of nonlinearity. Due to nonlinearity, the natural frequencies decrease by about 1%, which is smaller than the corresponding decrease for symmetric plates. The ply arrange- ment with at = 300 is stiffer than the arrangement with CL = 60° , yielding higher natural frequencies. The number of iterations required for convergence of the r.m.s nodal displace- ments are listed in Table 5.6, and increase with the excitation load level since the nonlin- earity becomes more pronounced. Unlike the symmetric ply, the number of iterations required for convergence for the the ply arrangement with (1 = 30 is equal to or more than the corresponding number with a = 600 for SO between 3x104stec to 6x104stec. This indicates that the nonlinearity is more significant for the ply arrange- ment with on = 30 . TABLE 5.5 First five natural frequencies from linear (first iteration) and nonlinear (last iteration) analysis for the excitation intensity S0 = 100,000 N25 Mode [30°/0°l-30°] , [60°/0°l—60°] First Iteration Last Iteration First Iteration Last Iteration (Hz) (Hz) (Hz) (Hz) 1 77.31 74.96 59.26 58.61 2 336.73 333.76 290.77 286.77 3 453.71 446.06 380.57 375.07 4 561.93 557.96 446.91 445.1 1 5 1039.00 1029.03 929.58 919.79 98 TABLE 5.6 Number of iteration required for convergence Load level No. of Iteration So ( 1031.12 . sec) [30°/0°/-30°] [60°/O°/-60°] 10 3 3 20 3 3 30 4 3 40 4 3 50 4 3 60 4 3 70 4 4 8O 4 4 90 4 4 100 4 4 Fig. 5.16 shows the variation of the absolute r.m.s shear strains in material coordi- nates at the center of element 2 with the excitation level. In-plane shear strains are com- puted at the top surface of the plate while out-of—plane shear strains are computed at the mid-depth where they are maximum. For any given excitation level, the in-plane shear strains for at = 30. are higher than the corresponding shear strains for a = 600 , while the out-of-plane shear strains for or = 300 are smaller than the corresponding shear 0 strains for 0t = 60 . The variation of the normalized r.m.s. z-displacement at the free right corner nodes of the plate with the excitation level is shown in Fig. 5.17. The displacement is normalized by dividing by the linear response The results show that the effect of nonlinearity increases the tip displacement by 4% for the [30°/0°l-30°] ply and 2% for the [60°/0°/- 99 0.006 r f . r . r — 7,2 forat=300 —- 7,3 forat=300 0.005 -—- mfora=30° ------ -- 0 ------ 7,2 forat=600 m —- 7,3 forat=600 ..... /. ' 50004 -- 723fora=60 /' fi . _ /. / b /' /// J V) /' // ’ ’ ' a ‘ ...«-— .«- .. "' 2 0.003 ///a: ’ ’ ’ .... .1 m ._."' ’ ’ ’ ’ 5 0.002 ......... / : .. .". / ’ L": /',/’,/’ ’ ’ ____._.. ~— -" "' " firfx” ,,._.._...-—-——-""' 0.001 E ..... .... ---- _ 00 4 l . 1 i J 4 l 100 150 200 250 300 Excitation Level (80)”2 (N seem) Figure 5.16 Variation of absolute r.m.s shear strains with excitation level for antisymmetric three- ply l .045 1.04 141;! l .035 l .03 l .025 1.02 Normalized RMS Displacement 100 150 200 250 300 Excitation Level (80)”2 (N. seem) Figure 5.17 Variation of normalized r.m.s displacement with excitation level for antisymmetric three-ply 100 60°] ply relative to the linear response at S0 = 105stec . This indicates that the effect of nonlinearity is more pronounced for the antisymmetric ply than for the symmetric lami- nates cases. The variation of the normalized r.m.s. normal and shear stresses in material direc- tions at the center of element 2 with the excitation level is shown in Fig. 5.18 and 5.19. The normalization is performed by dividing each stress by the corresponding stress from linear system. The results show that the effect of nonlinearity on the transverse normal r.m.s. stress for the [30°/0°/-30°] ply increases them by about 4.5% from the linear case at SO = 105stec. For longitudinal normal r.m.s. stresses, about a 2.5% increase is indi- cated. Fig. 5.19 shows that the in-plane r.m.s. shear stress for the [30°/0°l-30°] and [60°/ 0°/—60°] plies are the most affected by nonlinearity. 112 decreases by about 15% from the . 5 linear responses at S0 = 10 stec. T I ' I ' — allfora=30° // i 1.04 —- (Izzfora=300 // _ 0 --- Unfora=60 / o / 1.035 °°°° 022 for a = 60 // _ l .03 1 .025 Normalized RMS Normal Stresses 8 Excitation Level (80)” 2 (N seem) Figure 5.18 Variation of normalized r.m.s normal stresses with excitation level for antisymmetric three-ply 1.0 0.95 0.9 Normalized RMS Shear Stresses 0.85 101 _ — ‘ _—__ ——-—- ~._ -—~ '- — — o . — 1'12 fora=302 —- 7,3 fora=30 -—- rz3fora=300 ------ rlzfora=600 —- *rl3fora=600 -- rz3fora=60 1 I r A 100 1 50 200 250 Figure 5.19 Variation of normalized r.m.s shear stresses with excitation level for antisymmetric 1 .045 _. I-t _ o e o 9 8 fl E —o C N M Normalized RMS Normal Strains E e Excitation Level (80)” 2 (N seem) three-ply — e” forat=300 —- ezzfora=300 °—’ 6“ fora=600 ------ ezzfora=60 l "T 1 Y r V LL14. Excitation Level (80)” 2 (N seem) Figure 5.20 Variation of normalized r.m.s normal strains with excitation level for antisymmetric three-ply 102 The variation of the normalized r.m.s. normal and shear strains in material direc- tions at the center of element 2 with the excitation level are shown in Figs. 5.20 and 5.21. e“ and £22 for the [30°/0°l-30°] ply are most affected by nonlinearity. The out—of-plane shear strains shown in Fig. 5.21 display similar behavior for the [30°/0°l-30°] and [60°/0°l ~60°] plies. The 712 response for the [60°/0°/-60°] is least affected by nonlinearity. 5.5 Comparisons with Lower Order Theories Laminated plates have been analyzed by several theories, which mainly differ in the way they account for shear deformation. Two theories, which have been widely applied to plates, are classical laminated plate theory (CLPT) and first-order shear defor- mation theory (FSDT). The classical laminated plate theory based on the Kirchhoff assumption neglects transverse shear deformation, while the first-order shear deformation theory is the lowest-order improvement to the classical theory and yields constant values 1.07 r r r l ' l — 712 fora=300 /, —- 713 for a = 300 // a 1'06 -—- 723 forci=300 /’ '1 'a ------ rmfora=603 , ,« 51) 1.05 :‘L 713 foratf600 ,/ .. g 723 for a — 60 / ‘ ..i: V) 1.04 U) E 1.03 'o o .13. g 1.02 o z 1.01 1.0 1 r n 1 100 150 200 250 300 . . l/ 2 1/2 Excrtatron Level (So) (N sec ) Figure 5.21 Variation of normalized r.m.s shear strains with excitation level for antisymmetric three-ply 103 of transverse shear strains through the thickness of the plate. Results from these simple and popular theories are compared to results from Reddy’s third-order theory considered so far. 5.5.1 Symmetric three-ply laminated plates The response of the symmetric three-ply cantilevered plate with the same geome- try as shown in Fig. 5.8 was computed using the lower order theories. Two values of a, 30° and 60° were used and the plate was loaded at the free end in the z-direction. The exci- tation level So was increased from 10,000 stec to 100,000 stec in 10,000 stec inter- vals. Table 5.7 shows the first five undamped natural frequencies of the plate from the last iteration of the analysis at the excitation level S0 = 100, 000 stec for the first-order shear deformation theory (FSDT), classical laminated plate theory (CLPT) and Reddy’s theory. The natural frequencies resulting from the CLPT are larger than those resulting from the FSDT, which in turn are larger than those resulting from Reddy’s theory. The dif- ferences between CLPT and Reddy’s theory are as high as 32% for some frequencies. The difference in natural frequencies between the FSDT and Reddy’s theory are within 10%. Inclusion of shear deformation through progressively higher order theories therefore yields more flexible models. TABLE 5.7 Natural frequencies predicted by the three theories-Symmetric three-ply Mode [30°/0°l30°] [60°/0°/60°] Reddy FSDT CLPT Reddy FSDT CLPT l 84.97 84.66 88.02 59.63 59.15 60.09 2 331.66 335.56 376.79 260.47 273.64 299.20 3 490.56 493.75 493.85 381.26 358.52 377.06 4 493.95 548.36 654.29 440.02 441.84 441.92 5 934.98 939.84 1 106.18 840.25 872.58 993.46 104 Fig. 5.22 to 5.24 show the variation of the absolute r.m.s shear strains in material coordinates at the center of element 2 of the plate with the excitation level. Fig. 5.22 shows that Reddy’s theory produces the highest 712 values for (1 = 60° and the lowest for = 30°. Figs. 5.23 and 5.24 indicate that Reddy’s theory yields significantly higher transverse shear strains than the FSDT. The CLPT neglects transverse shear deformation and hence Figs. 5.23 and 5.24 do not contain transverse shear strains for the CLPT. For any given excitation level within any theories, the shear strains for the ply arrangement with or = 30 are smaller than the corresponding shear strains for the ply arrangement with or = 60 . 0.005 . I . I - l r /r . / I / zfflfi 0.0045 - //z ’ . ' / , I, /,/ 0'... o/ 0.004 r .. U) E 0.0035 — - m 5 0.003 - - .e: W i ‘0 0.0025 F - 052 /./-’ — 712 for a= 30° (Reddy) 0.002 . —- 712fora= 30° (FSDT) -—- 7,2fora= 30: (CLPT) -- rm fora: 60: (Reddy) 0'00” - - 712 fora: 60: (FSDT) ' ------ rm fora: 60o (CLPT) owl r l 1 1 i 100 150 200 250 300 Excitation Level (80)”2 (N secm) Figure 5.22 Variation of absolute r.m.s. 712 shear strain with excitation level for symmetric three-ply 105 0.005 r r l r r —— 7,3 fora: 30: (Reddy) . 0.0045 '—' 7,3 fora: 30: (FSDT) //’- —- 7,3 fora: 60: (Reddy) // . 0,004 """" 7,3 fora: 60O (FSDT) // q a ‘ / . a 0.0035 g 0.003 :3 . o g 0.0025 m . 5 0.002 0.0015 0.001 0.0005 P ' ‘ 4 1 A - I 100 150 200 250 300 Excitation Level (80)"2 (N seem) Figure 5.23 Variation of absolute r.m.s. 713 shear strain with excitation level for symmetric three-ply 0.004 r r r ' l — 723 fora: 30: (Reddy) . '_' 723 fora: 30: (FSDT) _ 00°35 —- 723 fora 60° (Reddy) /. ----- 723 fora: 60o (FSDT) // ‘ 9 § RMS Shear Strains a .. a 3 § 8 0.001 Excitation Level (80)“2 (N seem) Figure 5.24 Variation of absolute r.m.s. 723 shear strain with excitation level for symmetric three-ply 106 1.03 . , . , , — a=300(Reddy) —- a=300(FSDT) E 1-02 -—- a=300(CLPl‘) - GE) — a=600 (Reddy) J 3 - - a=60°(FSDT) 3 1.01 ------ a=60°(CLPT) - E l --------------------- '8 099 _ ........................................................................................................ 4 .5 —- —- —- —- —- —- —- —- —- —- —- —. _- _- _ '3 E 0.98 - _ O z _________________ 0.97 - ———————————————————————— r r l . l 1 l . l 100 150 200 250 300 1/2 Excitation Level (80)”2 (N sec ) Figure 5.25 Variation of normalized r.m.s. displacement with excitation level for symmetric three-ply The variations of the normalized r.m.s. z-displacement at the free right corner node of the plate with the excitation load level as predicted by the different theories are shown in Fig. 5 .25. The normalization has been performed by dividing each response by the non- linear response resulting from Reddy’s theory at the same load level. The figure indicates that the displacement from the CLPT is about 3% less than that from Reddy’s theory for (1 = 30°. At any given load level the displacements increase with the order of the theory, i.e. CLPT —) FSDT —-> Reddy. The variation of the normalized r.m.s. normal stresses in material directions at the center of element 2 of the plate with the excitation level as predicted by the different theo- ries are shown in Figs. 5.26 and 5.27. The behavior of the stresses is opposite to that of the displacements with the normal stresses resulting from CLPT and FSDT being larger than those resulting from Reddy’s theory. _o t—I . . r. u—n u—n — N Ih _ o 00 § Normalized RMS Normal Stresses '8 § 1'" O 100 Figure 5.26 Variation of normalized r.m.s. a” with excitation level for symmetric three-ply i—nu—I r—o—n NAGOO fl 0 ~ ~ Normalized RMS Normal Stresses '2 '8 8 H __ . ._ 'o 8 8 Figure 5.27 Variation of normalized r.m.s. 0'22 with excitation level for symmetric three-ply 107 a” for a: 30: (Reddy) 0“ for a: 30: (FSDT) - a” for a: 302(CL171') a“ for a: 60: (Reddy) a“ for a: 60: (FSDT) ------ on for a: 600 (CLPT) p—o—o—u—u— I—n—e-r— — ‘ '—.—I—-o—n_o—a_.—.—o—o— —.— . . o u— _ D Excitation Level (80)” 250 2 112 (N sec ) 022 for a: 30: (Reddy) 022 for a: 30: (FSDT) ' 022 for a: 30: (CLPT) on for a: 60:: (Reddy) 022 for a: 60: (FSDT) """ 022 for a: 600 (Clam I .—O‘ — - - _0~ .h . ~. ‘- ~. ~I‘ . _C‘. ‘- ~o—o_ _ — —'- . .—.—.—I 0 150 Excitation Level (80)” 200 250 2 m (N sec ) 108 The variation of the normalized r.m.s. shear stresses in material directions from the different theories at the center of element 2 of the plate with the excitation level are shown in Figs. 5.28, 5.29 and 5.30. Fig. 5.28 indicate that the normalized in-plane shear stress for 300 for the FSDT and CLPT is larger than that for Reddy’s theory, while for (1 CL = 600 the reverse is true. Fig. 5.29 and 5.30 show that the transverse shear stresses pre- dicted by FSDT are significantly lower than those predicted by Reddy’s theory. This is because the FSDT yields constant transverse shear stresses having approximately average values. Fig. 5.31 and 5.32 show that the variation of r.m.s. transverse shear stress through the thickness for S0 = 100, 000 st. Again it is apprant that the constant stresses in each layer predicted by FSDT are significantly smaller than the maximum value predicted by Reddy’s theory. 1.06 4 , L , , r p I I.“ - ~- ~' -——.- -I .—.—. —'—- ~O‘. ~g~ ._"-‘—— — o§.—- ‘ . ~u—-—. -.~I—. _o—n- . ‘ 1.0 J — 7,2 for a = 30: (Reddy) —- m for a = 30 (FSDT) 0 ' ' 1'12 for a = 300 (Clam L — 7,2 for a = 600 (Reddy) 0-9 - - 7,2 for a = 60 (FSDT) .o \o 4; .1 Normalized RMS Shear Stresses O O :3 8 l l l l I I I I I I I I I I I | L °°°°°° 7,; for a = 600 (CLPT) 088 1 l i l i 100 150 200 250 300 Excitation Level (80)“ 2 (N seem) Figure 5.28 Variation of normalized r.m.s. 1:12 with excitation level for symmetric three-ply 109 F I l 1-1 — 1'13 for a = 300 (Reddy) ‘ m -—- 7,3 for a = 30: (FSDT) a 1.05 — 7'13 for a = 600 (Reddy) - 3 """ 1,3 for a = 60 (FSD'D l: m 1.0 ... 5 r .c: 0.95 - _ m m o 5 0 9 - // a '8 / 1 N 0.85 - / ‘ fl * ./ ,/ g 0.8 — _/,/ 1 2 .,.,/ ............................. l 0'75 p -r-' . _—: ; 31.37.1171: ------------------ ‘1 Lm\\fl’""’ “.______._.. m). ..... ‘ 0.7 1 t 1 - . 100 150 200 250 300 Excitation Level (80)“ 2 (N seem) Figure 5.29 Variation of normalized r.m.s. 113 with excitation level for symmetric three-ply if I ' l 1.2 — 723 fora=300(Reddy) - -—- 12, fora: 30°(FSDT) , § 1.1 —— Tz3fora=600(Reddy) _ 8 """" 7’23 fora=600 (FSDT) l: m 1.0 _ is ’ . -== - _ m 0.9 m b 1 E 08 - - g 0.7 ~ . E 0.6 ~ - 2 t 0.5 L—’“'_"""——°—'—'"'—'—'—'-'-‘-°-—--—-—-—-—.—._.:.._ L ................................................................................................... 0.4 I 1 1 L n 100 150 200 250 300 Excitation Level (80)” 2 (N secm) Figure 5.30 Variation of normalized r.m.s. 1:23 with excitation level for symmetric three-ply 110 2.0 - r I r a , f -- n3 fora=30o (Reddy) 1.5 — — 1,3fora=300(FSDT) 1.0 ~ Thickness .0 O .5 LII .1. 'o -l.5 - _20 A L a . A L - . A J . . l O 2 4 6 8 10 12 l4 l6 RMS Shear Stresses for Three-Ply (MPa) Figure 5.31 Variation of RMS ‘I:l3 through the thickness for symmetric three-ply — 723 for a = 300 (Reddy) - — 723 for a = 30° (FSDT) l I ‘8 I E I .0 0.0 : E . -05 . I f——" I -1.0 ~ I < l -I.5 - . I . a I . A a . - . m . - -2.0 0 2 4 6 8 10 12 14 RMS Shear Stresses for Three-Ply (MPa) Figure 5.32 Variation of RMS 1:23 through the thickness for symmetric three-ply 111 5.5.2 Antisymmettic three-ply laminated plates The response of the antisymmetric three-ply cantilevered plate as shown in Fig. 5.15 was computed using the lower order theories. Two values of 0t, 30° and 60°, were used and the plate was loaded at the free end in the z-direction. The excitation level So was increased from 10,000 stec to 100,000 stec in 10,000 stec intervals. Table 5.8 shows the first five undamped natural frequencies of the plate from the last iteration of the analysis at the excitation level S0 = 100, 000 stec for the CLPT, FSDT and Reddy’s theory. The predicted natural frequencies are larger for the CLPT than for Reddy’s theory and FSDT. The differences between CLPT and Reddy’s theory for ct = 300 are as high as 22%. However the ply arrangement with CL = 60° produces predicted frequencies that differ by less than 10% for Reddy’s theory and CLPT. In general, the differences in natural frequencies predicted by Reddy’s theory and CLPT are larger for the symmetric three-ply than for the antisymmetric three-ply. As with the symmetric three-ply, the FSDT predicts frequencies that fall between those predicted by Reddy’s theory and CLPT. The variation of the absolute r.m.s shear strains with excitation level for the anti- symmetric case is similar to the corresponding variation for the symmetric case. At any TABLE 5.8 Natural frequencies predicted by the three theories-antisymmetric three-ply Mode [30°/0°/-30°] [60°/0°l-60°] Reddy FSDT CLPT Reddy FSDT CLPT 1 74.96 75.84 77.30 58.61 58.47 59.18 2 333.76 357.95 409.24 286.77 296.35 320.18 3 446.06 463.94 494.70 375.07 352.07 368.78 4 557.96 548.38 580.1 1 445.1 1 443.99 448.30 5 1029.03 1 102.35 1324.55 919.79 930.96 1060.66 112 given excitation level the r.m.s shear strains from Reddy’s theory are larger than those from CLPT and FSDT. The variations of the normalized r.m.s. z-displacements from the different theories also has a pattern similar to that obtained for the symmetric case. The variations of the normalized r.m.s. normal stresses in material directions from the different theories at the center of element 2 of the plate with the excitation level were examined. The difference between the responses for the symmetric and antisymmetric cases is that for the latter, the normalized responses are not much different from Reddy’s theory and FSDT for or = 300 and or = 600. However, the responses predicted by CLPT is significantly different from those predicted by the other two theories. For example, the normalized 011 for (1 = 600 is larger for the antisymmetric case than for or. = 600 in the symmetric case. As with the symmetric ply, all shear stresses predicted by FSDT and CLPT are smaller than those predicted by Reddy’s theory. 5.6 Effects on 2-3 Plane Isotropy Assumption From experimental results, the shear stress-strain response of a unidirectional lam- ina is clearly nonlinear, while the longitudinal and transverse stress-strain responses of most unidirectional laminae are linear or only mildly nonlinear. Different models for the nonlinear stress-strain response of composite laminae have been developed by different researchers. However, most models are usually limited to two-dimensional cases that con- sider only the in-plane shear stress-strain behavior. The Hahn and Tsai (1973) nonlinear model, which is also limited to two dimensions, has been used in this study, and for the 3- D case, the out-of-plane shear stress-strain relations were assumed to have the same form as the in-plane shear stress-strain relation. In thick FRP laminae, the fibers are usually ran- domly packed in the 2-3 plane as shown in Fig. 5.33. Therefore, it is resonable to assume that there is no preferred direction in the 2-3 plane. By considering the 2-3 plane to be iso- tropic, the transverse normal shear modulus of the lamina may be expressed as: 113 3 / ML . 151i ............ III-II II uuuuuuuuuuuuu IIIII I .................. II nnnnnnnnnnnnnnnnnn I uuuuuuuuuuuuuuuuuu II oooooooooooooooooo -- q. .I ‘I ‘I I I I I I I I I I I I I I I I I I I I I .......... . I C I I Figure 5.33 Material direction with 2-3 plane of isotropy 52 623 = 2(1+v23) (5.1) where 623 , E 2 and v23 are the transverse-normal shear modulus, transverse Young’s modulus and the Poisson’s ration in transverse normal (2-3) plane, respectively. Since the transverse stress-strain relation between 0'2 and 82 is almost linear or mildly nonlinear for most FRPs (e. g. boron/epoxy is mildly nonlinear and a graphite/epoxy is almost linear), it is assumed that E 2 is constant. From Eq. 5.1, if the Poisson’s ratio, v23 , is constant as com- monly assumed, 023 is also constant. Therefore, the transverse-normal shear stress-strain relation between 023 and 723 will be linear under these assumptions. Experimental results have also shown that the shear stress-strain law in the transverse normal plane is almost lin- ear or very mildly nonlinear for certain types of FRPs and ply arrangements (Gipple and Hoyns 1994). Their experiments were limited to unidirectional single lamina with 0° fiber orientation and symmetric balanced cross-ply laminae with [0°/90°]s fiber directions for graphite/epoxy FRP. The effect of assuming linear behavior in the 2-3 plane for a single- ply laminae is considered in this section. The constitutive equation for the normal longitu— 114 dinal (1-3 plane) shear stress—strain relation is assumed to have the same nonlinear form as the constitutive equation for the in-plane (1-2 plane) shear stress-strain relation. The trans- verse normal (2-3 plane) shear stress-strain relation is assumed to be linear. It should be noted that this assumption may be appropriate only for unidirectional lamina and laminae with balanced symmetric cross-plies, since no data appear to available for other cases. Nev- ertheless, this section compares responses computed assuming 2-3 plane isotropy (linear behavior in the 2-3 plane) with responses computed assuming nonlinear shear behavior in the 2-3 plane for general plies. 5.6.1 Single-Ply Plate with 2-3 Plane Isotropy A single-ply cantilevered plate with the same geometry and fiber orientations shown in Fig. 5.1 is examined with and without shear nonlinearity in the 2-3 plane. Fig. 5.34 and 5.35 show the variation of the absolute r.m.s shear strains in material coordinates at the center of the element 2 of the plate with the excitation level. The differences result- ing from the assumptions of linear vs. nonlinear shear behavior in the 2-3 plane are small and are not significant for the 713 responses. The difference is more noticeable for the [30°] ply than for the [60°] ply. 115 Excitation Level (S0)U2 (N secm) 0.005 r ' T ' — 7,2 for a: 30: — — 1,3 for 01:30: . "_’ 723 for a: 30: J 0004 ------ 712 for a: 30: (2 3 Plane Isotropy) w —- 7,3 for a: 30: (2- 3 Plane Isotropy) . .§ _ _ 723 for a: 300 (2- 3 Plane Isotropy) // / i / - % 01KB ///'/ /- / d / ‘ / s a? ’ .5 m m 0.002 A E ‘g/ “’ " = J 0001 {/’/ '.—-“_r_‘_,,_=-z-=‘-'- ...,. . ash—ud— ‘"‘ M, . 4"” .J’. "I’. F'” 0.0 1 L 1 l ‘ l l 100 150 200 250 300 . . 1/2 1/2 Exc1tation Level (So) (N sec ) Figure 5.34 Variation of absolute r.m.s. shear strains with excitation level for single-ply 0.007 r ' ‘ ' -— 7,2 for a: 60: _- 723 for a: 60: 0-006 -—- 7,3 for (1:60: - ...... .712 for a: 60: (2- 3 Plane Isotropy) —- for a: 60: (2- 3 Plane Isotropy) 00 723 g 0 5 - - 7,3 for a: 600 (2- 3 Plane Isotropy) j 1:! m 0.004 d a . a" .8 f" m 0.003 z/’ "‘5‘- m / Ins-’55— a "” ‘ 0.002 , ,— ‘ ’/ ’ ,-/',——/’ I / ’ I ’ 0001 =5," ‘ 0.0 t 1 1 l . l l 100 150 200 250 300 Figure 5.35 Variation of absolute r.m.s. shear strains with excitation level for single-ply 116 The variation of the normalized r.m.s. z—displacement at the free right corner nodes of the plate with the excitation load level is shown in Fig. 5.36. Again there is no signifi- cant differences between the 2-3 plane linear and nonlinear cases. The variation of the normalized r.m.s. normal and shear stresses in material direc— tions at the center of element 2 of the plate with the excitation level are shown in Figs. 5.37 to 5.40. Differences between the 2-3 plane linear and nonlinear cases are noticeable for the [30°] ply but less significant for the [60°] ply. This limited study indicates that for single ply FRP plates, it may be reasonable to assume linear behavior in the 2-3 plane. Computational cost is significantly reduced by this assumption. 1.06 0 I ' I T I : a i 300 (Shear) // ‘ a - 60 (Shear) / E 1,05 *—' a = 300 (2-3 Plane Isotropy) / - a ------ a = 600 (2-3 Plane Isotropy) // 0 is On .2 Q E 8 '—3 E 2 1.0 I l . I A I . l 100 150 200 250 300 Excitation Level (30)"2 (N. seem) Figure 5.36 Variation of normalized r.m.s. displacement with excitation level for single-ply l .06 l .05 1.04 l .03 Normalized RMS Normal Stresses 117 —— a” fora = 300 _- 022 for a = 300 -—- a“ for a = 300 (2-3 Plane Isotropy) ----- 022 for a = 300 (2-3 Plane Isotropy) b / p- g l l 200 250 Excitation Level (30)" 2 (N secm) l 50 300 Figure 5.37 Variation of normalized r.m.s. 0'“ with excitation level for single-ply .— O U) U! 2" O m I—I O N M § 1.015 1.01 Normalized RMS Normal Stresses l .005 — a” for a = 600 —- 022 for a = 60 - — - on for a = 600 (2-3 Plane Isotropy) ------ 022 for a = 600 (2-3 Plane Isotropy) l 250 Excitation Level (30)" 2 (N sec‘”) 150 200 Figure 5.38 Variation of normalized r.m.s. 022 with excitation level for single-ply 1.04 l .02 1.0 0.98 0.96 0.94 0.92 Normalized RMS Shear Stresses 0.9 0.88 118 b ~.‘ I I.. '- I ‘ '--§.~ "~Uc—I~ .~. ~P‘ . ‘.~O~O~I-. —-I—I I... c _— .— _— _— ..——— — .I— — _— .- _— _—————— ———- —.—---v.—. . — c I I I I I I I I I I I I I I I I I I I I I. I ' 7'23 fora = 300 7,2 for a = 300 7,3 for a = 30" m for a = 30° (2-3 Plane Isotropy) 0 m for a = 30 (2-3 Plane Isotropy) T23 for a e. 30° (2-3 Plane Isotropy) L 100 Figure 5.39 Variation of normalized r.m.s. shear stresses with excitation level for single-ply 1.02 W I ‘ l l ' 1,0 —- m=r——_' — - _—.:...—...' ='..—:.:_: :7: ':_-_ ———————————— q a: T~—_.§ —’ —"" “"‘°- '8. Er=-'-—_'_- 9 ‘“~~~ ----- i 8 0.93 L ‘ ~~ e * 1:1 ' - m 0.96 ‘ \ \ \ ‘ 5 “‘~- _c 0.94 _ U) . .I (D E 0.92 _ i '8 0.9 4 N a 0,33 -— 7'12 for a = 600 ‘ E —- 7,3fora=600 r o 0.86 —' 123 for a = 600 ' Z ----- 1'12 for a = 600 (2-3 Plane Isotropy) 0.84 —- n, for a = 60: (2-3 Plane Isotropy) ‘ - - 1'23 for a = 60 (2-3 Plane Isotropy) 0.82 ‘ ‘ 100 150 200 250 300 l 50 200 250 . . ll 2 1/2 Exc1tatlon Level (So) (N sec ) Excitation Level (80)” 2 (N secm) Figure 5.40 Variation of normalized r.m.s. shear stresses with excitation level for single—ply CHAPTER 6 Deterministic Nonlinear Dynamic Analysis 6.1 Introduction The only way to assess the accuracy of nonlinear random vibration analysis using the equivalent linearization method is through deterministic simulations of nonlinear dynamic response. Deterministic nonlinear dynamic analysis is often very costly and problems with stability and accuracy can be encountered especially if long responses are to be computed. A nonlinear dynamic analysis is usually performed using an incremental formulation, in which the variables such as displacements, accelerations, and velocities. are updated incrementally corresponding to time steps. An implicit unconditionally stable time integration scheme is usually employed for multi-DOF systems. The trapezoidal rule, more popularly known as the Newmark [3 method (1959), is suitable for the time integra- tion and offers unconditional stability. For acceptable accuracy, the iterative solution of the incremental equilibrium equations is required to assure equilibrium with iterations within each time step being terminated when appropriate convergence criteria are met. While the iterative scheme often yields sufficiently accurate solutions, difficulties such as slow con- vergence and divergence can sometimes occur. The most frequently used iteration scheme for the solution of nonlinear dynamic equations is Newton-Raphson iteration. The New- ton-Raphson iteration requires the tangent stiffness matrix [K T] to be formed and factor- ized at each iteration within each time step, and can be expensive for large-order systems. A modified Newton-Raphson iteration is often used to avoid the expensive repetitions of forming and factorizing the tangent stiffness matrix. This method differs from the New- ton-Raphson iteration only in that the tangent stiffness is either not updated or is updated infrequently. However, the modified method often needs more iteration to converge. In general it is efficient to update the stiffness matrix at the start of every time step and reuse it for all the equilibrium iterations within the time step (Bathe 1986). 119 120 6.2 Incremental Formulation of Equation of Motion The incremental equilibrium equations to be solved at time t + At in nonlinear dynamic analysis is [M]{ii}f'3m+[C]{d}filA,+[KTl,{Au}(i) = {R},.A,—{F}fi'1) (6.1) where [M], [C], and [K T] are the mass, damping and tangent stiffness matrices respec- tively, {Au} is the displacement increment, {R} is the externally applied load, {F} is a vector of nodal loads that correspond to the internal element stresses, and the superscripts (i) and (i — 1) represent quantities at equilibrium iteration i and i - 1 respectively. After convergence of the equilibrium iteration at each time step, the displacements are updated using {u},.A, = {u},+{Au} (6.2) where {u }, is the final converged displacements at time t and {Au} is the increment over time At. When a system has proportional damping, the mode superposition method can be used to obtain a damping matrix if a damping ratio can be specified for each mode. A pro- portional damping matrix can be constructed using [C] = [Mll¢l[diag(2if’)]I¢ITlMI (6.3) where [(0] is the mode shape matrix, M j is the jth modal mass, and C]- and to]- are the jth modal damping ratio and natural frequency, respectively. 6.2.1 Implicit Time Integration The incremental equilibrium equations given by Eq. 6.1 can be solved at each time step using a numerical time integration scheme. Various integration operators are used in practice. The trapezoidal rule, i.e. the Newmark B method (Newmark 1959) with constants B = 0.25 and Y = 0.5 is considered because of its relatively good stability and accuracy 121 characteristics. In the trapezoidal rule of time integration, the following assumptions are employed. 00.... = {a}.+%’<{a}.+ia}..a> A (6.4) tun... = {u}.+;‘<{u}.+{u}..a> Using Eqs. 6.4, Eq. 6.1 becomes (iIMl +-2—[C] + [K I){Au} = {R} - {F} Atz At T 1 HA: t+At 4 4 . .. - [MI(P({u},.A,— {“}r)'§{“}r‘ {u},) (6.5) —[C](A£t({u}t+At— {u}1) _ {u}t) It is solved for all time steps recursively. Although the time increment can be varied during the computation, a constant increment is typical. In general, the magnitude of the time step At is an important factor that affects the accuracy of the solution. Since it is often difficult to judge whether the chosen time step is small enough, iterations that assure equilibrium at the end of each time step are desirable. Equilibrium iteration procedures for the step-by- step solution of the equation of motion is described in Section 6.4. 6.2.2 Formation of Tangent Stiffness Matrix For the incremental formulation, the tangent stiffness matrix [KT] is required. For the constitutive law in Eq. 3.6 and Eq. 3.9, the incremental stress-strain relation in material coordinates is 122 _3___{G'} _3([Q1{€' 1) +3([f1{€'}) [CT]: a{e} a{e} 8{€} BIQ] 8{e'} alf] aw} =a—.—{8}{e'+} a—'_{c}[Q]+ a—{——8.'+}{e} a——{€.}[fl (6.6) _ 81f] —[0]+[Ql+a—{— .'}{e}+[f] in which [f]: [diag(0, 0, f6(712), f4('yz3), f5(713))]. Since [f] only contains shear components, each element of [CT] corresponding to shear 18 stated explicitly. at12 _ 8(Q66712) f)____(f..(rn))Y 9m 57;. “av— av—T . "no H4412” n . 21—1 21 = Q66 ‘1' 2 “61(20le (712) "' 2, “61'le (6.7) i=1 i=1 n . 2' = Qoo+ 2a6i(21+1)y,; i=1 Similarly 9.12 = Q44+ flag-(214 mg; (6-3) 8723 i=1 and 6113_Q+" 2'12'. 69 a— - 55 24511” )713 (.) 713 i: 1 Therefore [C 1] may be compactly written as [CT] = [Q] + [QN] (6.10) where [Q N] is the nonlinear part of [CT] in material coordinate, i.e. [QM]: [diag(0, 0, Z a6,(2i+ 1mg, 2 a4,.(2i+ 1172;, 2a,,(2i+1)yf§]] (6.11) i=1 i=1 i=1 123 From Eq. 6.10, the incremental constitutive law in material coordinates in Eq. 6.6 may be expressed as {46'} = [[Q1+[Q~]]{d€'} (6.12) The constitutive equation in global coordinate given by Eq. 3.10 may be extended as {o} = 1011s} + 1T1“1diag(f.(v.2>. f4(723), f5(713))1[T1'T{€} {Que}+[T:]{8}(f6(712))+[T;]{8}(f4(723))+[T;]{e}(f5(713)) (613) 1011s}+ITIIIe}(f.(IT311e})> + 11311210401411.2111 + [T;]{e}(f5([T5]{e})) where [Ti] is a matrix consisting of the ith row of [T]-T and [TI] = ITI“1diag(0.0. 1.0011111" 1T2] = 1T1“1diag10.0.o. 1. cilm‘T (614) 11;] = 1T1“1diag(0.0.0.0. 111m" The tangent constitutive matrix in global coordinates is - 3_{__G G} C l T]: 88{ } _ as . =1QI—E—1+ [T1]{8} a{f6([T3]{3})} a_{8} am +111 111161113112 21)}— 3{e } 8 [T ] 8 +[T2]{8} {f4(a{:}{ 1)} +[T;]{f4([T41{8})}a-{-—E} 3{s} 3{f5([T51{8})}+ + [T3]{8} 812} +[T;]{f5([T51{8 iii-3%; Recognizing that [T3]{8} = 712 , the derivative of the first nonlinear term on the RHS is £2! (6.15) 124 61120131121» _ 3{f6([T31{€})}3712 d{€} — a'Yiz a{€} _ a{f6(712)}a([T31{3}) - a'Y12 a{8} a n 21 = — 2061712)}[T3HI] {3712£= 1 (6.16) = [2 “61(2i)712_1][T3] i=1 = [ a6i(2i)([T31{3})2i—1)[T3] i=1 The other two nonlinear derivatives in Eq. 6.15 can be expressed in a similar manner, and thus [ET] = [Q] +1T111e11731(2 a6.(2i)(17311e}>2“‘)+11:1[2 a6.(1r3l{e})2‘] i=1 i=1 +1T§11211T41[2 a4,(2i)(1T.I{e})2“ ‘] + 1T;I[ 041([T4]{8})2i] (6.17) 1:] + [T;]{e}[T5][ 2 aSi(2i)([T5]{8})2i—1]+[T;][Z a5,-([T5]{e})2,-) . i=1 i=1 = [Eli-[EN] where [QN] is the nonlinear part of [ET] containing all but the first term on the RHS of Eq. 6.17. The incremental constitutive law in global coordinate is {dd} = [[§]+[§~]]{d8} (6.18) It can be shown that [EN] = [T]"1[QN][T]—T (see Appendix 13.). The tangent stiffness matrix of a finite element with N laminae is 125 N 1K.,TJ = Z [ IBIT11§I+IQ~IIIBIdv k=1V N N = Z Itmfiattmdw 2 [IBITIQNIIBIdV V k=l k=lV = [Kc] + 1K. N] (6.19) where [Kc] is the linear element stiffness matrix given in Eq. 3.56 and [K e, N] is the ma- trix accounting for nonlinearities. Using the relations 712 = [T3]{£} = [T3][B]{ue} 723 = [T4]{£} = [T4][B]{u¢} 1’13 = [T51{€} = [T5][B]{ue} and simplifying Eq. 6.17, [K e, N] may be expressed as N 1K.,~1 = 2 jiBl’Iévllmdv k=lV N = 2il{ki.,~}+{k2.,,v}+{k3e,N}1dV k = 1V where {kw} = [B]T[T:1[81{§: 2i+1)a.,-(1T31131{u }) 2"} 1an} = 131’1T3IiBi{i§n‘,(2i+1)a..(1T.I1BI{u1) } {k3e,~} = [B]T[T;][31{.1(2i+ 1)05.-([T5][B]{u D21} In terms of the natural E— n - i; coordinate system, Eq. 6.19 becomes (6.20) (6.21) (6.22) 126 N [Ke,N] = Z k=l N +21 kl- 1 l j lik..,~}IJId§dndt; Lb..— H b b.— C—u—‘I— k Jdd d { 2e,N}| 1 Es T1 (.1 (6.23) _- 1 fl 1 _I +k§l—Il-:[1 l j {k3.,~}|J|dE.dndC —l = [16“,] + [Kan/1 + [Kma N] in which each element of [K e, N] is evaluated using Gaussian quadrature in all three direc- tions. 6.2.3 Formation of Element Load Vector The element load vector {F e} in Eq. 6.1 contains the nodal loads that are consis- tent with the element internal stresses and is defined by N {Fe} = 2 J[B]T{o}dV (6.24) k=lV Using Eqs. 3.10 and 3.56, the element load vector {F e} may be expressed as N {F.} = 2 {[131’1011311u.}dv k=1V + j{BITITIIIBMu.}(f,(v.2))dv V + j[31711311811u.}(f.(vaa))dV (6.25) V + JIBITITQIIBI1u.})dv 2‘, {13171111 llB]{u 1(a,,y,;)dv (6.26) i=1 V =20 l IBITITIIIBIIu.}(1T.lIBI{u.})2"dv V and similarly, the other two terms are if). N} = 2 a4,j1BITIT;I[BI{u.}(1T.11BI{u.})2‘dV (6.27) i=1 n {f3.,N} = 2a5.leI 7,1T311311u }([T5][B]{u.})2idV (6.28) i=1 As with the stiffness matrix the integrations in Eq. 6.26 through Eq. 6.28 are performed nu- merically using three point Gauss quadrature. 6.3 Convergence Criteria Significant errors may be introduced into the solution of the incremental equation of motion in Eq. 6.1 if a time step marching procedure is used without correcting for devi- ations from the equilibrium state. The errors can become serious when a large time step is used and/or the system is highly nonlinear. The errors can result in instability of the solu- tion even though an integration operator that is unconditionally stable in linear analysis is used. In time history dynamic analysis, the solution at a specific time is dependent on the history of solutions and problems of slow convergence and divergence can be encountered for nonlinear systems, especially if the response is required for a long duration. The use of 128 equilibrium iterations at the end of each time step improves accuracy and appropriate con- vergence criteria are required to determine when the iterations should be terminated. Three measures are commonly used to check for the convergence: displacements, out-of- balance loads, and the incremental internal energy. 6.3.1 Convergence Criterion Based on Displacements Since the primary solution variable is displacement, it is easy to check the magni- tude of the displacement correction at the end of each iteration, and this commonly takes the form lliAui‘ll ————.— < tolerance (6.29) i llmwll where the superscript i denotes the quantity at the ith equilibrium iteration and the subscript t + At denotes the time at which the quantity is completed. Fig. 6.1 illustrates one dimen- {R}I+AI """"""""""""""""" ' """"""""" {F}(i) _______ {rim ---------- I {mm ----- I I {F}. : : : : 3 Ana» 3mm EAum: 1 <——N——>-<—v - A - {u}t {(1101) {u}(1){1:}(2) {Uigm Figure 6.1 Convergence criteria based on displacement increments and out-of-balance loads 129 sional equilibrium iteration for a static problem. Another way is to check convergence is to monitor the convergence factor = ———“{A“.}l” (6.30) i|{Au}"1” If the iteration is converging, q q2 > q3 > as is apparent in Fig. 6.1. However in general nonlinear analysis of multi-DOF systems the convergence factor may not behave in such a stable manner. 6.3.2 Convergence Criterion Based on Out-of-balance Load A more reliable convergence criterion is based on requiring the out-of—balance loads to be within a preset tolerance of the original load increment: out-of-balance load load increment = "1111....-mill..—1M11a1f11.-1611u}i:1.||(tolerance ||{R},,.A,-{F},-[M]{ii},-[C]{l1},ll (6.31) The out-of-balance load is shown in Fig. 6.1. The disadvantage in using an out-of-balance load check is that when the load vector contains different types of loads, such as forces and moments, inconsistencies in units can appear. 6.3.3 Convergence Criterion Based on Internal Energy In order to provide some indication of when both displacements and loads are near their equilibrium values, the increment of internal energy during each iteration, i.e. the amount of work done by the out-of-balance loads on the displacement increments, can be compared to the initial internal energy increment. This convergence criterion takes the form 130 {Aui}T({R}t+At"{F}i::3t‘[M]{a}i;lAt‘[C]{u}i;ist) l T at 1: . if .1 L L 1 l 1 l 1 I 1 I 1 J 4 l 1 I k l 1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Time (see) I I I I I I I I fi ------ Nonlinear ; -— Linear 5. .. MM 1 4 I 1 l 1 l 1 I 1 I 1 I 1 I 1 I 1 I 1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Time (see) Figure 6.2 Typical simulations of linear and nonlinear responses 135 Since the deterministic simulations contain a transient part at the beginning, r.m.s. values for comparison with stationary random vibration analysis results should be com- puted after truncating the transient portion. Fig. 6.3 shows the r.m.s. values computed from 20 nonlinear deterministic simulations after truncating progressively longer segments at the beginning. Truncation of segments longer than 0.45 seconds do not significantly affect the computed r.m.s. displacements indicating that stationary response is reached after about 0.45 seconds. Consequently, the first 0.45 seconds of the response was truncated for each linear and nonlinear deterministic simulation. In order to assess how many simula- tions are required to obtain stable r.m.s. values, the computed r.m.s. values based on the number of simulations used are shown in Fig. 6.4. The figure indicates that when more than about 8 simulations are used, the nonlinear r.m.s. displacement is larger than the lin- ear displacement. Within the limited number of simulations used, the r.m.s. displacements estimated from the simulations compare reasonably well with the results from the random vibration analysis. However, the fluctuation in the r.m.s. values even with 20 simulations indicates that a much larger number of simulations would be required to get reliable esti- mates. From the computational point of view, a single nonlinear random vibration analy- sis took approximately the same amount of computer time (about 6.6 hours on a Sun Ultra] Unix workstation) as six deterministic simulations at excitation level, 3x104stec . As the excitation level is increased, the computing time for the deterministic simulations increase because of the increased time for convergence iterations. For high excitation levels deterministic simulations can become unstable after some time. This is due to errors that build up over time due to stronger nonlinearity. The instability could perhaps be avoided by using a smaller time step or performing full New- ton-Raphson iterations during each equilibrium iteration. However, these refinements were not attempted owing to the significant computational cost. 136 0.017 0.0168 0.0166 0.0164 0.0162 0.016 RMS Displacement (m) 0.0158 I 0.0156 I 1 A 1 l — Simulation ------ RVA 0.0 0.1 0.2 0.3 0.4 Truncated Time (sec) 0.5 Figure 6.3 Variation of nonlinear r.m.s. displacement with truncated length of initial transient part T I r T f Y I r ' I 0.017 *- .1 ————————————————————— —1---—-*-I-—-1——‘7 E 0.016 _ v E o E 0015 E . 9.. fl Q d m 0.014 5 E a 0 013 _ ------ Nonlinear Simulation ' — Linear Simulation —- Nonlinear RVA °—- Linear RVA 0012 l 1 1 1 1 l 1 l 1 4 6 8 10 12 14 16 18 Number of Simulated Samples 20 Figure 6.4 Variation of linear and nonlinear r.m.s. displacements with number of simulations CHAPTER 7 Conclusions and Recommendations 7.1 Summary and Conclusions A random vibration analysis technique for laminated FRP plates modeled with finite elements and including material nonlinearity is presented. Since the shear stress- strain response of a FRP lamina is clearly nonlinear, this feature is considered in the anal- ysis. Nonlinear material models available are mostly limited to two-dimensional cases, i.e. to in-plane shear stress-strain behavior. A particular 2-D model proposed by Hahn and Tsai (1973) is used in this study and extended to 3-D based on some assumptions. Use of the shear strain-stress equation proposed by Hahn and Tsai is not tractable for the random vibration analysis, and therefore an approximate nonlinear shear stress-strain law expressed in terms of an odd-powered polynomial is used. Since transverse shear deformation is important for even moderately thick FRP plates, a high-order shear deformation theory proposed by Reddy (1984) which has cubic displacement fields, is adopted. The plate is discretized into finite elements with 4 nodes per element. Each node has 7 DOF, which include extensional deformations, shear defor- mations in the x and y directions, and flexural deformations in all three directions. Exact solutions of nonlinear random vibration problems can only be obtained for simple sys- tems. Therefore, an approximate nonlinear random vibration analysis of the laminated FRP plate is performed using the method of equivalent linearization. Equivalent lineariza- tion replaces the original nonlinear system by an equivalent linear system, and is capable of dealing with relatively complex structural models. The solutions can be obtained by using an iterative approach, where linear random vibration analysis is performed in each iteration. Within iterations, convergence of the solution is checked by using the nodal dis- placement variances. A computer program was developed and implemented for the non- linear random vibration analysis. 137 138 Several numerical examples are presented for single-ply and three-ply cantilevered FRP plates with rectangular geometry. All laminae are assumed to be made of Boron/ Epoxy Narmco 5505. Various types stacking sequences with fiber orientations of 0°, i300 and :60, are considered. Displacement, stress, and strain responses are computed. The results show that the effect of nonlinearity on the responses becomes more pro- nounced as the excitation level is raised. The number of iterations required for conver- gence also increases with the excitation level. The nonlinear shear stress-strain law produces significant nonlinearity in some dis- placement, strain, and stress responses, depending on the fiber orientations, stacking sequences, material properties and plate geometry. Ply arrangements with a fiber orienta- tion of 30° are stiffer than those with a fiber orientation of 60°. In terms of displacements, ply arrangements with laminae having 60° fiber orientations behave more nonlinearly than ply arrangements with laminae having 30° fiber orientations in symmetric laminated FRP plates, while in antisymmetric laminated FRP plates the reverse is true. The r.m.s. responses obtained from classical laminated plate theory (CLPT) and the first-order shear deformation theory (FSDT) are compared with the r.m.s. responses obtained from Reddy’s theory. The structural model becomes progressively more flexible as progressively higher-order shear theories are used. Of the three theories considered, Reddy’s theory yields the largest displacements while CLPT yields the smallest displace- ments. Transverse shear stresses predicted by FSDT are significantly smaller than those predicted by Reddy’s theory. CLPT neglects transverse shear stresses altogether. Deterministic simulations of nonlinear dynamic responses are performed for a can- tilevered plate with an antisymmetric ply arrangement and a fiber orientation of 30°. Lin- ear and nonlinear dynamic response are computed for a duration of 1 sec using a time step of 0.001 see. As expected, the amplitudes of nonlinear responses usually exceed the amplitudes of corresponding linear responses. The r.m.s. responses estimated from 20 139 simulated deterministic responses compare very well with the values computed through the approximate nonlinear random vibration method. A single nonlinear random vibration analysis took approximately the same amount of computer time (about 6.6 hours on a Sun Ultral Unix workstation) as six deterministic simulations at a moderate excitation level. As the excitation level is increased, the comput- ing time for the simulation increases because of slower convergence. Since hundreds of deterministic simulations are required to reliably estimate statistics of the random responses, the nonlinear random vibration analysis is much more cost effective than deter- ministic simulations. 7.2 Recommendations for Future Research Based on the study performed in this thesis, the following recommendations are suggested for future research on the nonlinear random vibration analysis of laminated FRP plates: 1. Global models such as those considered in this work are adequate for estimating global responses such as displacements and natural frequencies. However, if accu- rate interlaminar stresses are of interest then a more accurate local layerwise theory should be used. Layerwise theories consider each layer individually, and therefore the size of problems increase significantly. 2. In addition to nonlinear shear stress-strain behavior, some FRP lamina such as bo- ron/aluminum and carbon/carbon composites exhibit nonlinearity in normal stress- strain relations. The formulation used in this work needs to be extended to include this behavior. 3. The contribution of transverse normal stress and strain, 0, and 82 , is neglected in this work. The accuracy of the analysis can be increased by including these. 4. FRP laminates are usually used in either plate or shell configurations. The formula- tion of laminated shell elements are more complicated than that of laminated plate 140 elements, because of their curvature. The formulation should be extended to shell elements for more general applications. LIST OF REFERENCES Adams, D. F. (1970). "Inelastic analysis of a unidirectional composite subjected to trans- verse normal loading," Journal of Composite Materials, 4, 310-328. Amijima, S. and Adachi, T. (1979). "Nonlinear stress-strain response of laminated com- posites," Journal of Composite Materials, 13, 206-218. Atalik, T. S. and Utku, S. (1976). "Stochastic linearization of multi-degree-of-freedom non-linear systems," Earthquake Engineering and Structural Dynamics, 4, 411-420. Baber, T. T. and Noori, M. N. (1984). "Random vibration of pinching, hysteretic systems," Journal of Engineering Mechanics, ASCE, 110, 1036-1049. Bathe, K.-J. (1982). Finite element procedures in engineering analysis, Prentice-Hall, Englewood Cliffs. Bathe, K.-J. and Cimento, A. P. (1980). "Some practical procedures for the solution of nonlinear finite element equations," Computer Methods in Applied Mechanics and Engi— neering, 22, 59—85. Bathe, K.-J., Ozdemir, H. and Wilson, E. L. (1974). Static and dynamic geometric and material nonlinear analysis, Report UCSESM 74-4, Structural Engineering Lab., Univer- sity of California, Berkeley, California. Belytschko, T. (1976). "A survey of ’ numerical methods and computer programs for dynamic structural analysis," Nuclear Engineering and Design, 37(1), 23-34. Bendat, J. S. (1958). Principles and applications of random noise theory, Wiley, New York. Bendat, J. S. and Piersol, A. G. (1986). Random data: analysis and measurement proce- dures, 2nd ed., Wiley, New York. Caughey, T. K. (1959). "Response of a nonlinear string to random loading," Journal of Applied Mechanics, 26, 341-344. 141 142 Caughey, T. K. (1963). "Equivalent linearization techniques" Journal of the Acoustical Society of America, 35(11), 1706-171 1. Caughey, T. K. (1971). "Nonlinear theory of random vibrations," Advances In Applied Mechanics, 11, 209-253. Cederbaum, G., Elishakoff, I. and Librescu, L. (1989). "Random vibrations of laminated plates modeled within the first order shear deformation theory," Composite Structures, 12, 97-111. Cederbaum, G., Librescu, L. and Elishakoff, I. (1988). "Random vibration of laminated plates modeled within a high-order shear deformation theory," Journal of the Acoustical Society of America, 84(2), 660-666. Cook, R. D., Malkus, D. S. and Plesha, M. E. (1989). Concepts and applications of finite element analysis, 3rd ed., Wiley, New York. Crandall, S. T. (1963). "Perturbation techniques for random vibration of nonlinear sys- tems," Journal of the Acoustical Society of America, 35(11), 1700- 1705. DiSciuva, M. (1985). "Development of an anisotropic, multilayered, shear-deformable rectangular plate element," Computers and Structures, 21(4), 789-796. DiSciuva, M. (1987). "An improved shear-deformation theory for moderately thick multi- layered anisotropic shells and plates," Journal of Applied Mechanics, 54, 589-596. Dvorak, G. J. and Bahei-El-Din, Y. A. (1979). "Elastic-plastic behavior of fibrous compos- ites," Journal of the Mechanics and Physics of Solids, 27, 51-72. Foster, E. T. (1968). "Semi-linear random vibrations in discrete systems," Journal of Applied Mechanics, 35, 560-564. Foye, R. L. (1973). "Theoretical post-yielding behavior of composite laminates," Journal of Composite Materials, 7, 178-193. Gipple, K. L. and Hoyns, D. (1994). "Measurement of the out-of-plane shear response of thick section composite materials using the v-notched beam specimen," Journal of Com- posite Materials, 28, 543-572. Hahn, H. T. (1973). "Nonlinear behavior of laminated composites," Journal of Composite Materials, 7, 257-271. Hahn, H. T. and Tsai, S. W. (1973). "Nonlinear elastic behavior of unidirectional compos- ite laminae," Journal of Composite Materials, 7, 102-118. Harichandran, R. S. (1992). "Random vibration under propagating excitation: closed-form solutions," Journal of Engineering Mechanics, ASCE, 118(3), 575-586. Harichandran, R. S. and Hawwari, A. (1992). "Non-linear random vibration of filamentary composites," Computing Systems in Engineering, 3(1-4), 469-475. 143 Harichandran, R. S. and Naja, M. K. (1997). "Random vibration of laminated composite plates with material non-linearity," International Journal of Non-Linear Mechanics, 32(4), 707-720. Hashin, Z., Bagchi, D. and Rosen, W. (1974). "Non-linear behavior of fiber composites," NASA Contractor Report 2313. Iwan, W. D. and Mason, A. B. (1980). "Equivalent linearization for systems subjected to nonstationary random excitation," International Journal of Non-Linear Mechanics, 15, 71-82. Iwan, W. D. and Yang, I-M. (1972). "Application of statistical linearization techniques to nonlinear multidegree-of-freedom systems," Journal of Applied Mechanics, 39, 545-550. Iyenger, R. N. (1988). "Higher order linearization in nonlinear random vibration," Intema- tional Journal of Non-Linear Mechanics, 23(5/6), 385-391 . Izumi, M., Zaiming, L. and Kimura, M. (1989). "A stochastic linearization technique and its application to response analysis of nonlinear system based on weighted least-square minimization," Journal of Structural and Construction Engineering, 395, 72-81. Jones, R. M. (1975). Mechanics of composite materials, McGraw-Hill, New York. Jones, R. M. and Morgan, H. S. (1977). "Analysis of nonlinear stress-strain behavior fiber- reinforced composite materials," AIAA Journal, 15(12), 1669-1676. Jones, R. M. and Nelson, D. A. R (1975). "A new material model for the nonlinear biaxial behavior of ATJ-S graphite," Journal of Composite Materials, 9, 10-27. Jones, R. M. and Nelson, D. A. R (1975). "Further characteristics of a nonlinear material model for ATJ-S graphite," Journal of Composite Materials, 9, 251-265. Jones, R. M. and Nelson, D. A. R (1976a). "Material models for nonlinear deformation of graphite," AIAA Journal, 14(6), 709-717. Jones, R. M. and Nelson, D. A. R (1976b). "Theoretical-experimental correlation of mate- rial models for nonlinear deformation of graphite," AIAA Journal, 14(10), 1427-1435. Kapania, R. K. and Raciti, S. (1989). "Recent advances in analysis of laminated beams and plates, part 1: shear effects and buckling," AIAA Journal, 27(7), 923-934. Kazakov, I. E. (1965). "Generalization of the method of statistical linearization to multidi- mensional systems," Automation and Remote Control, 26, 1201-1206. Kenaga, D., Doyle, J. F. and Sun, C. T. (1987). "The characterization of boron/aluminum composite in the nonlinear range as an orthotr0pic elastic-plastic material," Journal of Composite Materials, 21, 516-531. Kimura, K. and Sakata, M. (1981). "Nonstationary response of a nonsymmetric nonlinear 144 system subjected to a wide class of random excitation," Journal of Sound and Vibration, 76(2), 261-271. Kuppusamy, T., Nanda, A. and Reddy, J. N. (1984). "Materially nonlinear analysis of lam- inated composite plates," Composite Structures, 2, 315-328. Librescu, L. (1986). "A general transverse shear deformation theory of anisotropic plates," In Refined Dynamical Theories of Beams, Plates and Shells and Their Application, Elish- koff, I. and Irretier, H. (eds), Springer-Verlag, Berlin, 32-43. Lo, K. H., Christensen, R. M. and Wu, E. M. (1977). "A high-order theory of plate defor- mation, part 2: laminated plates," Journal of Applied Mechanics, 44, 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, 59, 502-509. Mei, C. and Prasad, C. B. (1989). "Effects of large deflection and transverse shear on response of rectangular symmetric composite laminates subjected to acoustic excitation," Journal of Composite Materials, 23, 606-639. Mei, C. and Wentz, K. R. (1982). "Large-amplitude random response of angle-ply lami- nated composite plates," AIAA Journal, 20(10), 1450-1458. Mindlin, R. D. (1951). "Influence of rotary inertia and shear on flexural motions of isotro- pic elastic plates," Journal of Applied Mechanics, 18, 31-38. Murakami, H. (1986). "Laminated composite plate theory with improved in-plane responses," Journal of Applied Mechanics, 53, 661-666. Nanda, A. and Kuppusamy, T. (1991). "Three-dimensional elastic-plastic analysis of lami- nated composite plates," Composite Structures, 17, 213-225. Newland, D. E. (1993). An introduction to random vibrations, spectral and wavelet analy- sis, 3rd ed., Longman Scientific & Technical, New York. Newmark, N. M. (1959). "A method of computation for structural dynamics," Journal of Engineering Mechanics, ASCE, 85, 67-94. Nigam, N. C. (1983). Introduction to random vibration, The MIT Press, Cambridge. Noor, A. K. and Burton, W. S. (1989). "Assessment of shear deformation theories for mul- tilayered composite plates," Applied Mechanics Review, 42(1), 1-12. Parzen, E. (1962). Stochastic processes, Holden-Day, San Francisco. Petit, P. H. and Waddoups, M. E. (1969). "A method of predicting the nonlinear behavior of laminated composites," Journal of Composite Materials, 3, 2-19. Reissner, E. (1945). "The effects of transverse shear deformation on the bending of elastic plates," Journal of Applied Mechanics, 12, 69-77. 145 Reddy, J. N. (1984). "A simple higher-order theory for laminated composite plates," Jour- nal of Applied Mechanics, 51, 745-752. Reddy, J. N. (1987). "A Generalization of two-dimensional theories of laminated compos- ite plates," Communications in Applied Numerical Methods, 3, 173-180. Reddy, J. N. (1990). "A review of refined theories of laminated composite plates," The Shock and Vibration Digest, 22(7), 3-17. Roberts, J. B. and Spanos, P. D. (1990). Random vibration and statistical linearization, Wiley, New York. Sandhu, R. S. (1976). "Nonlinear behavior of unidirectional and angle ply laminates," Journal of Aircraft, 13(2), 104-111. Shinozuka, M. and Jan, C.-M. (1972). "Digital simulation of random processes and its applications," Journal of Sound and Vibration, 25(1), 111-128. Shokrieh, M. M. and Lessard, L. B. (1996). "Effects of material nonlinearity on the three- dimensional stress state of pin-loaded composite laminates," Journal of Composite Mate- rials, 30, 839-861. Singh, M. P., Khdeir, A. A., Maldonado, G. O. and Reddy, J. N. (1989). "Random response of antisymmetric angle-ply laminated plates," Structural Safety, 6, 115-127. Socha, L. and Soong, T. T. (1991). "Linearization in analysis of nonlinear stochastic sys- tems," Applied Mechanics Review, ASME, 44(10), 399-422. Spanos, P. D. (1980). "Formulation of stochastic linearization for symmetric or asymmet- ric MDOF nonlinear systems," Journal of Applied Mechanics, ASME, 47, 209-211. Spanos, P. D. and Iwan, W. D. (1978). "On the existence and uniqueness of solutions gen- erated by equivalent linearization," International Journal of Non-Linear Mechanics, 13, 71-78. Sun, C. T., Feng, W. H. and Koh, S. L. (1974). "A theory for physically nonlinear elastic fiber-reinforced composites," International Journal of Engineering Science, 12, 919-935. Sun, C. T. and Chen, J. L. (1989). "A simple flow rule for characterizing nonlinear behav- ior of fiber composites," Journal of Composite Materials, 23, 1009-1020. Vaziri, R., Olson, M. D. and Anderson, D. L. (1991). "A plasticity-based constitutive model for fibre-reinforced composite laminates," Journal of Composite Materials, 25, 512-535. Wang, S. S., Hu, H. T. and Hajali, R. (1995). "Effect of material nonlinearity on buckling and postbuckling of fiber composite laminated plates and cylindrical shells," Composite Structures, 33, 7-15. 146 Weaver, Jr. W. and Johnston, P. R. (1984). Finite elements for structural analysis, Pren- tice-Hall, Englewood Cliffs. Wen, Y. K. (1976). "Method for random vibration of hysteretic systems," Journal of Engi- neering Mechanics, ASCE, 102, 249-263. Whitney, J. M. and Pagano, N. J. (1970). "Shear deformation in heterogeneous anisotropic plates," Journal of Applied Mechanics, 37, 1031- 1036. Wirsching, P. H., Paez, T. L. and Ortiz, K. (1995). Random vibrations: theory and prac- tice, Wiley, New York. Witt, M. (1986). "Dynamic response and reliability of laminated plates to random load- ing," In Refined Dynamical Theories of Beams, Plates and Shells and Their Application, Elishkoff, I. and Irretier, H. (eds), Springer-Verlag, Berlin, 274-285. Witt, M. and Sobczyk, K. (1980). "Dynamic response of laminated plates to random load- ing," International Journal of Solids and Structures, 16, 231-238. Yang, P. C., Norris, C. H. and Stavsky, Y. (1966). "Elastic wave propagation in heteroge- neous plates," International Journal of Solids and Structures, 2, 665-684. Zienkiewicz, O. C. and Taylor, R. L. (1989). The finite element method, 4th ed., McGraw- Hill, New York. Zienkiewicz, O. C. and Cheung, Y. K. (1964). "The finite element method for analysis of elastic isotropic and orthotr0pic slabs," Proceedings of the Institute of Civil Engineers, 28, 471-488. Appendix A Higher-Order Moments of Two Gaussian Random Variables The nth moment of a random variable X, m" = E [Xn], is often difficult to evalu- ate directly using the definition E[X"] = H...jx"f(x)dx (A.1) depending on the form of the probability density function, f (x) . Fortunately an alternative method exists using the moment generating function. The moment generating function de- fined by m(s) = E[e‘x] = j esxf(x)dx (A.2) which is the Laplace transform of f (x) . It can shown that E[X"] = inmmIs = 0 (A3) ds It is relatively easy to evaluate higher-order moments if the moment generating function can be obtained analytically. For two random variables x and y, the moment generating function is defined as mm = E[e”“’l = j je‘“"f(x.y)dxdy (AA) —M —oo where f (x, y) is the joint probability density function of x and y. If x and y are zero-mean Gaussian variables, then 147 148 f(x, y) = [x2 - My ”21] (A5) 1 [— exp 2 Znoxoydl — p2 2010',“ - P ) 2 2 . . . . . where (SJr and G), are variances of x and y, respectively, and p 1S the correlation coefficrent between x and y (i.e. p = [y ). The joint moment generating function m(s, t) of x and x y y can be obtained analytically as 0000 m(s, t) = Elem’yl = j jesx+tyf(x.y)dxdy -99-” (A.6) l 2 2 2 2 _ exp[§(s ox + 2pstoxoy + t (5)] A general high-order moment may then be evaluated through a’" + "m(s, t) E[xmy"l = m n as at (A.7) s,t=0 If m + n is odd, then E [xm y"] = 0. Expressions are given here for a few cases in which m + n is even. 4 34 2 2 (1) E1712] = —, = 3Elvul as s,t=0 1 5 86m 2 2 (2) E[(le)(723)l = l 5 = 15517127315173] as a: s 1:0 5 5 8mm (3) E[(712)(723)1= 5 5 as at 3 i=0 2 2 2 = 225EI’Yf2] £1723] El:‘YIZ'Y231l +600E[712] E[stl 51712723] 0 o +120E['y?2] 51723] 13171272315 149 The various covariances in Section 4.4 can be evaluated using the technique outlined here. More detailed descriptions on higher-order moments of Gaussian random variables are giv- en by Parzen (1962) and Bendat and Piersol (1986). Appendix B Transformation of Incremental Elastic Constants The incremental elastic constants in material coordinates is defined in Eq. 5.6 and 5.10 by [CT] = lQl+lQNl (B-l) where [Q] is given by Eq. 3.5 and [QN] = [1112140, 0, 2 a6,(2i+1)yf;, 2 a4i(2i+1)y§;, 2 a5i(2i+ 1”ng (13.2) i=1 i=1 i=1 The incremental matrix of elastic constants in global coordinates, [CT] . was formally in Eq. 5.21 as [ET] = [Q] +{IT11{2}[T31[Z a,.-(2i)([T31{e})2“ ‘] + Hub: a,.(tT31{e})2‘]} i-l i=1 _ +{[T;]{8}1T41[ n a4,(2i)([T4]{8})2i’ ‘] + [7316 a4.-([T41{e})2’]} (3'3) i=1 '=1 +{[T§1{e}trsi[i a5,(2i)([T5]{e})2i’ ‘]+ [131E a,,([15]{e})2‘]} i=1 i=1 It is proved here that this matrix may also be obtained through the transformation [61] = [Tl-IICTHTI’T (13.4) Focusing on the in-plane nonlinear terms (second term on the RHS of Eq. B3), 150 151 [T1"[diag<0.0,1,0,0)][TJ‘TIITJ’T1"{sum [TIHeHm = [T]“[diag(0. 0, 1. 0. omnmrm = [Tl"l[diag(0,0,1.0.0)]{81.82.7121Y231713}T1T31 (3.5) = [T1’1{0.0,712,0,0}T[T31 = ”11712 Using Eq. B5 and noting that [T3]{8} = 712 , the second term on the RHS of Eq. B.3 be- comes [TI]{e}[T31[2 a,.-(2i)(tT,1{e})2"“]+ [1}][2 a,,([T31{e})2‘] i=1 i=1 = ”11(2) 06.121“ + will] (B.6) i=1 “6i(Y12)2i] '- 1 (— 2 1,4207%;- ‘] + [TIJ[ l=1 [ritual Using Eqs. B], B2 and 3.4, the above equation becomes 152 1511 = 1T1“1Q1111’T+1T1"1QN11T1‘T = [Q] + [1]“ [diag(O, 0, 1,0, onm‘TX (16,121+ mfg] i=1 + [T]_l[diag(0, 0, 0, 1, 0)][T1'T[X 04.12” 1W3; i=1 n (B.7) + [T]_1[diag(0, 0,0, 0, 1)][T]‘T( 2 115,121+ mfg] i=1 = [a] + [T:][2 061(25+1)Yfi] i=1 + [T;][ Z 114,121+ 1)y§§]+ [T;](2 115,121+ 1)yf§] i=1 i=1 The second term on the RHS of Eq. B.7 is identical to the RHS of Eq. B.6 which is the sec- ond term on the RHS of Eq. 33. Similarly, the third and fourth terms on the RHS of Eq. B.7 are equal to the third and fourth terms on the RHS of Eq. B.3. Thus Eqs. B3 and 3.4 are equivalent. The stresses in global and material coordinates are related through {0'} = {THO} _1 (B8) {6} = [T] {6'} and the strains in global and material coordinates are related through . _ —1 = —T {8} — [RllTllRl {8} [T] {8} (8.9) {a} = 1171"]"1e'1 Where [R] = [diag( 1, l, 2, 2, 2)] is the Reuter’s matrix and 153 ' T {0} = {0110221122123’113} ' T {e} = {€1,821Y121'Yz3’713} (B.10) T {a} = {0x9 6)” IX)” Tyz’ sz} r {e} = {epeyflxyflyzflxz} [T] is the rotational transformation matrix in Eq. 3.12 and [Ti] is the ith row matrix of LT; = [—2cosesin0 2cosGsin0 cosze— sin20 0 0] "I A I ' [0 0 0 cos0 —sin0] (B-ll) "I u. I ‘ [0 0 0 sin0 —cosO]