- Email: [email protected]
Development of an anisotropic, multilayered, shear-deformable rectangular plate element
Computers & Sfrucfures Printed in Great Britain.
Vol. 21. No. 4. pp. 789-l%.
0045-7949185 $3.00 t .oo 0 1985 Pergamon Press Ltd.
1985
DEVELOPMENT OF AN ANISOTROPIC, SHEAR-DEFORMABLE RECTANGULAR Department
MULTILAYERED, PLATE ELEMENT
MARCO Dr SCIUVA of Aerospace Engineering, Polytechnic of Turin, Italy (Received
23 July 1984)
Abstract-A
multilayered anisotropic flat plate element which includes the effects of the transverse shear deformation is developed by making use of the displacement formulation. The discrete element is a rectangle with 32 degrees of freedom which include extension, bending, and transverse shear deformation states. The formulation is based on an improved bi-dimensional transverse shear deformation plate theory which accounts for piecewise linear distribution across the thickness of the inplane displacements u and U, and allows the contact conditions at the interfaces between the layers to be satisfied. In order to demonstrate the accuracy and efficiency of the developed finite element, the solutions for two sample problems are given. Problems, for which analytical solutions are available, were selected so that a comparison between the analytical and numerical solutions could be made. The first example is the well-established cylindrical bending problem of cross-ply and angle-ply laminates; the second one is the bending of a symmetric cross-ply square plate simply-supported on all edges. The tests carried out show that the element is very efficient in predicting the responses of thick
and thin laminated plates. Also predicting, at very low span-to-thickness ratios and models, the warpage of the cross-section, which has been found to be unfeasible by other conventional finite-element displacement
approaches.
NOMENCLATURE
Cartesian coordinates with z-axis oriented in thickness-wise direction and measured from the top surface natural coordinates for the rectangular finite element in-plane dimensions of rectangular plate laminate thickness in-plane dimensions of rectangular finite element total number of layers in the laminate span-to-thickness ratio displacements in x-, y-, and z-direction, respectively values of the displacements on the reference surface values of the shear rotations in the (x, z) and (y. z) planes, respectively, ofline normal to the reference surface values of the total rotations in the (x, z) and (y, z) planes, respectively, of line normal to the reference surface values of the z-coordinate at the interface between the kth and (k + 1)th layers functions determined by satisfying the contact conditions on the transverse shearing stresses at interface k constants of continuity for the transverse shearing stresses heaviside unit function strain energy function stress components in Cartesian coordinates strain components in Cartesian coordinates elements of stiffness matrix for individual layer connecting stresses and strains elements of reduced stiffness matrix for 789 CAS 21:4-L
individual layer defining constitutive relations of plane stress type EL. ET elastic moduli of individual layer in lon-
GLT,
GTT
VLT. VTT L,
T
rY
lo-/ (.).*
gitudinal and transverse directions, respectively shear moduli of individual layer Poisson’s ratios longitudinal and transverse directions of unidirectional layer, respectively matrix transpose of matrix column matrix row matrix null matrix partial derivative with respect to the Xcoordinate
INTRODUCTION
Exact three-dimensional analysesll-31 have indicated that, for laminated anisotropic plates, the distortion of the deformed normal due to transverse shear is dependent, not only on the laminate thickness, but also on the orientation and degree of orthotropy of the individual layers. From the above, it is evident that any approximate effective bidimensional theory for multilayered anisotropic plates should account for the effects of transverse shear deformation. In this respect, the axiomatic approach has been employed either by assuming the distributions of the transverse shear stress uXxzand u,,[4-61, or by assuming a displacement field[7141. In Ambartsumyan’s shear deformation theory[4] for specially orthotropic, symmetric laminates, the distributions of the transverse shear stresses crxrand u,.~ in each layer are assumed to be known. Whit-
790
M. DI SCIUVA
hand, the finite element method provides a convenient and powerful technique for analyzing complex structures. There has been considerable effort of late directed toward the development of shear deformable. multilayered anisotropic plate finite elements[lS-191. As a basis for the development of a displacement plate element of this type, the previous authors have employed the shear deformation plate theory developed by Yang et a1.[7]. This approach requires an arbitrary shear correction factor which is dependent on the lamination parameters for obtaining Improvements of this formulation have been accurate results: furthermore. it does not model the warpage of the cross-section, which is very strong suggested by Srinivas[9] by assuming the previous ratio values. displacement field for each layer. High-order the- for low span-to-thickness Therefore, it seems desirable to develop a finite ories have been developed by Lo et al.[ lo] and Murelement displacement approach based on a refined thy[ll]. shear deformation plate bending theory which alFor plates consisting of layers perfectly bonded lows one to model, among other deformation states, together (the layers function concurrently without slippage), we know from elasticity theory that, at the warpage of the cross-section. The purpose of this article is to develop. on the the interface of the kth and (k + 1)th layers, the basis of the multilayered anisotropic plate theory following conditions must hold? proposed in [13], a rectangular flat finite element which includes the distortion of the deformed nor= bI’k+, M’k l’k = T’k+ I l(k = llk+l mal to the reference surface, the transverse shear (2) effects, as well as the coupling effects of stretching and bending. u;:,k = u;;.k+ I. uyz.k = u_vz.k+I Uxz.k = ux.-.k+ I Some preliminary numerical results show that (3) the developed finite element is very efficient in predicting gross behavior of thick and thin multilayered From eqn (2) it follows immediately that anisotropic plates.
ney[51 has extended this theory to symmetric laminates with arbitrary orientations. Refinements of the Ambartsumyan’s theory, including all transverse effects in a consistent manner, are presented by Voyiadjis and Baluch[6]. In the shear deformation theory developed by Yang et a/.[73 and applied to the specific problem of bending and flexural vibrations by Whitney and Pagano[8], the following displacement field is assumed for the whole plate behavior
PLATE
Continuity conditions cannot be imposed on the remaining stresses and strains; in general, these are discontinuous at the interfaces. Since the transverse shear strains derived by the displacement field given by eqn (1) are constant across the thickness, it is not difficult to realize that the contact conditions on the stresses cannot be satisfied if the layers have different mechanical properties. In a series of papers[l2-141, the writer has proposed a linear displacement model which allows the contact conditions to be satisfied simultaneously. Such a displacement model has been utilized to analyze the flexural behavior of the multilayered orthotropic[ 121 and anisotropic[ 13, 141 plates. Numerical tests carried out on the cylindrical bending of an infinite symmetrically laminated strip have demonstrated the accuracy and efficiency of the proposed model. Because of the complexity of the theories involved, the investigations discussed above are severely
limited
in their
applications.
Displacement
EQUATIONS
field
The laminate under consideration has constant thickness h and is composed of a finite number N of thin orthotropic layers of uniform thickness perfectly bonded together. The principal axes of elasticity of any individual layer are assumed to be oriented arbitrarily with respect to the laminate axes, and the reference surface is superposed on the z = 0 plane (see Fig. I). The material properties and the thickness of each layer may be entirely different. The constitutive relations for any individual layer are given byS
On the other
where the C, are the elastic coefficients t In the following the index k (k = 1, N - 1) will refer to interfaces and the index s (S = 1, N) to layers. N being the number of layers.
and the
$ It must be noted that, in the given constitutive relations, only nine of the elastic coefficients are independent.
791
Development of an anisotropic, multilayered, shear-deformable rectangular plate element
faces (for the undefined symbols, see Nomenclature) . In examining the relations (7)-(9), it is not difficult to realize that the displacements u and v are continuous functions of the z-coordinate for all values of the (pp(x, y) and x&r, y) functions. It follows that the expressions for (Pk and Xk may be found from the contact conditions for the transverse shearing stresses. Once these conditions have been satisfied, it is found that eqns (7)-(9) can be rewritten as1131
Fig. 1. Coordinate system and numbering of the layers and their interfaces.
x ?J
=
u”
+
dy,
-
@v,
+
(z x
k
usual notation for stresses and engineering strain components has been adopted. As usual, the normal stress oZZis assumed to be small in comparison with other normal stresses and is neglected. Taking into account this assumption, the constitutive relations assume the following contracted form
x w
=
(Z
-
zk)Yb
(dkrk
-
-
+
zk)
(7a)
Zk)
(84
bky,)
Zk)Y(Z
-
HP
@a)
where ak, bk, ck, and dk are known constants only depending on the mechanical properties of various layers. Strain
energy
The strain energy per unit area of the plate is given by or
@ = g Q(s) s=l
(10)
where and
@‘(s)= ;
j-:‘,{dr
(11)
{E}dz
(6) stands for the strain energy per unit area stored in the sth layer and
or Iat) = [Q2(s)l k,I
64
where Qij = Cij for i, j = 4, 5 and Qij = C, Ci3Cj3IC33 for i, j = 1, 2, 6. The definition of the matrices introduced in the relations (5) and (6a) follows directly from the relations (5) and (6). The displacement field is assumed to be of the following form[ 131
By making use of eqn (11) in conjunction with the eqns (7)-(9) and the constitutive relations, it is not difficult to show that eqn (10) assumes the following matrix form @ = 1 {6}r[E]{6}
(lOa)
where U = U0 + z(rx - $!r) + z (Pk(Z - Zk)Y(Z - Zk) k
W (7)
7.J= Y0 + z(y, - $,)
+ 2
Xk(Z
-
zk)y(z
-
= LupxG!X& $V b % YXJ Y.v.xYX.?Y?..” -
0 WJX
-
0 w,yy
-
w!ryJ
.
Zk)
Explicit expressions for the elements matrix are given in the Appendix.
k
of the [El-
(8) w
=
w”.
(9)
Here Y(z - zk) is the Heaviside unit function and the summation is extended over the N - 1 inter-
FINITE ELEMENT FORMULATION
As is known, the accuracy of the finite element solution depends on the ability of the assumed in-
M. DI SCILJVA
792
terpolation functions to accurately model the deformation modes of the structure. Since the coupling between extension, bending, and transverse shear deformation states affects the static and dynamic responses of the laminated anisotropic plates, a correct finite element displacement formulation should model such a coupling. The results obtained by Prior and Barker[lS] show that a flat rectangular plate element with seven degrees-of-freedom (DOF) at each nodal point should be accurate for this purpose. The nodal parameters consist of the three generalized displacements u’, Y’, and N”: the two shear rotations yr and y,., and the two total rotations WY)= 4,Yand wpY = 6,. This displacement model fully satisfies the inter-element compatibility conditions. However, it has a great disavantage in that it does not include the constant twist term. In order to improve the formulation, we have added the mixed second derivative uzp,,= 6,X,; thus, the proposed rectangular plate element has eight DOF at each nodal point, u”, zf”, w”, 6,, 6?, 6,, yI. yy , giving a total of 32 DOF. Linear polynomials have been used for approximating the inplane displacements [lo, zl”, and the shear rotations yX and Y?; cubic polynomials have been used for approximating the transverse displacement MI’. Let us call the nodes as in Fig. 2. It can be shown that the following expressions satisfy the inter-element compatibility and the constant strain conditions[21] (u stands for 11’. zf”, yX, y,,):
and similarly for a’)‘(q), H!(q), etc., replacing 5 with q and a, with b, in above. The vector (6)’ comparing into eqn f 10al can be written in terms of the nodal degrees of freedom by taking the appropriate derivatives of the interpolation polynomials [eqns ( 1211. When this is given, we obtain
where {cl}, is the vector of the degrees of freedom of the nodes surrounding element e, and [N] is a rectangular matrix (13 x 32) whose elements are products of the polynomials H by itself or by their derivatives. The strain energy of the element can be determined by substituting expression (13) into eqn (lOa) and integrating over the area of the element. The result is
where
is the rectangular plate element stiffness matrix (32 x 32). External loadings and boundary conditions are subject to the usual finite element analysis considerations. NUMERICAL EXAMPLES
In order to assess the accuracy and efficiency of the developed rectangular finite element, some preliminary numerical tests have been carried out. This evaluation was conducted on two sample problems, for which analytical solutions are available. In this work, for a definition of analytical solutions we where (5 = xla, and q = y/b,) mean the values obtained by solving exactly the system of the differential plate equations based on Hm = 1 - 5; H6Ykc) = 5; HIaS) = 1 - 35’ the displacement field given by eqns (7a)-(9a). + 26’; H6:‘(5) = S2(3 - 25) The first example is the well-established cylindrical bending problem of cross-ply and angle-ply fi:‘(Sl = a,W - 2[ + 11: @t?(Q = a,S*(E - 11 laminates; the second one is the bending of a symmetric, cross-ply square plate simply-supported on (2,2) all edges. (1,2) In addition to comparing the approximate solution with the analytical one. the convergence of the discretization error was also examined numerically . In the following, we consider a three-layered, symmetric laminate with layer properties a,
y~q~z be[ (1,lI
I
5
(2,l)
w
EJEr G&ET
= 25, = 0.2,
GLrlET
= 0.5,
VLT = VT7 = 0.25.
X
Fig. 2. Global and local coordinate systems and node numbering for rectangular olate element.
where L and Tare the directions parallel and normal to the fibers, respectively, and vLTis Poisson’s ratio
Development
measuring transverse parallel to the fibers.
of an anisotropic,
multilayered,
strain under uniaxial
shear-deformable
CROSS-PLY
stress
r0’/90”/0’,:
AnalytIcal 0
2.0
Cylindrical bending
Flnlte
Sotu tlon
Element
~WWF.E.S.
-
WA.S.YWA.S.
where Wr.n.s. stands for finite element solution and WA.s, is the analytical solution. It seems that in both laminate configurations the convergence is very rapid. Figure 4 gives the analytical and finite element values for W as a function of the plate thickness ratio
Cl21
Solution
/
1.0
SolutionC13.141
X Finite
Element
//
I’
ANGLE-PLYt30’/-30’/30): --Anatytlcal
The plates considered in this example are of infinite length in the y-direction, uniformly supported along the edges x = 0, a and subjected to a sinusoidal distributed transverse loading p;(x) = p. sin Ax (A = r/a) on the reference surface. The cross-ply laminate (0”/90”/0”) has three layers of equal thickness with the fibers in the outer layers oriented in the x-direction (0”) and those in the inner layer oriented in the y-direction (90”). The angle-ply laminate (30”/- 30”/30”) has ply thickness (h/4, h/2, H/4), respectively, with fibers in the outer layers oriented at +30” and those in the inner layer at - 30”, with respect to the x-axis. Under these conditions, the deflected surface is cylindrical; therefore, only a strip of the plate needs to be modeled for the computational purpose. Furthermore, due to symmetry, one only has to model a half span. Figure 3 shows the convergence patterns of the nondimensional transverse displacement, W, at the ratio) center of the span for p (span-to-thickness equal to 10. In this figure, as well as in Fig. 5, we have posed %E =
-
2.2
193
rectangular plate element
/
Solution
I’
/’ /
1.4 1.2
t t
I
I
6
I
8
I
10
I
12
I
14
I
16
I
16
I
2
Fig. 4. Typical results for the center deflection of a threelayered, symmetric, simply-supported laminate in cylindrical bending.
For both laminate configurations, and for all values of the thickness-to-span ratio, the approximate solution is very close to the analytical one given in Refs. [12-141. Close agreement between numerical results and analytical solutions is also obtained in the distribution of the in-plane displacements and flexural stresses across the plate thickness. Owing to this agreement, and also observing the fact that the form of the various distribution curves correspond closely to those given in [ 12-141, these curves are not shown here. Simply-supported
square plate
The numerical results given in this section refer to a three-layered, symmetric cross-ply, square
h/a.
The plotted finite element values were obtained with a 5 x 1 mesh which modeled the behavior of a half span.
1
IE
-10 $I
-20 M
CROSS-PLY
(0390”/0”,
X--3(
ANGLE-PLY(30”/-JO”f3@“. -30
-10
-50
N”lnb.x
Fig. 3. Cylindrical bending of an infinite three-layered, symmetric, simply-supported laminate under sinusoidal transverse loading POsin Ax. Maximum deflection error as a function of the number of elements.
of etementr. NE
Fig. 5. Bending of a three-layered, symmetric cross-ply (0”/90”/0”), simply-supported square plate under sinusoidal transverse loading PO sin Ax sin Xy (A = n/a). Maximum deflection error as a function of the number of elements.
M. Dt
f%XJVA
P
No”dimenlio”al
in-plane
displacement.
0 = ET”Lo,a/2,r,
/hp,
h
._ ; I B ;
Fig. 8. Variation of u through the thickness of a threelayered, symmetric cross-ply (0”/90”/0”), simply-supported square plate under sinusoidal transverse loading po
-.2 2
4
6
0
10
12
14
16
10
I
Fig. 6. Maximum deflection of a three-layered, symmetric cross-ply (O”/W/o”), simply-supported. square plate under sinusoidal transverse loading p. sin Ax sin hy (A = P/U).
plate simply-supported on all edges and subjected to a sinusoidal distributed transverse loading &(x, y) = p sin (TX/U) sin (phyla) on the reference surface. The laminate has three layers of equal thickness, with the fibers in the outer layers oriented in the xdirection (0”) and those in the inner layer oriented in the y-direction (90”). The analytical results are obtained by choosing u”, u’, w”, y,, and yv as u’ =A,cosAxsinhy: u” = A,,sinhrcoshy; M’O= A,, sin Ax sin hy ; (A = 57/u) ?/x = A; cos Ax sin Ay ; yv = A; sin Ax cos Ay
which satisfy the differential equilibrium equations and the simple support, without in-plane normal restraint boundary conditions given in [12]. Figure 5 shows the convergence pattern of the nondimensional deflection ti; at the center of the
-10
-12
CROSS-PLY ---Elasticity -.‘A
Thin Flnile
-4
plate for p = 10. Because of the symmetry, only one quadrant of the plate was modeled and NE stands for the number of the elements for a quarter of the plate. Fast convergence to the analytical solution is observed: for the 5 x 5 mesh the error in the approximate results is not higher than 2%. Thus, this mesh has been used to obtain the results given in the following. Figure 6 compares the exact, analytical, and finite element values for M:as a function of the plateto-thickness ratio h/u. The distribution of the flexural stress a,, as resulting from the present finite element approach is given in Fig. 7 in conjunction with those resulting from the elasticity solution[2], the approximate bidimensional theory developed in [ 121, and the classical lamination theory (C.L.T., Kirchhoff thin plate theory)[201. The distributions of the in-plane displacement U are plotted in Fig. 8. Close agreement of the analytical and finite elements results, with the exact elasticity solution]21 being obtained. As expected, at very low span-tothickness ratio values the thin plate theory is in-
-2
(0~/9O~/V): Solu~lon
Adytic.l
-6
-13
sin !LXsin Ay CA= T/U).
;
/I/
sd”tlon
Plate Theory Element
(C.L.T.)
Solution
INTERFACES
Nondimensional
bending
stress,
zXX
= bXX(a/2,a12,z)
Ip,
Fig. 7. Variation of c,, through the thickness of a three-layered, symmetric cross-ply (0”/90”/0”) simplysupported square plate under sinusoidal transverse loading po sin AXsin hy (A = n/n).
79.5
Development of an anisotropic, multilayered, shear-deformable rectangular plate element adequate to model the correct behavior of the multilayered composite plates.
deformationtheoryformultilayeredorthotropicplates. Proc. of the VII Congress0 Naz. AIDAA, Oct. 1983, Naples, Italy. In press in 1’Aerotecnica Missiii e Spazio.
13. M. Di Sciuva, A refinement of the transverse shear deformation theory for multilayered anisotropic
CONCLUDING REMARKS
Based upon a refined approximate shear deformation theory proposed by the author in a series of previous papers[ 12-141, an effective multilayered anisotropic rectangular plate element is developed. The finite element is a rectangle with 32 degrees of freedom which include stretching, bending, and transverse shear deformation states. Numerical tests carried out on two sample problems show that the element is very efficient in predicting responses of thick and thin laminated plates under static loading. also at very low span-to-thickness ratios, and models the warpage of the crosssection, which has not been found to be feasible by other conventional finite element displacement approaches[lS191. It is believed that the approach can conveniently be extended to perform static, vibration, and buckling analyses for a large class of laminates with varied loading and boundary conditions. This is in the planning for a future article. Acknowledgment-This
research has been worked out by the financial support of the Italian Ministry of Education (MPI).
plates. Department
(1981).
18, J. N. Reddy and W. C. Chao. Nonlinear oscillations
Mater. 4, 20 (1970).
3. N. J. Pagano, Influence of shear coupling in cylindrical bending of anisotropic laminates. J. Compos.
Vibration of thick rectancrular mates of bimodulus composites materials. J. Aipl. iech. 48, 371 (1981). 20 R. M. Jones, Mechanics of Composites Materials. McGraw-Hill (1975). 21 F. K. Bogner, R. L. Fox and L. A. Schmit, The generation of inter-element-compatible stiffness and mass matrices by the use of interpolation formulas. AFFDL-TR-66-80, Proc. of Conference on Matrix Methods in Structural Mechanics, p. 397 (1965).
Mater. 3, 534 (1969).
6. G. Z. Voyiadjis and M. H. Baluch, “Refinements in the Bending of Plates with one Plane of Elastic Svmmetry”. P&c. of the Intern. Symposium on the Mechanical Behaviour of Structured Media, 507. May 1981, Ottawa, Canada. P. C. Yang, C. H. Norris and Y. Stavsky, Elastic wave propagation in heterogeneous plates. Int. J. Solids and Struct. 2, 665 (1966).
J. M. Whitney and N. J. Pagano, Shear deformation in heterogeneous anisotropic plates. J. Appi. Mech. 37(4), 1031 (1970). S. Srinivas. A refined analysis of composite laminates. J. Sound and Vib. 30(4), 495 (1973).
10. K. H. Lo, R. M. Christensen and E. M. Wu, A highorder theory of plate deformation. Part 2: Laminated plates. Trans. ASME, Ser. E, J. Appl. Mech. 44, 669 (1977). 11. M. V. V. Murthy, An improved transverse shear deformation theory for laminated anisotropic plates. of the transverse
By following the same procedure outlined in [12], it is readily shown that eqn (11) takes the compact form
where
shear
101 lEl3O)l l-&(S)1 [Ezz(s)l WI 101
]Ellts)l
[E(s)1 =
Plates.
Translated from the Russian by T. Cheron and edited by J. E. Ashton, Technomic Publishing Co. (1969). 5. J. M. Whitney, The effect of transverse shear deformation on the bending of laminated plates. J. Compos.
NASA TP I903 (1981). 12. M. Di Sciuva, A refinement
rectangular plates. J. Appl.
of laminated, anisotropic,
Mech. 48, 396 (1982). I,~ *,, C. W. Bert, J. N. Reddy, W. C. Chao and V. S. Reaay, IY.
APPENDIX
1. N. J. Pagano, Exact solutions for composite laminates in cylindrical bending. J. Compos. Mater. 3, 398 (1969). 2. N. J. Pagano, Exact solutions for rectangular bidirectional composites and sandwich plates. J. Compos.
Theory of Anisotropic
Engineering-Poly-
M. Di Sciuva, A refined transverse shear deformation theory for multilayered anisotropic plates. In press in Atti Accad. delle Scienze di Torino (Italy). IS. C. W. Pryor and R. M. Barker, A finite-element analysis including transverse shear effects for applications to laminated plates. AIAA J. 9(5). 912 (1971). 16. A. K. Noor and M. D. Mathers, Shear-flexible finiteelement models of laminated composite plates and shells. NASA TND-8044 (19751. 17. .I. N. Reddv and W. C. Chao. Laree deflection and large ampliiude free vibrations’ of laminated composite-material plates. Comput. Structures 13(2), 341
REFERENCES
Mater. 4, 330 (1970). 4. S. A. Ambartsumyan,
of Aerospace
,4. technic of Turin (Italv)-ReDort n.5 (1983).
[E33(s)l
sYk.?&C
[E34(S)l LE44Cs)l I
with
lE130)l = z[Ei~(sIl
[EIICSII= [AITIQ,WllAl: + lAITIQi(s)llB*l [EM(
= zlAITIQ,tsIl;
[E340)1
= z([E14(~)1 x [,544(s)]
LE33b)I
= z*[EIIWI
+
[W[Q,b)i
[&z(s)1 = lQz(IIlt1Il + [A*])
+ [B*lTIQ~(sN); = zz[Q,(s)l +
zt[Al*lQ,(s)l IB*l
w*il? + [B*i%21(~~i
LB*].
The undefined matrix symbols which appear in the previous relations stand for rl000i
[A]=
0001 ; I ni in I
TT T7
[A*]=
;; 1
J
M. Dr SCIUVA
796
+ 2D& E(10, 10) = 022 + 2042 + 2Dk, + D 2626 + DG.; E(9, 12) = Dxj + Dh + D$; E(l1. II) = D,,
where
--__
(4, B. c, 0) _ _ (A, B, C, D) 1
a-1
E(9. 13) = Da + D& + D$,; E(l0. 11) = D,z + D$ + Dfe; E(10,12) = Dzz + D$ + Dk:
Substitution of the relation (1A) into eqn (IO) gives the following expressions for the elements of the symmetric matrix [El [see eqn (lOa)]: E(1. 1) = A,,: E(1,2) = E(l, 3) = A,Q(l, = Aiz; E(1. 7) = B,l + By, + B,h,
4)
E(1, 8) = 816 + Bib + BY,; E(1. 9) = 816+ BQ6 + BP2; E(l, 10) = B,? + B?z + BfCj: E(1, 11) = B,, E(1, 12) = B,z; E(1, 13) = 816: E(2. 2) = E(2, 3) = AM; E(2, 4) = Ax:
E(2,7) = B,c, + Bik + B&6; E(2, 8) = Be6 + B$ + IS;
E(2, 9) = Bfjb + I?&, + B$:
E(10,13) = 026 + D$ + D&i; E(l1, 12) = D,*; E(11, 13) = Dlh: E(12, 12) = D22; E(12, 13) = Dxj: where [Ail, B,j, Dijl = < 11, Z, ~‘1 Qtj’xs) > for i, j = 1, 2. 6 and A;j= hQij(l) for i, j = 4, 5 &I&2&6&%&6
E(4, 9) = BZb + B% + B$; E(4, 10) = B22 + B$* + B&: E(4. 11) = B,z: E(4, 12) = &*; E(4. 13) = Bx,: E(5,5) = AM + A& + A&; E(6,6) = Ass + A& + A&; E(5. 6) = Aa + At + Ai4 = A45 + A% + A&: E(7. 7) = D,, + 2DQ, + 2D$ + DQf
+ 0:: + 207:;
=
<[I:
B%2Bf2&‘6B:6@6
~Q,,~~~Q,~~~~Q,~~~~Q~6oa66(s)1AA ’
1
=
[ D$2Dt2DtbD46D$
<:
[Q22(~)Q,z(~~Q,~(~)Q26(s)e6603 Bk ’
[1
&,BT2&6&sB& [ D~IDTJWX.S~&
<[l;
E(3, 9) = 866+ B&i+ B&k: E(3, 10) = 826 + B% + B&6; E(3, 11) = B,c; E(3. 12) = 826 E(3, 13) = 866; E(4, 4) = Azz; E(4, 7) = B,2 + By* + B$; E(4, 8) = BZ6 + B& + BC*:
1
LDfl,D?zDf16Dy6D&
E(2. 10) = 826 + B$ + B&; E(2, 11) = 816; E(2, 12) = 826; E(2, 13) = I&,$,; E(3, 3) = A6,5; E(3, 4) = Ax,; E(3, 7) = B,6 + B;16 + B&6; E(3, 8) = B@ + B$ + Bie;
E(13. 13) = DM,
1
=
[Q,,(s)Q,~(s)Q,s(s)Q26(S)e6601
c/, ’
I DA > <[If [Qzz~~)~I~~~)~
&2Bi’2Bi’6&!6&7
=
L&2Df2Df6’6D$D$
[D8DiWP,“D~l
= < [Q,,(s)Q,~(s)Q~~(s)Q~~(~)I &‘b
’
U%D~:D~iX%l
= < [Qzz(s)Q,s(s)Q~~(s)Q~~(~)~ &B, ’
[DYiD(‘6D%D%I = < [QI ~(S)Q,~(S)Q~~(S)Q,(S)~ ckc, > [D%Di’$!D#D&‘] = < [Qz~(S)Q,~(S)Q~~(S)Q~~O~DkD, > [DfiDYtDtiD@l
= < [Q,z(s)Q,s(s)Q26(s)Q6601 &Br ’
E(7,8)
= D,6 + Dyh + D& + Di', + D& + Dyi + Dd + D$ + Den b“; E(7, 13) = D,6 + D& + 0th:
[D3Di’kD%D~‘l = < [Q,,(s)Q,b(s)Q26(s)Q66(s)l A&r >
E(7.9)
= D,c, + 2D& + D$ + D& + Dy: + D(;f + D&' + D@; E(8. 8) = 0% + 2D$ + 2D&, + 0';; + 2Df; + Ott,";
[Db’D%D&D&] = < [Q,2(s)Q,b(s)Q26(s)Q6601 B&r >
E(7. 10) = D,2 + D$ + D& + Dye + 0% + Dr$' + D?: + 046"+ D&; E(7. II) = D,, + Dy, + D$,;
[DdDi%%@d]
= < [Q,z(s)Q,h(s)Q26(s)Qb6(s)l A/A
[D%D%!D&!D&‘]= < [Q~z(s)Q,~(S)Q~~(S)Q~~(S)~ BkDr > [DFfDf6do’i66D& = < [Q,z(s)Q,~(s)Q~~(S)Q~~(~)~ CkDr >
E(7, 12) = D,? + D& + D&6: E(8. 9) = Dhb + D&
+ D&, + Die + D$*,+D Ebb + DQ; + 0% + D,z. E(8, 10) = Dzb + 2D& + Dyz + D&, + D$'k+ Df$ + Dig + D&; E(9. 9) = D66 + 2046 + 2D& + Dg: + D‘$ + 2D$ E(8. 11) = D,6 + Dt6 + DC,; E(8. 12) = Ox, + D$ + Diz; E(8, 13) = Db6 + D& + DL,; E(9, 10) = Dzb + D26 + D$ + D;;, + D$ + 04; + D$ + D$$ + Dg;; E(9, 11) = D,6 + Dye + D:,;
>
For notational convenience,
we have posed
(..) = \$,(..)dr:
G, = x&(z
where t-l =OL= =:;
-z,),
Vol. 21. No. 4. pp. 789-l%.
0045-7949185 $3.00 t .oo 0 1985 Pergamon Press Ltd.
1985
DEVELOPMENT OF AN ANISOTROPIC, SHEAR-DEFORMABLE RECTANGULAR Department
MULTILAYERED, PLATE ELEMENT
MARCO Dr SCIUVA of Aerospace Engineering, Polytechnic of Turin, Italy (Received
23 July 1984)
Abstract-A
multilayered anisotropic flat plate element which includes the effects of the transverse shear deformation is developed by making use of the displacement formulation. The discrete element is a rectangle with 32 degrees of freedom which include extension, bending, and transverse shear deformation states. The formulation is based on an improved bi-dimensional transverse shear deformation plate theory which accounts for piecewise linear distribution across the thickness of the inplane displacements u and U, and allows the contact conditions at the interfaces between the layers to be satisfied. In order to demonstrate the accuracy and efficiency of the developed finite element, the solutions for two sample problems are given. Problems, for which analytical solutions are available, were selected so that a comparison between the analytical and numerical solutions could be made. The first example is the well-established cylindrical bending problem of cross-ply and angle-ply laminates; the second one is the bending of a symmetric cross-ply square plate simply-supported on all edges. The tests carried out show that the element is very efficient in predicting the responses of thick
and thin laminated plates. Also predicting, at very low span-to-thickness ratios and models, the warpage of the cross-section, which has been found to be unfeasible by other conventional finite-element displacement
approaches.
NOMENCLATURE
Cartesian coordinates with z-axis oriented in thickness-wise direction and measured from the top surface natural coordinates for the rectangular finite element in-plane dimensions of rectangular plate laminate thickness in-plane dimensions of rectangular finite element total number of layers in the laminate span-to-thickness ratio displacements in x-, y-, and z-direction, respectively values of the displacements on the reference surface values of the shear rotations in the (x, z) and (y. z) planes, respectively, ofline normal to the reference surface values of the total rotations in the (x, z) and (y, z) planes, respectively, of line normal to the reference surface values of the z-coordinate at the interface between the kth and (k + 1)th layers functions determined by satisfying the contact conditions on the transverse shearing stresses at interface k constants of continuity for the transverse shearing stresses heaviside unit function strain energy function stress components in Cartesian coordinates strain components in Cartesian coordinates elements of stiffness matrix for individual layer connecting stresses and strains elements of reduced stiffness matrix for 789 CAS 21:4-L
individual layer defining constitutive relations of plane stress type EL. ET elastic moduli of individual layer in lon-
GLT,
GTT
VLT. VTT L,
T
rY
lo-/ (.).*
gitudinal and transverse directions, respectively shear moduli of individual layer Poisson’s ratios longitudinal and transverse directions of unidirectional layer, respectively matrix transpose of matrix column matrix row matrix null matrix partial derivative with respect to the Xcoordinate
INTRODUCTION
Exact three-dimensional analysesll-31 have indicated that, for laminated anisotropic plates, the distortion of the deformed normal due to transverse shear is dependent, not only on the laminate thickness, but also on the orientation and degree of orthotropy of the individual layers. From the above, it is evident that any approximate effective bidimensional theory for multilayered anisotropic plates should account for the effects of transverse shear deformation. In this respect, the axiomatic approach has been employed either by assuming the distributions of the transverse shear stress uXxzand u,,[4-61, or by assuming a displacement field[7141. In Ambartsumyan’s shear deformation theory[4] for specially orthotropic, symmetric laminates, the distributions of the transverse shear stresses crxrand u,.~ in each layer are assumed to be known. Whit-
790
M. DI SCIUVA
hand, the finite element method provides a convenient and powerful technique for analyzing complex structures. There has been considerable effort of late directed toward the development of shear deformable. multilayered anisotropic plate finite elements[lS-191. As a basis for the development of a displacement plate element of this type, the previous authors have employed the shear deformation plate theory developed by Yang et a1.[7]. This approach requires an arbitrary shear correction factor which is dependent on the lamination parameters for obtaining Improvements of this formulation have been accurate results: furthermore. it does not model the warpage of the cross-section, which is very strong suggested by Srinivas[9] by assuming the previous ratio values. displacement field for each layer. High-order the- for low span-to-thickness Therefore, it seems desirable to develop a finite ories have been developed by Lo et al.[ lo] and Murelement displacement approach based on a refined thy[ll]. shear deformation plate bending theory which alFor plates consisting of layers perfectly bonded lows one to model, among other deformation states, together (the layers function concurrently without slippage), we know from elasticity theory that, at the warpage of the cross-section. The purpose of this article is to develop. on the the interface of the kth and (k + 1)th layers, the basis of the multilayered anisotropic plate theory following conditions must hold? proposed in [13], a rectangular flat finite element which includes the distortion of the deformed nor= bI’k+, M’k l’k = T’k+ I l(k = llk+l mal to the reference surface, the transverse shear (2) effects, as well as the coupling effects of stretching and bending. u;:,k = u;;.k+ I. uyz.k = u_vz.k+I Uxz.k = ux.-.k+ I Some preliminary numerical results show that (3) the developed finite element is very efficient in predicting gross behavior of thick and thin multilayered From eqn (2) it follows immediately that anisotropic plates.
ney[51 has extended this theory to symmetric laminates with arbitrary orientations. Refinements of the Ambartsumyan’s theory, including all transverse effects in a consistent manner, are presented by Voyiadjis and Baluch[6]. In the shear deformation theory developed by Yang et a/.[73 and applied to the specific problem of bending and flexural vibrations by Whitney and Pagano[8], the following displacement field is assumed for the whole plate behavior
PLATE
Continuity conditions cannot be imposed on the remaining stresses and strains; in general, these are discontinuous at the interfaces. Since the transverse shear strains derived by the displacement field given by eqn (1) are constant across the thickness, it is not difficult to realize that the contact conditions on the stresses cannot be satisfied if the layers have different mechanical properties. In a series of papers[l2-141, the writer has proposed a linear displacement model which allows the contact conditions to be satisfied simultaneously. Such a displacement model has been utilized to analyze the flexural behavior of the multilayered orthotropic[ 121 and anisotropic[ 13, 141 plates. Numerical tests carried out on the cylindrical bending of an infinite symmetrically laminated strip have demonstrated the accuracy and efficiency of the proposed model. Because of the complexity of the theories involved, the investigations discussed above are severely
limited
in their
applications.
Displacement
EQUATIONS
field
The laminate under consideration has constant thickness h and is composed of a finite number N of thin orthotropic layers of uniform thickness perfectly bonded together. The principal axes of elasticity of any individual layer are assumed to be oriented arbitrarily with respect to the laminate axes, and the reference surface is superposed on the z = 0 plane (see Fig. I). The material properties and the thickness of each layer may be entirely different. The constitutive relations for any individual layer are given byS
On the other
where the C, are the elastic coefficients t In the following the index k (k = 1, N - 1) will refer to interfaces and the index s (S = 1, N) to layers. N being the number of layers.
and the
$ It must be noted that, in the given constitutive relations, only nine of the elastic coefficients are independent.
791
Development of an anisotropic, multilayered, shear-deformable rectangular plate element
faces (for the undefined symbols, see Nomenclature) . In examining the relations (7)-(9), it is not difficult to realize that the displacements u and v are continuous functions of the z-coordinate for all values of the (pp(x, y) and x&r, y) functions. It follows that the expressions for (Pk and Xk may be found from the contact conditions for the transverse shearing stresses. Once these conditions have been satisfied, it is found that eqns (7)-(9) can be rewritten as1131
Fig. 1. Coordinate system and numbering of the layers and their interfaces.
x ?J
=
u”
+
dy,
-
@v,
+
(z x
k
usual notation for stresses and engineering strain components has been adopted. As usual, the normal stress oZZis assumed to be small in comparison with other normal stresses and is neglected. Taking into account this assumption, the constitutive relations assume the following contracted form
x w
=
(Z
-
zk)Yb
(dkrk
-
-
+
zk)
(7a)
Zk)
(84
bky,)
Zk)Y(Z
-
HP
@a)
where ak, bk, ck, and dk are known constants only depending on the mechanical properties of various layers. Strain
energy
The strain energy per unit area of the plate is given by or
@ = g Q(s) s=l
(10)
where and
@‘(s)= ;
j-:‘,{dr
(11)
{E}dz
(6) stands for the strain energy per unit area stored in the sth layer and
or Iat) = [Q2(s)l k,I
64
where Qij = Cij for i, j = 4, 5 and Qij = C, Ci3Cj3IC33 for i, j = 1, 2, 6. The definition of the matrices introduced in the relations (5) and (6a) follows directly from the relations (5) and (6). The displacement field is assumed to be of the following form[ 131
By making use of eqn (11) in conjunction with the eqns (7)-(9) and the constitutive relations, it is not difficult to show that eqn (10) assumes the following matrix form @ = 1 {6}r[E]{6}
(lOa)
where U = U0 + z(rx - $!r) + z (Pk(Z - Zk)Y(Z - Zk) k
W (7)
7.J= Y0 + z(y, - $,)
+ 2
Xk(Z
-
zk)y(z
-
= LupxG!X& $V b % YXJ Y.v.xYX.?Y?..” -
0 WJX
-
0 w,yy
-
w!ryJ
.
Zk)
Explicit expressions for the elements matrix are given in the Appendix.
k
of the [El-
(8) w
=
w”.
(9)
Here Y(z - zk) is the Heaviside unit function and the summation is extended over the N - 1 inter-
FINITE ELEMENT FORMULATION
As is known, the accuracy of the finite element solution depends on the ability of the assumed in-
M. DI SCILJVA
792
terpolation functions to accurately model the deformation modes of the structure. Since the coupling between extension, bending, and transverse shear deformation states affects the static and dynamic responses of the laminated anisotropic plates, a correct finite element displacement formulation should model such a coupling. The results obtained by Prior and Barker[lS] show that a flat rectangular plate element with seven degrees-of-freedom (DOF) at each nodal point should be accurate for this purpose. The nodal parameters consist of the three generalized displacements u’, Y’, and N”: the two shear rotations yr and y,., and the two total rotations WY)= 4,Yand wpY = 6,. This displacement model fully satisfies the inter-element compatibility conditions. However, it has a great disavantage in that it does not include the constant twist term. In order to improve the formulation, we have added the mixed second derivative uzp,,= 6,X,; thus, the proposed rectangular plate element has eight DOF at each nodal point, u”, zf”, w”, 6,, 6?, 6,, yI. yy , giving a total of 32 DOF. Linear polynomials have been used for approximating the inplane displacements [lo, zl”, and the shear rotations yX and Y?; cubic polynomials have been used for approximating the transverse displacement MI’. Let us call the nodes as in Fig. 2. It can be shown that the following expressions satisfy the inter-element compatibility and the constant strain conditions[21] (u stands for 11’. zf”, yX, y,,):
and similarly for a’)‘(q), H!(q), etc., replacing 5 with q and a, with b, in above. The vector (6)’ comparing into eqn f 10al can be written in terms of the nodal degrees of freedom by taking the appropriate derivatives of the interpolation polynomials [eqns ( 1211. When this is given, we obtain
where {cl}, is the vector of the degrees of freedom of the nodes surrounding element e, and [N] is a rectangular matrix (13 x 32) whose elements are products of the polynomials H by itself or by their derivatives. The strain energy of the element can be determined by substituting expression (13) into eqn (lOa) and integrating over the area of the element. The result is
where
is the rectangular plate element stiffness matrix (32 x 32). External loadings and boundary conditions are subject to the usual finite element analysis considerations. NUMERICAL EXAMPLES
In order to assess the accuracy and efficiency of the developed rectangular finite element, some preliminary numerical tests have been carried out. This evaluation was conducted on two sample problems, for which analytical solutions are available. In this work, for a definition of analytical solutions we where (5 = xla, and q = y/b,) mean the values obtained by solving exactly the system of the differential plate equations based on Hm = 1 - 5; H6Ykc) = 5; HIaS) = 1 - 35’ the displacement field given by eqns (7a)-(9a). + 26’; H6:‘(5) = S2(3 - 25) The first example is the well-established cylindrical bending problem of cross-ply and angle-ply fi:‘(Sl = a,W - 2[ + 11: @t?(Q = a,S*(E - 11 laminates; the second one is the bending of a symmetric, cross-ply square plate simply-supported on (2,2) all edges. (1,2) In addition to comparing the approximate solution with the analytical one. the convergence of the discretization error was also examined numerically . In the following, we consider a three-layered, symmetric laminate with layer properties a,
y~q~z be[ (1,lI
I
5
(2,l)
w
EJEr G&ET
= 25, = 0.2,
GLrlET
= 0.5,
VLT = VT7 = 0.25.
X
Fig. 2. Global and local coordinate systems and node numbering for rectangular olate element.
where L and Tare the directions parallel and normal to the fibers, respectively, and vLTis Poisson’s ratio
Development
measuring transverse parallel to the fibers.
of an anisotropic,
multilayered,
strain under uniaxial
shear-deformable
CROSS-PLY
stress
r0’/90”/0’,:
AnalytIcal 0
2.0
Cylindrical bending
Flnlte
Sotu tlon
Element
~WWF.E.S.
-
WA.S.YWA.S.
where Wr.n.s. stands for finite element solution and WA.s, is the analytical solution. It seems that in both laminate configurations the convergence is very rapid. Figure 4 gives the analytical and finite element values for W as a function of the plate thickness ratio
Cl21
Solution
/
1.0
SolutionC13.141
X Finite
Element
//
I’
ANGLE-PLYt30’/-30’/30): --Anatytlcal
The plates considered in this example are of infinite length in the y-direction, uniformly supported along the edges x = 0, a and subjected to a sinusoidal distributed transverse loading p;(x) = p. sin Ax (A = r/a) on the reference surface. The cross-ply laminate (0”/90”/0”) has three layers of equal thickness with the fibers in the outer layers oriented in the x-direction (0”) and those in the inner layer oriented in the y-direction (90”). The angle-ply laminate (30”/- 30”/30”) has ply thickness (h/4, h/2, H/4), respectively, with fibers in the outer layers oriented at +30” and those in the inner layer at - 30”, with respect to the x-axis. Under these conditions, the deflected surface is cylindrical; therefore, only a strip of the plate needs to be modeled for the computational purpose. Furthermore, due to symmetry, one only has to model a half span. Figure 3 shows the convergence patterns of the nondimensional transverse displacement, W, at the ratio) center of the span for p (span-to-thickness equal to 10. In this figure, as well as in Fig. 5, we have posed %E =
-
2.2
193
rectangular plate element
/
Solution
I’
/’ /
1.4 1.2
t t
I
I
6
I
8
I
10
I
12
I
14
I
16
I
16
I
2
Fig. 4. Typical results for the center deflection of a threelayered, symmetric, simply-supported laminate in cylindrical bending.
For both laminate configurations, and for all values of the thickness-to-span ratio, the approximate solution is very close to the analytical one given in Refs. [12-141. Close agreement between numerical results and analytical solutions is also obtained in the distribution of the in-plane displacements and flexural stresses across the plate thickness. Owing to this agreement, and also observing the fact that the form of the various distribution curves correspond closely to those given in [ 12-141, these curves are not shown here. Simply-supported
square plate
The numerical results given in this section refer to a three-layered, symmetric cross-ply, square
h/a.
The plotted finite element values were obtained with a 5 x 1 mesh which modeled the behavior of a half span.
1
IE
-10 $I
-20 M
CROSS-PLY
(0390”/0”,
X--3(
ANGLE-PLY(30”/-JO”f3@“. -30
-10
-50
N”lnb.x
Fig. 3. Cylindrical bending of an infinite three-layered, symmetric, simply-supported laminate under sinusoidal transverse loading POsin Ax. Maximum deflection error as a function of the number of elements.
of etementr. NE
Fig. 5. Bending of a three-layered, symmetric cross-ply (0”/90”/0”), simply-supported square plate under sinusoidal transverse loading PO sin Ax sin Xy (A = n/a). Maximum deflection error as a function of the number of elements.
M. Dt
f%XJVA
P
No”dimenlio”al
in-plane
displacement.
0 = ET”Lo,a/2,r,
/hp,
h
._ ; I B ;
Fig. 8. Variation of u through the thickness of a threelayered, symmetric cross-ply (0”/90”/0”), simply-supported square plate under sinusoidal transverse loading po
-.2 2
4
6
0
10
12
14
16
10
I
Fig. 6. Maximum deflection of a three-layered, symmetric cross-ply (O”/W/o”), simply-supported. square plate under sinusoidal transverse loading p. sin Ax sin hy (A = P/U).
plate simply-supported on all edges and subjected to a sinusoidal distributed transverse loading &(x, y) = p sin (TX/U) sin (phyla) on the reference surface. The laminate has three layers of equal thickness, with the fibers in the outer layers oriented in the xdirection (0”) and those in the inner layer oriented in the y-direction (90”). The analytical results are obtained by choosing u”, u’, w”, y,, and yv as u’ =A,cosAxsinhy: u” = A,,sinhrcoshy; M’O= A,, sin Ax sin hy ; (A = 57/u) ?/x = A; cos Ax sin Ay ; yv = A; sin Ax cos Ay
which satisfy the differential equilibrium equations and the simple support, without in-plane normal restraint boundary conditions given in [12]. Figure 5 shows the convergence pattern of the nondimensional deflection ti; at the center of the
-10
-12
CROSS-PLY ---Elasticity -.‘A
Thin Flnile
-4
plate for p = 10. Because of the symmetry, only one quadrant of the plate was modeled and NE stands for the number of the elements for a quarter of the plate. Fast convergence to the analytical solution is observed: for the 5 x 5 mesh the error in the approximate results is not higher than 2%. Thus, this mesh has been used to obtain the results given in the following. Figure 6 compares the exact, analytical, and finite element values for M:as a function of the plateto-thickness ratio h/u. The distribution of the flexural stress a,, as resulting from the present finite element approach is given in Fig. 7 in conjunction with those resulting from the elasticity solution[2], the approximate bidimensional theory developed in [ 121, and the classical lamination theory (C.L.T., Kirchhoff thin plate theory)[201. The distributions of the in-plane displacement U are plotted in Fig. 8. Close agreement of the analytical and finite elements results, with the exact elasticity solution]21 being obtained. As expected, at very low span-tothickness ratio values the thin plate theory is in-
-2
(0~/9O~/V): Solu~lon
Adytic.l
-6
-13
sin !LXsin Ay CA= T/U).
;
/I/
sd”tlon
Plate Theory Element
(C.L.T.)
Solution
INTERFACES
Nondimensional
bending
stress,
zXX
= bXX(a/2,a12,z)
Ip,
Fig. 7. Variation of c,, through the thickness of a three-layered, symmetric cross-ply (0”/90”/0”) simplysupported square plate under sinusoidal transverse loading po sin AXsin hy (A = n/n).
79.5
Development of an anisotropic, multilayered, shear-deformable rectangular plate element adequate to model the correct behavior of the multilayered composite plates.
deformationtheoryformultilayeredorthotropicplates. Proc. of the VII Congress0 Naz. AIDAA, Oct. 1983, Naples, Italy. In press in 1’Aerotecnica Missiii e Spazio.
13. M. Di Sciuva, A refinement of the transverse shear deformation theory for multilayered anisotropic
CONCLUDING REMARKS
Based upon a refined approximate shear deformation theory proposed by the author in a series of previous papers[ 12-141, an effective multilayered anisotropic rectangular plate element is developed. The finite element is a rectangle with 32 degrees of freedom which include stretching, bending, and transverse shear deformation states. Numerical tests carried out on two sample problems show that the element is very efficient in predicting responses of thick and thin laminated plates under static loading. also at very low span-to-thickness ratios, and models the warpage of the crosssection, which has not been found to be feasible by other conventional finite element displacement approaches[lS191. It is believed that the approach can conveniently be extended to perform static, vibration, and buckling analyses for a large class of laminates with varied loading and boundary conditions. This is in the planning for a future article. Acknowledgment-This
research has been worked out by the financial support of the Italian Ministry of Education (MPI).
plates. Department
(1981).
18, J. N. Reddy and W. C. Chao. Nonlinear oscillations
Mater. 4, 20 (1970).
3. N. J. Pagano, Influence of shear coupling in cylindrical bending of anisotropic laminates. J. Compos.
Vibration of thick rectancrular mates of bimodulus composites materials. J. Aipl. iech. 48, 371 (1981). 20 R. M. Jones, Mechanics of Composites Materials. McGraw-Hill (1975). 21 F. K. Bogner, R. L. Fox and L. A. Schmit, The generation of inter-element-compatible stiffness and mass matrices by the use of interpolation formulas. AFFDL-TR-66-80, Proc. of Conference on Matrix Methods in Structural Mechanics, p. 397 (1965).
Mater. 3, 534 (1969).
6. G. Z. Voyiadjis and M. H. Baluch, “Refinements in the Bending of Plates with one Plane of Elastic Svmmetry”. P&c. of the Intern. Symposium on the Mechanical Behaviour of Structured Media, 507. May 1981, Ottawa, Canada. P. C. Yang, C. H. Norris and Y. Stavsky, Elastic wave propagation in heterogeneous plates. Int. J. Solids and Struct. 2, 665 (1966).
J. M. Whitney and N. J. Pagano, Shear deformation in heterogeneous anisotropic plates. J. Appi. Mech. 37(4), 1031 (1970). S. Srinivas. A refined analysis of composite laminates. J. Sound and Vib. 30(4), 495 (1973).
10. K. H. Lo, R. M. Christensen and E. M. Wu, A highorder theory of plate deformation. Part 2: Laminated plates. Trans. ASME, Ser. E, J. Appl. Mech. 44, 669 (1977). 11. M. V. V. Murthy, An improved transverse shear deformation theory for laminated anisotropic plates. of the transverse
By following the same procedure outlined in [12], it is readily shown that eqn (11) takes the compact form
where
shear
101 lEl3O)l l-&(S)1 [Ezz(s)l WI 101
]Ellts)l
[E(s)1 =
Plates.
Translated from the Russian by T. Cheron and edited by J. E. Ashton, Technomic Publishing Co. (1969). 5. J. M. Whitney, The effect of transverse shear deformation on the bending of laminated plates. J. Compos.
NASA TP I903 (1981). 12. M. Di Sciuva, A refinement
rectangular plates. J. Appl.
of laminated, anisotropic,
Mech. 48, 396 (1982). I,~ *,, C. W. Bert, J. N. Reddy, W. C. Chao and V. S. Reaay, IY.
APPENDIX
1. N. J. Pagano, Exact solutions for composite laminates in cylindrical bending. J. Compos. Mater. 3, 398 (1969). 2. N. J. Pagano, Exact solutions for rectangular bidirectional composites and sandwich plates. J. Compos.
Theory of Anisotropic
Engineering-Poly-
M. Di Sciuva, A refined transverse shear deformation theory for multilayered anisotropic plates. In press in Atti Accad. delle Scienze di Torino (Italy). IS. C. W. Pryor and R. M. Barker, A finite-element analysis including transverse shear effects for applications to laminated plates. AIAA J. 9(5). 912 (1971). 16. A. K. Noor and M. D. Mathers, Shear-flexible finiteelement models of laminated composite plates and shells. NASA TND-8044 (19751. 17. .I. N. Reddv and W. C. Chao. Laree deflection and large ampliiude free vibrations’ of laminated composite-material plates. Comput. Structures 13(2), 341
REFERENCES
Mater. 4, 330 (1970). 4. S. A. Ambartsumyan,
of Aerospace
,4. technic of Turin (Italv)-ReDort n.5 (1983).
[E33(s)l
sYk.?&C
[E34(S)l LE44Cs)l I
with
lE130)l = z[Ei~(sIl
[EIICSII= [AITIQ,WllAl: + lAITIQi(s)llB*l [EM(
= zlAITIQ,tsIl;
[E340)1
= z([E14(~)1 x [,544(s)]
LE33b)I
= z*[EIIWI
+
[W[Q,b)i
[&z(s)1 = lQz(IIlt1Il + [A*])
+ [B*lTIQ~(sN); = zz[Q,(s)l +
zt[Al*lQ,(s)l IB*l
w*il? + [B*i%21(~~i
LB*].
The undefined matrix symbols which appear in the previous relations stand for rl000i
[A]=
0001 ; I ni in I
TT T7
[A*]=
;; 1
J
M. Dr SCIUVA
796
+ 2D& E(10, 10) = 022 + 2042 + 2Dk, + D 2626 + DG.; E(9, 12) = Dxj + Dh + D$; E(l1. II) = D,,
where
--__
(4, B. c, 0) _ _ (A, B, C, D) 1
a-1
E(9. 13) = Da + D& + D$,; E(l0. 11) = D,z + D$ + Dfe; E(10,12) = Dzz + D$ + Dk:
Substitution of the relation (1A) into eqn (IO) gives the following expressions for the elements of the symmetric matrix [El [see eqn (lOa)]: E(1. 1) = A,,: E(1,2) = E(l, 3) = A,Q(l, = Aiz; E(1. 7) = B,l + By, + B,h,
4)
E(1, 8) = 816 + Bib + BY,; E(1. 9) = 816+ BQ6 + BP2; E(l, 10) = B,? + B?z + BfCj: E(1, 11) = B,, E(1, 12) = B,z; E(1, 13) = 816: E(2. 2) = E(2, 3) = AM; E(2, 4) = Ax:
E(2,7) = B,c, + Bik + B&6; E(2, 8) = Be6 + B$ + IS;
E(2, 9) = Bfjb + I?&, + B$:
E(10,13) = 026 + D$ + D&i; E(l1, 12) = D,*; E(11, 13) = Dlh: E(12, 12) = D22; E(12, 13) = Dxj: where [Ail, B,j, Dijl = < 11, Z, ~‘1 Qtj’xs) > for i, j = 1, 2. 6 and A;j= hQij(l) for i, j = 4, 5 &I&2&6&%&6
E(4, 9) = BZb + B% + B$; E(4, 10) = B22 + B$* + B&: E(4. 11) = B,z: E(4, 12) = &*; E(4. 13) = Bx,: E(5,5) = AM + A& + A&; E(6,6) = Ass + A& + A&; E(5. 6) = Aa + At + Ai4 = A45 + A% + A&: E(7. 7) = D,, + 2DQ, + 2D$ + DQf
+ 0:: + 207:;
=
<[I:
B%2Bf2&‘6B:6@6
~Q,,~~~Q,~~~~Q,~~~~Q~6oa66(s)1AA ’
1
=
[ D$2Dt2DtbD46D$
<:
[Q22(~)Q,z(~~Q,~(~)Q26(s)e6603 Bk ’
[1
&,BT2&6&sB& [ D~IDTJWX.S~&
<[l;
E(3, 9) = 866+ B&i+ B&k: E(3, 10) = 826 + B% + B&6; E(3, 11) = B,c; E(3. 12) = 826 E(3, 13) = 866; E(4, 4) = Azz; E(4, 7) = B,2 + By* + B$; E(4, 8) = BZ6 + B& + BC*:
1
LDfl,D?zDf16Dy6D&
E(2. 10) = 826 + B$ + B&; E(2, 11) = 816; E(2, 12) = 826; E(2, 13) = I&,$,; E(3, 3) = A6,5; E(3, 4) = Ax,; E(3, 7) = B,6 + B;16 + B&6; E(3, 8) = B@ + B$ + Bie;
E(13. 13) = DM,
1
=
[Q,,(s)Q,~(s)Q,s(s)Q26(S)e6601
c/, ’
I DA > <[If [Qzz~~)~I~~~)~
&2Bi’2Bi’6&!6&7
=
L&2Df2Df6’6D$D$
[D8DiWP,“D~l
= < [Q,,(s)Q,~(s)Q~~(s)Q~~(~)I &‘b
’
U%D~:D~iX%l
= < [Qzz(s)Q,s(s)Q~~(s)Q~~(~)~ &B, ’
[DYiD(‘6D%D%I = < [QI ~(S)Q,~(S)Q~~(S)Q,(S)~ ckc, > [D%Di’$!D#D&‘] = < [Qz~(S)Q,~(S)Q~~(S)Q~~O~DkD, > [DfiDYtDtiD@l
= < [Q,z(s)Q,s(s)Q26(s)Q6601 &Br ’
E(7,8)
= D,6 + Dyh + D& + Di', + D& + Dyi + Dd + D$ + Den b“; E(7, 13) = D,6 + D& + 0th:
[D3Di’kD%D~‘l = < [Q,,(s)Q,b(s)Q26(s)Q66(s)l A&r >
E(7.9)
= D,c, + 2D& + D$ + D& + Dy: + D(;f + D&' + D@; E(8. 8) = 0% + 2D$ + 2D&, + 0';; + 2Df; + Ott,";
[Db’D%D&D&] = < [Q,2(s)Q,b(s)Q26(s)Q6601 B&r >
E(7. 10) = D,2 + D$ + D& + Dye + 0% + Dr$' + D?: + 046"+ D&; E(7. II) = D,, + Dy, + D$,;
[DdDi%%@d]
= < [Q,z(s)Q,h(s)Q26(s)Qb6(s)l A/A
[D%D%!D&!D&‘]= < [Q~z(s)Q,~(S)Q~~(S)Q~~(S)~ BkDr > [DFfDf6do’i66D& = < [Q,z(s)Q,~(s)Q~~(S)Q~~(~)~ CkDr >
E(7, 12) = D,? + D& + D&6: E(8. 9) = Dhb + D&
+ D&, + Die + D$*,+D Ebb + DQ; + 0% + D,z. E(8, 10) = Dzb + 2D& + Dyz + D&, + D$'k+ Df$ + Dig + D&; E(9. 9) = D66 + 2046 + 2D& + Dg: + D‘$ + 2D$ E(8. 11) = D,6 + Dt6 + DC,; E(8. 12) = Ox, + D$ + Diz; E(8, 13) = Db6 + D& + DL,; E(9, 10) = Dzb + D26 + D$ + D;;, + D$ + 04; + D$ + D$$ + Dg;; E(9, 11) = D,6 + Dye + D:,;
>
For notational convenience,
we have posed
(..) = \$,(..)dr:
G, = x&(z
where t-l =OL= =:;
-z,),