- Email: [email protected]
A high-order shear element for nonlinear vibration analysis of composite layered plates and shells
International Journal of Mechanical Sciences 41 (1999) 461—486
A high-order shear element for nonlinear vibration analysis of composite layered plates and shells O. Attia *, A. El-Zafrany Imperial College of Science, Technology and Medicine, Mechanical Engineering Department, London SW7 2BX, UK School of Mechanical Engineering, Cranfield University, Cranfield, Bedford MK43 0AL, UK Received 30 June 1998
Abstract This paper introduces a family of high-order facetted shell elements for linear and nonlinear stress and vibration analysis of composite layered plate and shell structures. Engineering slope angles are employed in element equations, and transverse stresses are expanded over the thickness. The lateral deflection is modelled by conforming or non-conforming Hermitian shape functions, within rectangular or paralellogrammic and triangular elements. Nonlinear terms associated with geometrical nonlinearity are also derived using a practical approach based upon the actual components of strain. A finite element programming package was designed employing the newly developed elements. Several case studies have been investigated and package results were compared with existing theoretical and/or experimental results. It has been proved that the developed elements can lead to accurate estimations of natural frequencies. The effect of fibre angles on natural frequencies has also been investigated with some case studies, and the results proved that the package can be a useful tool for the design optimization of composite layered plate and shell structures. 1998 Elsevier Science Ltd. All rights reserved. Keywords: High-order shear element; Nonlinear vibration analysis; Composite layered plates and shell
Notation a, b A B D D
side lengths of a rectangular plate matrix of displacement derivatives strain-displacement matrix plate flexural rigidity stress-strain matrix
* Corresponding author. 0020-7403/99/$ — see front matter 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 4 0 3 ( 9 8 ) 0 0 0 7 7 - 0
462
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
E F F,G,H G G G h I K l, m, n M N G u, v, w x, y, z x , y , z c d e h ,h V W q '
, t u k l o r
Young’s modulus nodal loading vector Hermitian shape functions plate thickness unit matrix stiffness matrix directional cosines mass matrix Lagrangian shape functions displacement components in x, y, z directions, respectively Cartesian coordinates Coordinates based on material principal axes (see Fig. 1) transverse shear strain vector nodal displacement vector strain vector average slope angles in the x and y directions, respectively transverse shear stress vector non-dimensional frequency parameter (,c , c ) transverse shear strains at the plate midplane VX WX natural frequency shear modulus Poisson’s ratio density stress vector
1. Introduction The use of composite materials in automotive and aircraft structures has increased in recent years, due to the high strength to weight ratio of composites and their many other advantages. Fibrous composites can optimally be designed with the tuning of fibre angles in different layers so as to obtain the required properties. Nevertheless, the sensitivity of composite layered plates to shear stresses necessitates the need to estimate transverse shear stresses as accurately as possible. Many industrial applications may also have a geometrical nonlinearity in the elastic range, which should be considered in the analysis. There are a number of finite element publications for composite materials reported in the literature, most of which are limited to static analysis. Epstein and Huttelmaier [1] presented a finite element formulation for the analysis of multilayered thick plates, where the transverse shear and normal strains were considered. Widera and Moumene [2] used an eight-node isoparametric quadrilateral shell element, which included the effect of transverse shear strain in the analysis of laminated plates and shells. In the area of dynamic analysis of laminated structures, Has and Wang [3] presented an analytical approach to the analysis of laminated cylindrical shells. They did not employ the Kirchhoff hypothesis of non-deformable normals by including the inertia terms in the equilibrium equations. Witt and Sobczyk [4] discussed the dynamic response of laminated plates to random
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
463
loading. Pillasch et al. [5] analysed thick laminated plates with a finite element model based on a three-dimensional element rather than plate theory. The assumption that displacements are linear functions of z, the coordinate in the thickness direction, has proved to be inadequate for predicting gross laminate response. Therefore, many higher-order theories incorporating transverse shear strains, have been proposed for better accuracy. Nelson and Lorch [6] used up to second-order terms of z in their displacement field. Reddy [7, 8] introduced a refined nonlinear theory of plates with transverse shear deformation, and a simple higher-order theory for laminated composite plates. The deduction of two-dimensional theories, from the three-dimensional theories, in a transversely isotropic body was presented by Wang [9]. Chang et al. [10] presented a Lagrangian formulation for large translation and rotation analysis of deformable rectangular homogeneous plates. They included the coupling effects between bending and stretching, in the finite element analysis of homogeneous plates. Wung and Reddy [11] presented a first-order shear/fourth-order transverse deformation theory for laminated composite plates. A recent theory was introduced by Tessler [12], who derived a two-dimensional laminated theory for linear elastic analysis of thick composite plates with the equivalent single-layer assumptions for the displacements, transverse shear strain, and transverse normal stress. Cho and Parmerter [13] developed a three-node non-conforming triangular bending element, with five degrees of freedom, for symmetric laminated composites. The derivation of this element is based on a higher order plate theory evolved earlier by the same authors. Many of the high-order element derivations presented in the literature are based on hypothetical nodal parameters, which are difficult to handle with different types of boundary conditions. This paper introduces a family of high-order facetted shell elements for linear and nonlinear stress and vibration analysis of composite layered plate and shell structures. Engineering slope angles are employed in element equations, and parabolic distributions of transverse shear strains over the thickness are considered. The lateral deflection is modelled by conforming or nonconforming Hermitian shape functions, within rectangular or parallelogrammic and triangular elements. The degrees-of-freedom have a maximum total of nine per node, but they are automatically reduced for coplanar nodes, so as to make it easy and efficient to use the elements for the analysis of box structures and curved shells. Non-linear terms associated with geometrical nonlinearity are also derived using a practical approach based upon the actual components of strain.
2. Deformation theory 2.1. Transverse strain modelling Consider a composite layered plate, at an instant of time t consisting of a number of orthotropic layers N , as shown in Fig. 1. Let the midplane of the plate be the Cartesian x-y plane, and the total J thickness at any point (x, y) on the midplane is h (x, y). Hence, the equations of the lower and upper surfaces of the plate are: z "!h/2 and z "h/2, respectively. If the values of c and c represent VX WX
464
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Fig. 1. Composite layered plate.
the transverse shear strains at the midplane, then from the boundary conditions the following values of transverse shear strains are defined as follows: at z "!h/2, and at z "h/2: c "c "0 VX WX at z "0: c "c , c "c . VX VX WX WX The composite layers are considered to be firmly bonded together, hence the distributions of displacement components represent continuous functions with respect to the z-direction. Considering also the infinitesimal shear strain—displacement equations: *u *w c " # VX *z *x
(1)
*v *w c " # WX *z *y
(2)
the transverse shear strains can be approximated as continuous functions with respect to z. Hence, the distributions of the transverse shear strains across the plate thickness can be deduced by applying Lagrangian interpolation to their given values, i.e. c (x, y, z, t),c (x, y, t) [1!4z/h] VX VX c (x, y, z, t)"c (x, y, t) [1!4z/h] WX WX
(3) (4)
2.2. Displacement and velocity components Substituting from Eqs. (3) and (4) into Eqs. (1) and (2) and integrating over the thickness, assuming that the lateral deflection w can be approximated to be independent of z, then
4z *w u(x, y, z, t)"u (x, y, t)!z #c (x, y, t) z! VX 3h *x *w 4z v (x, y, z, t)"v (x, y, t)!z #c (x, y, t) z! WX *y 3h
where u and v represent the values of u and v at z"0.
(5) (6)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
465
The corresponding velocity components can be obtained by differentiating displacement equations with respect to time, i.e.
*wR 4z uR (x, y, z, t)"uR !z #cR z! VX *x 3h *wR 4z vR (x, y, z, t)"vR !z #cR z! WX *y 3h
(7) (8)
2.3. Strain components The x-y strain components will be considered finite and by using a Lagrangian frame of reference, they can be expressed in terms of displacement components [14]: *u 1 e" # V *x 2 *v 1 e" # W *y 2
*u *v *w # # *x *x *x
*u *v *w # # *y *y *y
*u *v *u *u *v *v *w *w c " # # # # VW *y *x *x *y *x *y *x *y
(9) (10) (11)
Defining the following strain vector: e "+e e c , VW V W VW then it can be deduced from Eqs. (9)—(11) that
(12)
e "e #e VW * where *u *x *v e " *y *u *v # *y *x
(13)
(14)
which represents the infinitesimal terms, and
1 *u 1 *v 1 *w # # 2 *x 2 *x 2 *x 1 *u 1 *v 1 *w # # e " * 2 *y 2 *y 2 *y *u *u *v *v *w *w # # *x *y *x *y *x *y
which represents the effect of finite strains.
e* V , e* W e* VW
(15)
466
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Using Eqs. (5) and (6), the infinitesimal strain vector can be expressed as follows: e (x, y, z, t)"e (x, y, t)!ze (x, y, t)#f (z) e (x, y, t) K @ Q
(16)
where *u *x *v e " K *y *u *v # *y *x
*w *x *w , e" @ *y *w 2 *x *y
*c VX *x *c WX , e" Q *y *c *c VX# WX *y *x
(17)
and f (z)"z!4z/(3h) . Substituting from Eqs. (5) and (6) into Eq. (15) it can be proved that: 1 e*" V 2
*u *v *w # # #z *x *x *x
*w *w # *x *x *y
#f (z)
*c *c *u *w *v *w VX # WX !2z # *x *x *x *x *x *x *y
!2zf (z)
*v *c *c *w *c *w *u *c VX VX# WX # WX #2 f (z) *x *x*y *x *x *x *x *x *x
1 e*" W 2
*u *v *w # # *y *y *y
*w *w #f (z) # *y*x *y
#z
*c *c VX # WX *y *y
*u *w *u *c *v *w *v *c VX# WX # #2 f (z) *y *y*x *y *y *y *y *y *y
!2z
!2z f (z)
*c *w *c *w VX # WX *y *y*x *y *y
.
*u *u *v *v *w *w c* " # # VW *x *y *x *y *x *y !z
(18a)
*u *w *u *w *v *w *v *w # # # *x *y*x *y *x *x *y *y *x *y
(18b)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
!z f (z) #f (z)
#z
*c *w *c *w *c *w *c *w VX # VX # WX # WX *x *x*y *y *x *x *y *y *x*y
*u *c *v *c *v *c *u *c VX# VX# WX# VX *y *x *x *y *y *x *x *y
467
*w *w *c *c *w *w *c *c VX VX# WX WX # #f (z) *x *x *y *x*y *y *x *y *x *y
(18c)
3. Finite element derivation 3.1. Displacement and velocity interpolation The composite-layered plate element is defined in terms of the midplane (z"0), where the nodal parameters are also defined. The geometry of the plate midplane is defined in terms of nodal coordinates with respect to a Cartesian x—y—z system of coordinates. Most of the derivatives described here are with respect to a local system of axes, where the x—y plane is in the midplane, and the z-axis is normal to that midplane. Physically, the element is constructed of a number of layers (N ), similar to that shown in Fig. 1. J The l layer’s lower and upper surfaces are described by the following equations: z"ZJ (x, y), z"ZJ (x, y) * 3 and the layer thickness is hJ(x, y)"ZJ! ZJ . 3 * The layer has its own material axes, and the stress—strain relationships are defined with respect to the element local x, y, z axes, as follows: rJ V sJ rJ , rJ "D(l) e , sJ, VX "lJ VW VW W sJ WX sJ VW
c VW . c WX
(19)
The matrices D(l ) and l(l ) are stress—strain matrices rotated from material axes to element axes. The midplane displacement components in the x and y directions will lead to continuous displacement components, if they are represented by at least C-continuous functions. Hence, Lagrangian interpolation is suitable for their modelling, and for an n-node finite element: L L u (x, y, t)" u (t) N (x, y), v (x, y, t)" v (t) N (x, y) G G G G G G
(20)
where u , v are the values of the x—y displacement components, at the midplane node i, and time t, G G and N (x, y) represents a Lagrangian shape function. G
468
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Similarly, the transverse shear components, at the midplane can be interpolated in terms of Lagrangian shape functions as follows: L c (x, y, t), (x, y, t)" (t) N (x, y) , VX G G G L c (x, y, t),t (x, y, t)" t (t) N (x, y) WX G G G where
(21a) (21b)
(t),c (x , y , t), t (t),c (x , y , t) (22) G VX G G G WX G G The lateral deflection w must be C continuous over the plate, and it should be interpolated by means of a Hermitian interpolation. Two types of interpolation, non-conforming and conforming, based upon El-Zafrany and Cookson [15, 16] are employed, where the lateral deflection is interpolated as follows: L w(x, y, t)" +w (t)F (x, y)#w (t) G (x, y)#w (t) H (x, y)#w (t) Q (x, y), , (23) G G G V G G W G G VW G G where w , w , w , and w represent the values of w, *w/*x, *w/*y and *w/*x*y at node i and G G V G W G VW time t, and F , G , H and Q are Hermitian shape functions. Notice also that Q does not exist for G G G G G the non-conforming Hermitian interpolation. From the previous analysis it is clear that the displacement parameters required at each node i are u , v , , t , w , w , w , w , and a practical representation of nodal values will be described G G G G G G V G W G VW in Section 3.8. To simplify the derivation, the local nodal displacement vector for an n-node element will be defined at time t as follows: d(t)"+d (t) d (t) d (t), K @ Q where
(24)
d (t),+u (t) v (t) u (t) v (t) 2 u (t) v (t), , (25) K L L d (t),+w (t) w (t) w (t) w (t) 2 w (t) w (t) w (t) w (t), , (26) @ V W VW L L V L W L VW d (t),+ (t) t (t) 2 (t) t (t), (27) Q L L The velocity components at any point (x, y, z) inside the plate at time t can be represented vectorially in terms of nodal values and shape functions as follows:
q (x, y, z, t),
q (x, y, z, t) VW q (x, y, z, t) X
(28)
where
uR "N (x, y) d (t)!zN (x, y) d (t)#f (z) N (x, y) d (t) q , * K @ @ * Q VW vR
(29)
q ,[wR ]"N (x, y) d (t) X & @
(30)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
469
and
(31)
,
(32)
N 0 N 0 2 N 0 L N (x, y)" * 0 N 0 N 2 0 N L *F *G *H G G G *x *x *x
*Q G 2 *x
*F *G *H G G G 2 *y *y *y
*Q G 2 *y
2 N (x, y)" @
N (x, y)"[2 F G H Q 2 ] . & G G G G
(33)
3.2. Infinitesimal strain components Substituting from Eqs. (20) into Eq. (17), it can be deduced that:e (x, y, t)"B (x, y) d (t) K K K where *N G 2 0 *x B (x, y)" 2 K 2
0
(34)
2
*N G 2 *y
(35)
*N *N G 2 G *x *y
Similarly it can be proved that:e (x, y, t)"B (x, y) d (t) , Q K Q
(36)
and e (x, y, t)"B (x, y) d (t) @ @ @ where
(37)
2
*F G *x
*G G *x
*H G *x
*Q G *x
2
B (x, y)" 2 @
*F G *y
*G G *y
*H G *y
*Q G *y
2
2 2
*F *G *H *Q G 2 G 2 G 2 G 2 *x *y *x *y *x *y *x *y
(38)
470
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Substituting from Eqs. (34), (36), (37) into Eq. (16) then it can be shown that e (x, y, z, t)"B (x, y) d (t)!zB (x, y) d (t)#f (z) B (x, y) d (t) . K K @ @ K Q
(39)
Note also that the transverse shear strain components, as defined by Eqs. (3) and (4), can be represented vectorially as follows:
c (x, y, z, t),
c c (x, y, t) VX "f (z) VX c c (x, y, t) WX WX
(40)
where
4z f (z)" 1! h
(41)
Substituting from Eq. (21) into Eq. (40), it can be deduced that c (x, y, z, t)"f (z) B (x, y) d (t) , Q Q
(42)
where B (x, y),N (x, y) Q *
(43)
3.3. Modelling of finite strain contributions Defining the following rotation vectors:
*u *v *u *u h " K *x *x *y *y * *t * *t h" Q *x *x *y *y
*w *w *w *w h" @ *x *x*y *y*x *y *w *w h " U *x *y
(44)
(45)
(46)
(47)
then it can be proved from Eq. (15) that e "+A h #A h #zA h #f (z) A h * K K U U @ @ Q Q !z f (z) [A h #A h ]!z [A h #A h ]#f (z)[A h #A h ],, @ Q Q @ @ K K @ K Q Q K
(48)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
where
*u *v *x *x
0
*u *v A " 0 0 K *y *y *u *v *u *v *y *y *x *x
A" @
*w *x
*w *x*y
0
0
*w *y *x
*w *y
* *t *x *x
0
0
0
471
0
* *t , A" 0 0 Q *y *y * *t * *t *y *y *x *x *w *x
0
*w *w *y *x *y *w *w *x *x *y
0
*w , A " 0 U *y *w *w *y *x
Using the interpolation equations for displacement components, the following differential expressions can be obtained: dh "G dd , dh "G dd , dh "G dd , dh "G dd K K K Q K Q @ @ @ U U @ where *N G 2 0 2 *x *N G 2 2 0 *x G " , K *N G 2 0 2 *y *N G 2 2 0 *y
G" @
*G G *x
*H G *x
*Q G *x
2
*F G *x
2
*G *H *Q *F G G G 2 G *x *y *x *y *x *y *x *y
G ,N . U @
*F G *y
*G G *y
*H G *y
*Q G *y
(50)
2
*F *G *H *Q G G G G 2 2 *y *x *y *x *y *x *y *x 2
(49)
(51)
2 (52)
472
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Due to the non-linearity of terms in Eq. (48), it will be represented in the derivations as an integration of its differential, and it can be deduced from Eqs. (48) and (49) that: de "A G dd #A G dd #f (z) A G dd #z A G dd * K K K U U U Q K Q @ @ @ !z f (z) [A G dd #A G dd ]!z [A G dd #A G dd ] @ K Q Q @ @ @ K K K @ @ #f (z) [A G dd #A G dd ] K K Q Q K K
(53)
3.4. Element stiffness matrix for infinitesimal strains For the case of infinitesimal strains, it can be deduced from Eqs. (19), (39) and (40) that the corresponding stress vectors at the lth layer are as follows:rJ "DJ [B d !zB d #f (z) B d ] , VW VW K K @ @ K Q
(54)
sJ"lJ [ f (z) B d ] . Q Q
(55)
The strain energy density or the strain energy per unit volume at any point inside the element can be expressed as follows: ºM (x, y, z, t)" e r,e rxy#c s .
(56)
Hence, for the case of infinitesimal strains, the strain energy density at a point inside the lth layer will be: ºM J" [d B !zd B #f (z) d B ] DJ [B d !zB d #f (z) B d ] @ @ Q K VW K K @ @ K Q K K # [ d B f (z) lJ B d ] . Q Q Q Q
(57)
Notice that all the B matrices are functions of (x, y), and only the functions of z and the D matrices, which may be different for different layers, need to be integrated over the thickness. Writing the integration of every function of z over the total thickness in terms of integrations over different layers, then the following expressions of modified D matrices, at any point (x, y) on the midplane, can be deduced from the integration of Eq. (57) with respect to z:
D (x, y)" KK D (x, y)" @@ l (x, y)" QQ
F
,J DJ dz" [ZJ!ZJ ] DJ , VW VW 3 * \F J
F
,J (ZJ)!(ZJ) 3 * zDJ dz" DJ VW VW 3 \F J
F
\F
f (z) lJ dz ,
(58)
(59)
(60)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
D (x, y)"D (x, y)" K@ @K
F
\F F
zDJ dz , V W
473
(61)
f (z) DJ dz, (62) VW \F F D (x, y)"D (x, y)" z f (z) DJ dz . (63) @Q Q@ VW \F Integrating Eq. (57) with respect to the volume of the element, then the strain energy of the element can be expressed at any time t as follows: D (x, y)"D (x, y)" KQ QK
º(t)" + d K d #d K d !d K d K KK K @ @@ @ K K@ @ #d [KV W#KI] d !dR KR d #dR K d Q QQ QQ Q @ K@ K K KQ Q #d K d !d K d !d K d , Q KQ K @ @Q Q Q @Q @ where
(64)
CJCKCLR BK DKK BK dx dy, K@@"CJCKCLR B@ D@@ B@ dx dy, KV W" QQ CJCKCLR BK DQQ BK dx dy, KIQQ "CJCKCLR BQ lQQ BQ dx dy, K " B D B dx dy, K " B D B dx dy , K@ K K@ @ KQ K KQ K CJCKCLR CJCKCLR K " B D B dx dy . @Q @ @Q K CJCKCLR K " KK
The integrations over elements are carried out numerically in a way similar to that employed in two-dimensional elements (see for e.g. Ref. [17]), i.e.
\ \
K " KK
(x, y) B D B J dm dg , K KK K (m, g)
etc., and a double Gaussian quadrature is employed for the numerical evaluation. Notice that Eq. (64) may also be rewritten as follows: º(t)" d (t) K d (t), (65) where d (t) is the nodal displacement vector of the element at time t, as defined by Eq. (24), and K is the element stiffness matrix for the case of infinitesimal strains and can be constructed as follows: K
KK K" !K K@ K KQ
!K K
K@
@@ !K @Q
K
KQ . !K @Q K V W#K I QQ QQ
(66)
474
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
3.5. Element stiffness matrix for finite strain contributions The virtual or differential value of the finite strain contributions, due to a virtual or differential displacement field is de (x, y, z, t)"+de* de* dc* , * V W VW assuming negligible finite strain contributions for the transverse shear strains. The corresponding increment of strain energy density at any point (x, y, z) in layer l is dºM J (x, y, z, t)"de rJ * VW * Substituting from Eq. (53) into Eq. (67) then:
(67)
dºM "+dd G A #dd G A #z dd G A #f (z) dd G A * K K K @ U U @ @ @ Q K Q !z [dd G A #dd G A ]#f (z) [dd GR A #dd G A ] K K @ @ @ K K K Q Q K K !z f (z) [dd G A #dd G A ], rJ @ @ Q Q @ @ VW By direct multiplication of matrices, the following equations can be proved:
(68)
A r ,Sh "SG d , A r ,Sh "SG d , K VW K K K Q VW Q K Q A r ,Sh "SG d , A r ,Sw h "Sw G d , @ VW @ @ @ U VW U U @ where p
0
q
0 VW p 0 q V VW , S " pV qVW S" U q p q 0 p 0 VW W VW W 0 q 0 p VW W and Eq. (68) can be rewritten as follows: V 0
dºM J"dd G SJG d #dd G SJG d #zdd G SJG d #f (z)dd G SJG d K K K K @ U U U @ @ @ @ @ Q K K Q * !z[dd G SJG d #dd G SJG d ]#f (z) [dd G SJG d #dd G SJG d ] K K @ @ @ @ K K K K K Q Q K K K !zf (z) [dd G SJG d #dd G SJG d ] (69) @ @ K Q Q K @ @ Notice that the G matrices are independent of z, but S and S may vary with z even within the same U layer. To simplify the derivations, the average value of the stress vector inside each layer will be employed as a constant parameter within the thickness of the layer, i.e. SJ+ [S(ZJ )#S(ZJ )] * 3 SJ + [Sw (ZJ)#Sw (ZJ )] U * 3
(70) (71)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
475
The following integrated z functions can also be defined:
F
,J S dz" [ZJ!ZJ] SJ , 3 * \F J
S " KK
F
,J S dz" [ZJ!ZJ ] SJ , U U 3 * \F J
S " UU
F
S " @@ S " QQ
zS dz ,
(72)
(73)
(74)
\F F
\F
S "S " K@ @K
(75)
zS dz ,
(76)
4z z! S dz, 3h F
\F
S "S (x, y)" KQ QK S (x, y)"S " @Q Q@
F
(77)
4z S dz. z z! 3h
(78)
\F
F
\F
4z z! S dz, 3h
The corresponding increment of the strain energy is obtained by integrating Eq. (69) over the volume of the element, i.e. dº (t)"dd KN d #dd (KN #KN )d * K KK K @ @@ UU @ #dd KN d !dd KN d !dd (KN )d Q QQ Q K K@ K @ K@ K #dd KN d #dd (KN )d !dd KN d !dd (KN )d K KQ Q Q KQ K @ @Q Q Q @Q @ where
CJCKCLR GK SKK GK dx dy, KN@@"CJCKCLR G@ S@@ G@ dx dy KN " G S G dx dy, KN " G S G dx dy UU U UU U QQ K QQ K CJCKCLR CJCKCLR G S G dx dy, KN " G S G dx dy KN " K K@ @ KQ K KQ K K@ CJCKCLR CJCKCLR G S G dx dy KN " @ @Q K @Q CJCKCLR KN " KK
(79)
476
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Eq. (79) can also be rewritten as follows: dº (t)"dd (t) KNd (t) * where
(80)
KN !KN KN KK K@ KQ (81) KN" !(KN )R KN #KN !KN K@ @@ UU @Q (KN )R !(KN )R KN KQ @Q QQ and the total increment of the strain energy due to total elastic stresses and strains will be dº(t)"dd (t) [K#K p ] d (t)
(82)
3.6. Element mass matrix The kinetic energy per unit volume, at any point inside the element is: KE" o [uR #vR #wR ], o [q q #q q ] VW VW X X Substituting from Eqs. (29) and (30) into Eq. (83), it can be shown that:
(83)
KE" o +[d N !zd N #f (z)d N ] K * @ @ Q * ;[N d !zN d #f (z)N d ]#d N N d , (84) * K @ @ * Q @ & & @ The N matrices, in Eq. (84), are independent of z, and density functions can be defined at any point (x, y) on the element midplane, in a way similar to that employed in Eqs. (58)—(63), i.e. ,J o dz, [ZJ (x, y)!ZJ (x, y)] oJ 3 * \F J F ,J (ZJ (x, y))!(ZJ (x, y)) 3 * zodz" oJ , etc. o (x, y)" @@ 3 \F J Integrating Eq. (84) over the volume of the element, then the kinetic energy of the element at a time t can be expressed as follows:
o (x, y)" KK
F
KE (t)" +d M d #d (M #M ) d @ @@ && @ K KK K #d M d !d M d !d M d Q QQ Q K K@ @ @ K@ K #d M d #d M d !d M d !d M d , , K KQ Q Q KQ K @ @Q Q Q @Q @ where
CJCKCLR oKK N* N* dx dy,
M " KK
CJCKCLR o@@ N@ N@ dx dy,
M " @@
(85)
etc.
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
477
Eq. (85) can also be simplified as follows: KE (t)" d (t) Md (t) where
(86)
d (t)"+d (t) d (t) d (t), K @ Q and M is the mass matrix of the element which is constructed as follows: M !M M KK K@ KQ M #M !M M" !M K@ @@ && @Q M !M M KQ @Q QQ
(87)
3.7. Element dynamic equations Applying the principle of minimum total potential energy at an instant of time t, then ds (t)"dKE(t)#dº (t)!d¼(t)"0
(88)
where d¼ is the work done by the applied forces at time t, and can be expressed as follows: d¼(t)"dd (t) F (t)
(89)
where F (t) represents a nodal loading system equivalent to the actual applied force. Substituting from Eqs. (82), (86), (89) into (88), then it can be deduced that M d® (t)#[K#K p] d (t)"F (t)
(90)
which represents the dynamic matrix equation of the element. 3.8. Practical nodal values and rotated matrices The use of derivatives of w as nodal values causes some difficulties with the specification of boundary conditions for thick plates. Engineering slope angles can be defined so as to lead to displacement components linearized over thickness [18], and they are defined as follows: *w *w h " !0.8c , h "! #0.8c WX W VX V *y *x
(91)
*w h" for conforming elements, "0 for non-conforming elements. X *x*y The average values of transverse shear strains are used as nodal values, and can also be expressed as follows: cN " c ,u , cN "c ,!u , W WX WX V VX VX
(92)
478
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Defining a new nodal displacement vector based on such practical values as follows: d "+ 2 u v w (h ) (h ) (h ) y(u ) (u ) (u ) 2 , LCU G G G VG WG XG VG WG XG then it can be shown that
(93)
d"Pd
, (94) LCU where P is a transformation matrix defined from the previous relationships. The element local axes are defined with respect to structural global axes by means of directional cosines, and it can be shown that d "Rd , (95) LCU E where R is the element rotation matrix, and d is the global nodal displacement vector. Hence, it can E be deduced that: d"PRdg,Ldg
(96)
and the element matrices defined with respect to the global axes can be obtained as follows: K "LKL, E
KN"LK pL, M "LML E E
(97)
4. Case studies Based upon the theory described in previous sections, a modular finite element computer programming package was designed. The package was coded in FORTRAN 77, and tested on VAX/VMS main frame computers, Unix workstations, and DOS/WINDOWS operating PCs. The package consists of several modules, where each module is a file containing a group of subroutines. It can be employed for finite element static and dynamic analysis of plates and shells made of isotropic or composite layered materials. A number of case studies, with published analytical and/or experimental results, have been analysed so as to assess the accuracy and reliability of different aspects of the developed theory. 4.1. Dynamic analysis of isotropic square plate case This case represents a square plate, as shown in Fig. 2, made of isotropic material with the following properties: E"1.28;10 Pa, l"0.333, o"1.5;10 kg/m . The plate is fixed at one edge, whilst the other three edges remain free. The package results were compared with a published analytical solution carried out by Narita and Leissa [19], as demonstrated in Table 1, where non-dimensional frequency parameters were displayed. The nondimensional frequency parameter U is defined as follows: '"ua (oh/D
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
479
Fig. 2. Composite layered square plate.
Table 1 Natural frequency parameters for isotropic square plate Mode
1
2
3
4
5
Conforming Q. Non-conforming Q. Non-conforming T. Analytical Solution [19]
3.461 3.456 3.453 3.459
8.361 8.362 8.393 8.358
21.10 21.34 21.34 21.09
27.04 26.81 26.82 27.07
30.52 30.47 30.67 30.56
It is clear from Table 1 that natural frequency parameters, for different modes of vibration, as obtained by means of different finite elements, agree very well with corresponding analytical solutions. 4.2. Composite layered square plate case This case has overall dimensions similar to the previous case shown in Fig. 2, but it has eight layers of a fibrous composite material. Each layer has a thickness of 0.13 mm, and is made of a Hercules type AS/3501-6 Graphite/Epoxy composite with the following properties: E "1.28;10 Pa, E "0.11;10 Pa, k "4.48;10 Pa VY WY VYWY k "k "1.53;10 Pa, o"1.5;10 kg/m, l "0.25 . WYXY XYVY WYVY
480
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486 Table 2 Natural frequency parameters for eight-layer composite square plate Mode
1
2
3
4
5
Conforming-Q Non-conforming-Q Nonconforming-T Analytical Solution [19] Published FEM [20] Experimental [20]
1.812 1.807 1.820 1.813 1.792 1.692
6.521 6.512 6.506 6.553 6.443 6.089
10.44 10.70 10.90 10.48 10.38 10.20
17.11 17.00 16.98 17.29 17.11 15.07
21.23 21.31 21.40 21.49 21.26 19.17
The fibre angles h, in the eight layers are 45, !45, !45, 45, 45, !45, !45, 45°. The non-dimensional frequency parameter for the composite-layered case, is defined as follows:E h VY where D " 12(1!l l ) VYWY WYVY and h represents the total thickness of the plate. The non-dimensional frequency parameters obtained by means of the four-node conforming quadrilateral element, four-node non-conforming quadrilateral element, and three-node non-conforming triangular element, developed in this work, were compared with the analytical solution published by Narita and Leissa [19], and also compared with an experimental work published by Crawley [20], and his finite element results which were based on a moderately thick quadrilateral shallow shell element developed by Lee and Pian [21], as shown in Table 2. The package results have generally proved to be closer to the analytical solution than those based on Lee and Pian’s element. However, some deviation (+8%) between all theoretical results and corresponding experimental results was observed. Crawley [20] mentioned that the experimental results represented an average of the measured natural frequencies of nominally identical samples. Since there was a slight variation in the finished thickness from plate to plate, the measured frequencies were linearly corrected to a reference thickness before averaging. No other correction such as for fibre volume fraction, has been applied. '"ua (oh/D
4.3. Composite layered rectangular plate case This case represents a rectangular plate with aspect ratio a/b"2 : 1, where a is the length in the x-direction and b is the width in the y-direction. The plate has 8 layers of composite material with properties, fibre angles and thickness similar to the previous case but with dimensions: a"162.4 mm, b"76.2 mm. The results obtained were compared with the published experimental and finite element results of Crawley [20], as shown in Table 3. All finite element results of the package were very close to each other, but with a consistent deviation of about 8% from the experimental results. 4.4. Composite layered cylindrical shell cases It is useful to validate the package for the analysis of a curved shell, and a cylindrical shell with published results has been chosen to achieve that. The geometry and boundary conditions of the
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
481
Table 3 Natural frequencies for eight-layer composite rectangular plate Mode
Experimental results [20]
Published FEM [20]
Conf. Q.
Non-conf. Q.
Non-conf. T.
1 2 3 4 5
31.3 185.8 214.0 533.0 653.0
31.9 191.3 228.2 565.3 708.3
31.85 191.16 227.36 562.22 703.78
31.86 192.96 227.20 571.35 703.84
29.50 180.85 224.81 540.12 696.66
Fig. 3. Composite layered cylindrical shell.
cylindrical curved shell are shown in Fig. 3, and the material properties are the same as in the previous two cases. Two different shells were attempted: the first has the fibres lie in angles 0, 0, 30, !30, !30, 30°, 0, 0 and the angles of fibres in the second shell are 0, 45, !45, 90, 90, !45, 45°, 0. The results were compared with published experimental results of Crawley [20] and his finite element results based on the element of Lee and Pian [21] as shown in Tables 4 and 5. It is clear from those tables, that acceptable results have been obtained by using the developed elements. However, the results of the conforming element are more accurate than those of non-conforming elements. This case study represents a shallow shell, where the lack of stiffness values corresponding to h in the equations of the non-conforming elements will lead to ill-conditioning of elemental X rotated matrices and can affect the accuracy of the results.
482
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486 Table 4 Natural frequencies for eight-layer composite cylindrical first shell section h of layers (0, 0, 30, !30, !30, 30, 0, 0°) Mode
Experimental results [20]
Published FEM [20]
Conf. Q.
Non-conf. Q.
Non-conf. T.
1 2 3 4 5
161.0 245.1 555.6 670.0 794.0
165.7 289.6 597.1 718.5 833.3
168.72 295.09 606.06 713.01 816.94
126.34 141.74 452.03 793.62 964.10
141.27 195.37 539.76 982.72 1156.19
Table 5 Natural frequencies for eight-layer composite cylindrical second shell section h of layers (0, 45, !45, 90, 90, 45, !45, 0°) Mode
Experimental results [20]
Published FEM [20]
Conf. Q.
Non-conf. Q.
Non-conf. T.
1 2 3 4 5
177.0 201.8 645.0 754.0 884.8
192.4 236.1 705.8 808.2 980.6
196.53 246.74 720.91 813.96 985.06
190.24 244.77 674.14 1091.27 1376.14
210.59 315.18 773.32 1196.09 1441.25
4.5. Case with variation of fibre angles This case represents a rectangular plate with the length along the x-direction being a"152.4 mm and the length in the y-direction b"76.2 mm. Two plates were considered one with five layers, and the thickness of each layer is 0.13 mm, where the angles of fibres are arranged as h, !h, h, !h, h, and the second is a single layered plate with a thickness equal to the total thickness of the five-layer plate, and the fibres at an angle h. To demonstrate the effect of the fibre angle, the natural frequency analysis has been carried out on plates with different values of h, i.e. h"0, 15, 30, 45, 60, 75, 90°. The first two non-dimensional frequencies as obtained by finite element analysis with the conforming four-node quadrilateral, non-conforming four-node quadrilateral, and non-conforming three-node triangular elements, together with the corresponding analytical solution of Narita and Leissa [19], were plotted against fibre angles, for the five-layer as well as the single- layer plates, as shown in Figs. 4 and 5. Numerical results for the five-layer case showing the first five modes of vibration are demonstrated in Table 6. The following points can be deduced from the given results: (i) Results obtained using different finite elements for the two plates at different fibre angles are in a very close agreement to those obtained by the corresponding analytical solutions.
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
483
Fig. 4. Effect of fibre angle on the first natural frequency for single- and five-layer plates.
Fig. 5. Effect of fibre angle on the second natural frequency for single- and five-layer plates.
(ii) Increasing the fibre angle h, changes the natural frequency, both for the single layer and the five-layer plates, at the natural modes displayed. (iii) The change in the fibre angle h have a greater effect on the results of the five-layer plate than on those of the single layer plate. This may indicate the advantage of using more layers if the fibre angle effect is to be fully exploited. To investigate the thickness effect on the natural frequencies of a multi-layered plate, the five-layer plate was employed with two values of plate thickness being considered, the first is the same as before, but the second thickness was ten times greater. The first four non-dimensional natural
484
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Table 6 Non-dimensional natural frequencies for five-layer rectangular plate at different fibre angle h Fibres angle (deg) h"0
h"15
h"30
h"45
h"60
h"75
h"90
Mode 1
Conf-Q Nonconf-Q Nonconf-T Analytical
3.511 3.499 3.501 3.513
3.152 3.137 3.145 3.153
2.445 2.423 2.443 2.442
1.703 1.690 1.701 1.700
1.176 1.172 1.177 1.176
0.942 0.942 0.945 0.946
0.894 0.894 0.896 0.897
Mode 2
Conf-Q Nonconf-Q Nonconf-T Analytical
7.057 7.061 7.056 7.068
8.452 8.433 8.455 8.474
10.260 10.230 10.260 10.310
10.070 10.120 10.060 10.080
7.263 7.345 7.367 7.259
5.889 5.905 5.971 5.888
5.388 5.389 5.432 5.390
Mode 3
Conf-Q Nonconf-Q Nonconf-T Analytical
21.910 21.270 21.590 21.970
19.830 20.040 19.680 19.910
15.200 15.840 15.270 15.260
11.130 11.430 11.270 11.190
9.431 9.485 9.485 9.476
6.921 6.938 6.990 6.931
5.601 5.606 5.701 5.602
Mode 4
Conf-Q Nonconf-Q Nonconf-T Analytical
25.230 25.150 26.260 25.840
27.110 26.570 27.480 27.490
31.410 30.920 31.220 31.470
28.980 29.290 28.960 28.980
20.6 21.05 21.2 20.62
16.5 16.6 17.1 16.5
15.7 15.70 16.3 15.7
Mode 5
Conf-Q Nonconf-Q Nonconf-T Analytical
27.760 26.030 27.360 27.840
30.850 30.010 31.27 31.340
36.520 36.190 36.850 37.120
33.870 34.610 34.020 33.950
29.010 29.260 29.190 29.090
21.450 21.540 21.790 21.490
16.930 16.940 17.300 16.950
Fig. 6. Effect of thickness on the first natural frequency for five-layer plate.
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
485
Fig. 7. Effect of thickness on the second natural frequency for five-layer plate.
frequencies for the thin and thick plates were plotted against fibre angles h, as shown in Figs. 6 and 7. It is clear from those figures, that the transverse shear effect tends to reduce the non-dimensional frequency, but this reduction is greater at h values between 15 and 60°, approximately. It is also clear that the effect of the fibre angle on natural frequency may be reduced by increasing the plate thickness.
5. Conclusions It is clear from the previous case studies that conforming and non-conforming elements developed in this work have led to accurate and reliable results. Although, the difference between their corresponding results was insignificant for flat plate cases, the conforming elements have shown more accurate results for the curved shell case. The fact that the conforming element has stiffness and mass contributions corresponding to h , h , and h will make them well conditioned, V W X and provide more accurate results for shallow shell cases than those obtained from non- conforming elements. The conforming elements are also more flexible in use with curved shells and box structures, because they require boundary conditions to be in terms of slope angles with respect to global axes, whilst non-conforming elements require them in terms of local axes, which are not easy for the user to know. Although the developed elements consider transverse shear effects, they do not suffer from shear locking as do Mindlin-type elements, and reduced integration techniques are not required for the elements developed in this work. Transverse shear effects, which are increased by increasing thickness, tend to reduce the non-dimensional frequencies. The fibre angle h can have a significant effect on the stiffness matrix of composite layered plates and shells, and this will result in a corresponding effect on their natural frequencies. This effect can be increased by increasing the number of layers, while keeping the same thickness. The h effect may
486
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
provide a tool for the designers to tune the natural frequencies away from resonance due to excitation frequency, by simply changing the fibre angles without having to modify the mechanical design.
References [1] Epstein M, Huttelmaier HP. A finite element formulation for multilayered and thick plates. Computers and Structures 1983;16:645. [2] Widera, OEG, Moumene, M. Cutout in laminated plates and shells. Advances in Macro-Mechanics. Composite Material Vessels Components 1984;146:155. [3] Has, TM, Wang JTS. A theory of laminated cylindrical shells consisting of orthotropic laminae. AIAA Journal 1970;8:2141. [4] Witt M, Sobczyk K. Dynamic response of laminated plates to random loading. International Journal of Solids Structures 1980;16:231. [5] Pillasch DW, Majerus JN, Zak AR. Dynamic finite element model for laminated structures. Computers and Structures 1983;16:449. [6] Nelson RB, Lorch DRA. Refined theory for laminated orthotropic plates. Journal of Applied Mechanics 1974;41:177. [7] Reddy JN. A refined nonlinear theory of plates with transverse shear deformation. International Journal of Solids and Structures 1984;20:881. [8] Reddy, JN. A simple higher order theory for laminated composite plates. Journal of Applied Mechanics 1985;51:745. [9] Wang FY. Two dimensional theories deduced from three-dimension theory for a transversely isotropic body—I. Plate problems. International Journal of Solids and Structures 1990;26:445. [10] Chang B, Chen DC, Shabana AA. Effect of the coupling between stretching and bending in the large displacement analysis of plates. International Journal of Numerical Methods in Engineering 1990;30:1233. [11] Wung PM, Reddy JN. A transverse deformation theory of laminated composite plates. Computers and Structures 1991;41:821. [12] Tessler A. An improved plate theory of +1,2,-order for thick composite laminates. International Journal of Solids and Structures 1993;30:981. [13] Cho M, Parmerter R. Finite element for composite plate bending based on efficient higher order theory. AIAA Journal 1994;11:2241. [14] Crisfield MA. Non-linear finite element analysis of solids and structures. New York: Wiley, 1991. [15] El-Zafrany AM, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for triangular elements. International Journal of Numerical Methods in Engineering 1986;23:275. [16] El-Zafrany AM, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for quadrilateral elements. International Journal of Numerical Methods in Engineering 1986;23:1939. [17] Zienkiewicz OC. The finite element method, 3rd ed. New York: McGraw-Hill, 1977. [18] El-Zafrany AM, Fadhil S, Debbih M. An efficient approach for boundary element bending analysis of thin and thick plates. International Journal of Computers and Structures 1995;56:565. [19] Narita Y, Leissa AW. Frequencies and mode shapes of cantilevered laminated composite plates. Journal of Sound and Vibration 1992;54:161. [20] Crawley EF. The natural modes of graphite/epoxy cantilever plates and shells. Journal of Composite Materials 1979;3:195. [21] Lee SW, Pian TH. Improvements of plate and shell finite elements by mixed formulations. AIAA Journal 1978:16.
A high-order shear element for nonlinear vibration analysis of composite layered plates and shells O. Attia *, A. El-Zafrany Imperial College of Science, Technology and Medicine, Mechanical Engineering Department, London SW7 2BX, UK School of Mechanical Engineering, Cranfield University, Cranfield, Bedford MK43 0AL, UK Received 30 June 1998
Abstract This paper introduces a family of high-order facetted shell elements for linear and nonlinear stress and vibration analysis of composite layered plate and shell structures. Engineering slope angles are employed in element equations, and transverse stresses are expanded over the thickness. The lateral deflection is modelled by conforming or non-conforming Hermitian shape functions, within rectangular or paralellogrammic and triangular elements. Nonlinear terms associated with geometrical nonlinearity are also derived using a practical approach based upon the actual components of strain. A finite element programming package was designed employing the newly developed elements. Several case studies have been investigated and package results were compared with existing theoretical and/or experimental results. It has been proved that the developed elements can lead to accurate estimations of natural frequencies. The effect of fibre angles on natural frequencies has also been investigated with some case studies, and the results proved that the package can be a useful tool for the design optimization of composite layered plate and shell structures. 1998 Elsevier Science Ltd. All rights reserved. Keywords: High-order shear element; Nonlinear vibration analysis; Composite layered plates and shell
Notation a, b A B D D
side lengths of a rectangular plate matrix of displacement derivatives strain-displacement matrix plate flexural rigidity stress-strain matrix
* Corresponding author. 0020-7403/99/$ — see front matter 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 4 0 3 ( 9 8 ) 0 0 0 7 7 - 0
462
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
E F F,G,H G G G h I K l, m, n M N G u, v, w x, y, z x , y , z c d e h ,h V W q '
, t u k l o r
Young’s modulus nodal loading vector Hermitian shape functions plate thickness unit matrix stiffness matrix directional cosines mass matrix Lagrangian shape functions displacement components in x, y, z directions, respectively Cartesian coordinates Coordinates based on material principal axes (see Fig. 1) transverse shear strain vector nodal displacement vector strain vector average slope angles in the x and y directions, respectively transverse shear stress vector non-dimensional frequency parameter (,c , c ) transverse shear strains at the plate midplane VX WX natural frequency shear modulus Poisson’s ratio density stress vector
1. Introduction The use of composite materials in automotive and aircraft structures has increased in recent years, due to the high strength to weight ratio of composites and their many other advantages. Fibrous composites can optimally be designed with the tuning of fibre angles in different layers so as to obtain the required properties. Nevertheless, the sensitivity of composite layered plates to shear stresses necessitates the need to estimate transverse shear stresses as accurately as possible. Many industrial applications may also have a geometrical nonlinearity in the elastic range, which should be considered in the analysis. There are a number of finite element publications for composite materials reported in the literature, most of which are limited to static analysis. Epstein and Huttelmaier [1] presented a finite element formulation for the analysis of multilayered thick plates, where the transverse shear and normal strains were considered. Widera and Moumene [2] used an eight-node isoparametric quadrilateral shell element, which included the effect of transverse shear strain in the analysis of laminated plates and shells. In the area of dynamic analysis of laminated structures, Has and Wang [3] presented an analytical approach to the analysis of laminated cylindrical shells. They did not employ the Kirchhoff hypothesis of non-deformable normals by including the inertia terms in the equilibrium equations. Witt and Sobczyk [4] discussed the dynamic response of laminated plates to random
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
463
loading. Pillasch et al. [5] analysed thick laminated plates with a finite element model based on a three-dimensional element rather than plate theory. The assumption that displacements are linear functions of z, the coordinate in the thickness direction, has proved to be inadequate for predicting gross laminate response. Therefore, many higher-order theories incorporating transverse shear strains, have been proposed for better accuracy. Nelson and Lorch [6] used up to second-order terms of z in their displacement field. Reddy [7, 8] introduced a refined nonlinear theory of plates with transverse shear deformation, and a simple higher-order theory for laminated composite plates. The deduction of two-dimensional theories, from the three-dimensional theories, in a transversely isotropic body was presented by Wang [9]. Chang et al. [10] presented a Lagrangian formulation for large translation and rotation analysis of deformable rectangular homogeneous plates. They included the coupling effects between bending and stretching, in the finite element analysis of homogeneous plates. Wung and Reddy [11] presented a first-order shear/fourth-order transverse deformation theory for laminated composite plates. A recent theory was introduced by Tessler [12], who derived a two-dimensional laminated theory for linear elastic analysis of thick composite plates with the equivalent single-layer assumptions for the displacements, transverse shear strain, and transverse normal stress. Cho and Parmerter [13] developed a three-node non-conforming triangular bending element, with five degrees of freedom, for symmetric laminated composites. The derivation of this element is based on a higher order plate theory evolved earlier by the same authors. Many of the high-order element derivations presented in the literature are based on hypothetical nodal parameters, which are difficult to handle with different types of boundary conditions. This paper introduces a family of high-order facetted shell elements for linear and nonlinear stress and vibration analysis of composite layered plate and shell structures. Engineering slope angles are employed in element equations, and parabolic distributions of transverse shear strains over the thickness are considered. The lateral deflection is modelled by conforming or nonconforming Hermitian shape functions, within rectangular or parallelogrammic and triangular elements. The degrees-of-freedom have a maximum total of nine per node, but they are automatically reduced for coplanar nodes, so as to make it easy and efficient to use the elements for the analysis of box structures and curved shells. Non-linear terms associated with geometrical nonlinearity are also derived using a practical approach based upon the actual components of strain.
2. Deformation theory 2.1. Transverse strain modelling Consider a composite layered plate, at an instant of time t consisting of a number of orthotropic layers N , as shown in Fig. 1. Let the midplane of the plate be the Cartesian x-y plane, and the total J thickness at any point (x, y) on the midplane is h (x, y). Hence, the equations of the lower and upper surfaces of the plate are: z "!h/2 and z "h/2, respectively. If the values of c and c represent VX WX
464
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Fig. 1. Composite layered plate.
the transverse shear strains at the midplane, then from the boundary conditions the following values of transverse shear strains are defined as follows: at z "!h/2, and at z "h/2: c "c "0 VX WX at z "0: c "c , c "c . VX VX WX WX The composite layers are considered to be firmly bonded together, hence the distributions of displacement components represent continuous functions with respect to the z-direction. Considering also the infinitesimal shear strain—displacement equations: *u *w c " # VX *z *x
(1)
*v *w c " # WX *z *y
(2)
the transverse shear strains can be approximated as continuous functions with respect to z. Hence, the distributions of the transverse shear strains across the plate thickness can be deduced by applying Lagrangian interpolation to their given values, i.e. c (x, y, z, t),c (x, y, t) [1!4z/h] VX VX c (x, y, z, t)"c (x, y, t) [1!4z/h] WX WX
(3) (4)
2.2. Displacement and velocity components Substituting from Eqs. (3) and (4) into Eqs. (1) and (2) and integrating over the thickness, assuming that the lateral deflection w can be approximated to be independent of z, then
4z *w u(x, y, z, t)"u (x, y, t)!z #c (x, y, t) z! VX 3h *x *w 4z v (x, y, z, t)"v (x, y, t)!z #c (x, y, t) z! WX *y 3h
where u and v represent the values of u and v at z"0.
(5) (6)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
465
The corresponding velocity components can be obtained by differentiating displacement equations with respect to time, i.e.
*wR 4z uR (x, y, z, t)"uR !z #cR z! VX *x 3h *wR 4z vR (x, y, z, t)"vR !z #cR z! WX *y 3h
(7) (8)
2.3. Strain components The x-y strain components will be considered finite and by using a Lagrangian frame of reference, they can be expressed in terms of displacement components [14]: *u 1 e" # V *x 2 *v 1 e" # W *y 2
*u *v *w # # *x *x *x
*u *v *w # # *y *y *y
*u *v *u *u *v *v *w *w c " # # # # VW *y *x *x *y *x *y *x *y
(9) (10) (11)
Defining the following strain vector: e "+e e c , VW V W VW then it can be deduced from Eqs. (9)—(11) that
(12)
e "e #e VW * where *u *x *v e " *y *u *v # *y *x
(13)
(14)
which represents the infinitesimal terms, and
1 *u 1 *v 1 *w # # 2 *x 2 *x 2 *x 1 *u 1 *v 1 *w # # e " * 2 *y 2 *y 2 *y *u *u *v *v *w *w # # *x *y *x *y *x *y
which represents the effect of finite strains.
e* V , e* W e* VW
(15)
466
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Using Eqs. (5) and (6), the infinitesimal strain vector can be expressed as follows: e (x, y, z, t)"e (x, y, t)!ze (x, y, t)#f (z) e (x, y, t) K @ Q
(16)
where *u *x *v e " K *y *u *v # *y *x
*w *x *w , e" @ *y *w 2 *x *y
*c VX *x *c WX , e" Q *y *c *c VX# WX *y *x
(17)
and f (z)"z!4z/(3h) . Substituting from Eqs. (5) and (6) into Eq. (15) it can be proved that: 1 e*" V 2
*u *v *w # # #z *x *x *x
*w *w # *x *x *y
#f (z)
*c *c *u *w *v *w VX # WX !2z # *x *x *x *x *x *x *y
!2zf (z)
*v *c *c *w *c *w *u *c VX VX# WX # WX #2 f (z) *x *x*y *x *x *x *x *x *x
1 e*" W 2
*u *v *w # # *y *y *y
*w *w #f (z) # *y*x *y
#z
*c *c VX # WX *y *y
*u *w *u *c *v *w *v *c VX# WX # #2 f (z) *y *y*x *y *y *y *y *y *y
!2z
!2z f (z)
*c *w *c *w VX # WX *y *y*x *y *y
.
*u *u *v *v *w *w c* " # # VW *x *y *x *y *x *y !z
(18a)
*u *w *u *w *v *w *v *w # # # *x *y*x *y *x *x *y *y *x *y
(18b)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
!z f (z) #f (z)
#z
*c *w *c *w *c *w *c *w VX # VX # WX # WX *x *x*y *y *x *x *y *y *x*y
*u *c *v *c *v *c *u *c VX# VX# WX# VX *y *x *x *y *y *x *x *y
467
*w *w *c *c *w *w *c *c VX VX# WX WX # #f (z) *x *x *y *x*y *y *x *y *x *y
(18c)
3. Finite element derivation 3.1. Displacement and velocity interpolation The composite-layered plate element is defined in terms of the midplane (z"0), where the nodal parameters are also defined. The geometry of the plate midplane is defined in terms of nodal coordinates with respect to a Cartesian x—y—z system of coordinates. Most of the derivatives described here are with respect to a local system of axes, where the x—y plane is in the midplane, and the z-axis is normal to that midplane. Physically, the element is constructed of a number of layers (N ), similar to that shown in Fig. 1. J The l layer’s lower and upper surfaces are described by the following equations: z"ZJ (x, y), z"ZJ (x, y) * 3 and the layer thickness is hJ(x, y)"ZJ! ZJ . 3 * The layer has its own material axes, and the stress—strain relationships are defined with respect to the element local x, y, z axes, as follows: rJ V sJ rJ , rJ "D(l) e , sJ, VX "lJ VW VW W sJ WX sJ VW
c VW . c WX
(19)
The matrices D(l ) and l(l ) are stress—strain matrices rotated from material axes to element axes. The midplane displacement components in the x and y directions will lead to continuous displacement components, if they are represented by at least C-continuous functions. Hence, Lagrangian interpolation is suitable for their modelling, and for an n-node finite element: L L u (x, y, t)" u (t) N (x, y), v (x, y, t)" v (t) N (x, y) G G G G G G
(20)
where u , v are the values of the x—y displacement components, at the midplane node i, and time t, G G and N (x, y) represents a Lagrangian shape function. G
468
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Similarly, the transverse shear components, at the midplane can be interpolated in terms of Lagrangian shape functions as follows: L c (x, y, t), (x, y, t)" (t) N (x, y) , VX G G G L c (x, y, t),t (x, y, t)" t (t) N (x, y) WX G G G where
(21a) (21b)
(t),c (x , y , t), t (t),c (x , y , t) (22) G VX G G G WX G G The lateral deflection w must be C continuous over the plate, and it should be interpolated by means of a Hermitian interpolation. Two types of interpolation, non-conforming and conforming, based upon El-Zafrany and Cookson [15, 16] are employed, where the lateral deflection is interpolated as follows: L w(x, y, t)" +w (t)F (x, y)#w (t) G (x, y)#w (t) H (x, y)#w (t) Q (x, y), , (23) G G G V G G W G G VW G G where w , w , w , and w represent the values of w, *w/*x, *w/*y and *w/*x*y at node i and G G V G W G VW time t, and F , G , H and Q are Hermitian shape functions. Notice also that Q does not exist for G G G G G the non-conforming Hermitian interpolation. From the previous analysis it is clear that the displacement parameters required at each node i are u , v , , t , w , w , w , w , and a practical representation of nodal values will be described G G G G G G V G W G VW in Section 3.8. To simplify the derivation, the local nodal displacement vector for an n-node element will be defined at time t as follows: d(t)"+d (t) d (t) d (t), K @ Q where
(24)
d (t),+u (t) v (t) u (t) v (t) 2 u (t) v (t), , (25) K L L d (t),+w (t) w (t) w (t) w (t) 2 w (t) w (t) w (t) w (t), , (26) @ V W VW L L V L W L VW d (t),+ (t) t (t) 2 (t) t (t), (27) Q L L The velocity components at any point (x, y, z) inside the plate at time t can be represented vectorially in terms of nodal values and shape functions as follows:
q (x, y, z, t),
q (x, y, z, t) VW q (x, y, z, t) X
(28)
where
uR "N (x, y) d (t)!zN (x, y) d (t)#f (z) N (x, y) d (t) q , * K @ @ * Q VW vR
(29)
q ,[wR ]"N (x, y) d (t) X & @
(30)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
469
and
(31)
,
(32)
N 0 N 0 2 N 0 L N (x, y)" * 0 N 0 N 2 0 N L *F *G *H G G G *x *x *x
*Q G 2 *x
*F *G *H G G G 2 *y *y *y
*Q G 2 *y
2 N (x, y)" @
N (x, y)"[2 F G H Q 2 ] . & G G G G
(33)
3.2. Infinitesimal strain components Substituting from Eqs. (20) into Eq. (17), it can be deduced that:e (x, y, t)"B (x, y) d (t) K K K where *N G 2 0 *x B (x, y)" 2 K 2
0
(34)
2
*N G 2 *y
(35)
*N *N G 2 G *x *y
Similarly it can be proved that:e (x, y, t)"B (x, y) d (t) , Q K Q
(36)
and e (x, y, t)"B (x, y) d (t) @ @ @ where
(37)
2
*F G *x
*G G *x
*H G *x
*Q G *x
2
B (x, y)" 2 @
*F G *y
*G G *y
*H G *y
*Q G *y
2
2 2
*F *G *H *Q G 2 G 2 G 2 G 2 *x *y *x *y *x *y *x *y
(38)
470
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Substituting from Eqs. (34), (36), (37) into Eq. (16) then it can be shown that e (x, y, z, t)"B (x, y) d (t)!zB (x, y) d (t)#f (z) B (x, y) d (t) . K K @ @ K Q
(39)
Note also that the transverse shear strain components, as defined by Eqs. (3) and (4), can be represented vectorially as follows:
c (x, y, z, t),
c c (x, y, t) VX "f (z) VX c c (x, y, t) WX WX
(40)
where
4z f (z)" 1! h
(41)
Substituting from Eq. (21) into Eq. (40), it can be deduced that c (x, y, z, t)"f (z) B (x, y) d (t) , Q Q
(42)
where B (x, y),N (x, y) Q *
(43)
3.3. Modelling of finite strain contributions Defining the following rotation vectors:
*u *v *u *u h " K *x *x *y *y * *t * *t h" Q *x *x *y *y
*w *w *w *w h" @ *x *x*y *y*x *y *w *w h " U *x *y
(44)
(45)
(46)
(47)
then it can be proved from Eq. (15) that e "+A h #A h #zA h #f (z) A h * K K U U @ @ Q Q !z f (z) [A h #A h ]!z [A h #A h ]#f (z)[A h #A h ],, @ Q Q @ @ K K @ K Q Q K
(48)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
where
*u *v *x *x
0
*u *v A " 0 0 K *y *y *u *v *u *v *y *y *x *x
A" @
*w *x
*w *x*y
0
0
*w *y *x
*w *y
* *t *x *x
0
0
0
471
0
* *t , A" 0 0 Q *y *y * *t * *t *y *y *x *x *w *x
0
*w *w *y *x *y *w *w *x *x *y
0
*w , A " 0 U *y *w *w *y *x
Using the interpolation equations for displacement components, the following differential expressions can be obtained: dh "G dd , dh "G dd , dh "G dd , dh "G dd K K K Q K Q @ @ @ U U @ where *N G 2 0 2 *x *N G 2 2 0 *x G " , K *N G 2 0 2 *y *N G 2 2 0 *y
G" @
*G G *x
*H G *x
*Q G *x
2
*F G *x
2
*G *H *Q *F G G G 2 G *x *y *x *y *x *y *x *y
G ,N . U @
*F G *y
*G G *y
*H G *y
*Q G *y
(50)
2
*F *G *H *Q G G G G 2 2 *y *x *y *x *y *x *y *x 2
(49)
(51)
2 (52)
472
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Due to the non-linearity of terms in Eq. (48), it will be represented in the derivations as an integration of its differential, and it can be deduced from Eqs. (48) and (49) that: de "A G dd #A G dd #f (z) A G dd #z A G dd * K K K U U U Q K Q @ @ @ !z f (z) [A G dd #A G dd ]!z [A G dd #A G dd ] @ K Q Q @ @ @ K K K @ @ #f (z) [A G dd #A G dd ] K K Q Q K K
(53)
3.4. Element stiffness matrix for infinitesimal strains For the case of infinitesimal strains, it can be deduced from Eqs. (19), (39) and (40) that the corresponding stress vectors at the lth layer are as follows:rJ "DJ [B d !zB d #f (z) B d ] , VW VW K K @ @ K Q
(54)
sJ"lJ [ f (z) B d ] . Q Q
(55)
The strain energy density or the strain energy per unit volume at any point inside the element can be expressed as follows: ºM (x, y, z, t)" e r,e rxy#c s .
(56)
Hence, for the case of infinitesimal strains, the strain energy density at a point inside the lth layer will be: ºM J" [d B !zd B #f (z) d B ] DJ [B d !zB d #f (z) B d ] @ @ Q K VW K K @ @ K Q K K # [ d B f (z) lJ B d ] . Q Q Q Q
(57)
Notice that all the B matrices are functions of (x, y), and only the functions of z and the D matrices, which may be different for different layers, need to be integrated over the thickness. Writing the integration of every function of z over the total thickness in terms of integrations over different layers, then the following expressions of modified D matrices, at any point (x, y) on the midplane, can be deduced from the integration of Eq. (57) with respect to z:
D (x, y)" KK D (x, y)" @@ l (x, y)" QQ
F
,J DJ dz" [ZJ!ZJ ] DJ , VW VW 3 * \F J
F
,J (ZJ)!(ZJ) 3 * zDJ dz" DJ VW VW 3 \F J
F
\F
f (z) lJ dz ,
(58)
(59)
(60)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
D (x, y)"D (x, y)" K@ @K
F
\F F
zDJ dz , V W
473
(61)
f (z) DJ dz, (62) VW \F F D (x, y)"D (x, y)" z f (z) DJ dz . (63) @Q Q@ VW \F Integrating Eq. (57) with respect to the volume of the element, then the strain energy of the element can be expressed at any time t as follows: D (x, y)"D (x, y)" KQ QK
º(t)" + d K d #d K d !d K d K KK K @ @@ @ K K@ @ #d [KV W#KI] d !dR KR d #dR K d Q QQ QQ Q @ K@ K K KQ Q #d K d !d K d !d K d , Q KQ K @ @Q Q Q @Q @ where
(64)
CJCKCLR BK DKK BK dx dy, K@@"CJCKCLR B@ D@@ B@ dx dy, KV W" QQ CJCKCLR BK DQQ BK dx dy, KIQQ "CJCKCLR BQ lQQ BQ dx dy, K " B D B dx dy, K " B D B dx dy , K@ K K@ @ KQ K KQ K CJCKCLR CJCKCLR K " B D B dx dy . @Q @ @Q K CJCKCLR K " KK
The integrations over elements are carried out numerically in a way similar to that employed in two-dimensional elements (see for e.g. Ref. [17]), i.e.
\ \
K " KK
(x, y) B D B J dm dg , K KK K (m, g)
etc., and a double Gaussian quadrature is employed for the numerical evaluation. Notice that Eq. (64) may also be rewritten as follows: º(t)" d (t) K d (t), (65) where d (t) is the nodal displacement vector of the element at time t, as defined by Eq. (24), and K is the element stiffness matrix for the case of infinitesimal strains and can be constructed as follows: K
KK K" !K K@ K KQ
!K K
K@
@@ !K @Q
K
KQ . !K @Q K V W#K I QQ QQ
(66)
474
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
3.5. Element stiffness matrix for finite strain contributions The virtual or differential value of the finite strain contributions, due to a virtual or differential displacement field is de (x, y, z, t)"+de* de* dc* , * V W VW assuming negligible finite strain contributions for the transverse shear strains. The corresponding increment of strain energy density at any point (x, y, z) in layer l is dºM J (x, y, z, t)"de rJ * VW * Substituting from Eq. (53) into Eq. (67) then:
(67)
dºM "+dd G A #dd G A #z dd G A #f (z) dd G A * K K K @ U U @ @ @ Q K Q !z [dd G A #dd G A ]#f (z) [dd GR A #dd G A ] K K @ @ @ K K K Q Q K K !z f (z) [dd G A #dd G A ], rJ @ @ Q Q @ @ VW By direct multiplication of matrices, the following equations can be proved:
(68)
A r ,Sh "SG d , A r ,Sh "SG d , K VW K K K Q VW Q K Q A r ,Sh "SG d , A r ,Sw h "Sw G d , @ VW @ @ @ U VW U U @ where p
0
q
0 VW p 0 q V VW , S " pV qVW S" U q p q 0 p 0 VW W VW W 0 q 0 p VW W and Eq. (68) can be rewritten as follows: V 0
dºM J"dd G SJG d #dd G SJG d #zdd G SJG d #f (z)dd G SJG d K K K K @ U U U @ @ @ @ @ Q K K Q * !z[dd G SJG d #dd G SJG d ]#f (z) [dd G SJG d #dd G SJG d ] K K @ @ @ @ K K K K K Q Q K K K !zf (z) [dd G SJG d #dd G SJG d ] (69) @ @ K Q Q K @ @ Notice that the G matrices are independent of z, but S and S may vary with z even within the same U layer. To simplify the derivations, the average value of the stress vector inside each layer will be employed as a constant parameter within the thickness of the layer, i.e. SJ+ [S(ZJ )#S(ZJ )] * 3 SJ + [Sw (ZJ)#Sw (ZJ )] U * 3
(70) (71)
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
475
The following integrated z functions can also be defined:
F
,J S dz" [ZJ!ZJ] SJ , 3 * \F J
S " KK
F
,J S dz" [ZJ!ZJ ] SJ , U U 3 * \F J
S " UU
F
S " @@ S " QQ
zS dz ,
(72)
(73)
(74)
\F F
\F
S "S " K@ @K
(75)
zS dz ,
(76)
4z z! S dz, 3h F
\F
S "S (x, y)" KQ QK S (x, y)"S " @Q Q@
F
(77)
4z S dz. z z! 3h
(78)
\F
F
\F
4z z! S dz, 3h
The corresponding increment of the strain energy is obtained by integrating Eq. (69) over the volume of the element, i.e. dº (t)"dd KN d #dd (KN #KN )d * K KK K @ @@ UU @ #dd KN d !dd KN d !dd (KN )d Q QQ Q K K@ K @ K@ K #dd KN d #dd (KN )d !dd KN d !dd (KN )d K KQ Q Q KQ K @ @Q Q Q @Q @ where
CJCKCLR GK SKK GK dx dy, KN@@"CJCKCLR G@ S@@ G@ dx dy KN " G S G dx dy, KN " G S G dx dy UU U UU U QQ K QQ K CJCKCLR CJCKCLR G S G dx dy, KN " G S G dx dy KN " K K@ @ KQ K KQ K K@ CJCKCLR CJCKCLR G S G dx dy KN " @ @Q K @Q CJCKCLR KN " KK
(79)
476
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Eq. (79) can also be rewritten as follows: dº (t)"dd (t) KNd (t) * where
(80)
KN !KN KN KK K@ KQ (81) KN" !(KN )R KN #KN !KN K@ @@ UU @Q (KN )R !(KN )R KN KQ @Q QQ and the total increment of the strain energy due to total elastic stresses and strains will be dº(t)"dd (t) [K#K p ] d (t)
(82)
3.6. Element mass matrix The kinetic energy per unit volume, at any point inside the element is: KE" o [uR #vR #wR ], o [q q #q q ] VW VW X X Substituting from Eqs. (29) and (30) into Eq. (83), it can be shown that:
(83)
KE" o +[d N !zd N #f (z)d N ] K * @ @ Q * ;[N d !zN d #f (z)N d ]#d N N d , (84) * K @ @ * Q @ & & @ The N matrices, in Eq. (84), are independent of z, and density functions can be defined at any point (x, y) on the element midplane, in a way similar to that employed in Eqs. (58)—(63), i.e. ,J o dz, [ZJ (x, y)!ZJ (x, y)] oJ 3 * \F J F ,J (ZJ (x, y))!(ZJ (x, y)) 3 * zodz" oJ , etc. o (x, y)" @@ 3 \F J Integrating Eq. (84) over the volume of the element, then the kinetic energy of the element at a time t can be expressed as follows:
o (x, y)" KK
F
KE (t)" +d M d #d (M #M ) d @ @@ && @ K KK K #d M d !d M d !d M d Q QQ Q K K@ @ @ K@ K #d M d #d M d !d M d !d M d , , K KQ Q Q KQ K @ @Q Q Q @Q @ where
CJCKCLR oKK N* N* dx dy,
M " KK
CJCKCLR o@@ N@ N@ dx dy,
M " @@
(85)
etc.
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
477
Eq. (85) can also be simplified as follows: KE (t)" d (t) Md (t) where
(86)
d (t)"+d (t) d (t) d (t), K @ Q and M is the mass matrix of the element which is constructed as follows: M !M M KK K@ KQ M #M !M M" !M K@ @@ && @Q M !M M KQ @Q QQ
(87)
3.7. Element dynamic equations Applying the principle of minimum total potential energy at an instant of time t, then ds (t)"dKE(t)#dº (t)!d¼(t)"0
(88)
where d¼ is the work done by the applied forces at time t, and can be expressed as follows: d¼(t)"dd (t) F (t)
(89)
where F (t) represents a nodal loading system equivalent to the actual applied force. Substituting from Eqs. (82), (86), (89) into (88), then it can be deduced that M d® (t)#[K#K p] d (t)"F (t)
(90)
which represents the dynamic matrix equation of the element. 3.8. Practical nodal values and rotated matrices The use of derivatives of w as nodal values causes some difficulties with the specification of boundary conditions for thick plates. Engineering slope angles can be defined so as to lead to displacement components linearized over thickness [18], and they are defined as follows: *w *w h " !0.8c , h "! #0.8c WX W VX V *y *x
(91)
*w h" for conforming elements, "0 for non-conforming elements. X *x*y The average values of transverse shear strains are used as nodal values, and can also be expressed as follows: cN " c ,u , cN "c ,!u , W WX WX V VX VX
(92)
478
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Defining a new nodal displacement vector based on such practical values as follows: d "+ 2 u v w (h ) (h ) (h ) y(u ) (u ) (u ) 2 , LCU G G G VG WG XG VG WG XG then it can be shown that
(93)
d"Pd
, (94) LCU where P is a transformation matrix defined from the previous relationships. The element local axes are defined with respect to structural global axes by means of directional cosines, and it can be shown that d "Rd , (95) LCU E where R is the element rotation matrix, and d is the global nodal displacement vector. Hence, it can E be deduced that: d"PRdg,Ldg
(96)
and the element matrices defined with respect to the global axes can be obtained as follows: K "LKL, E
KN"LK pL, M "LML E E
(97)
4. Case studies Based upon the theory described in previous sections, a modular finite element computer programming package was designed. The package was coded in FORTRAN 77, and tested on VAX/VMS main frame computers, Unix workstations, and DOS/WINDOWS operating PCs. The package consists of several modules, where each module is a file containing a group of subroutines. It can be employed for finite element static and dynamic analysis of plates and shells made of isotropic or composite layered materials. A number of case studies, with published analytical and/or experimental results, have been analysed so as to assess the accuracy and reliability of different aspects of the developed theory. 4.1. Dynamic analysis of isotropic square plate case This case represents a square plate, as shown in Fig. 2, made of isotropic material with the following properties: E"1.28;10 Pa, l"0.333, o"1.5;10 kg/m . The plate is fixed at one edge, whilst the other three edges remain free. The package results were compared with a published analytical solution carried out by Narita and Leissa [19], as demonstrated in Table 1, where non-dimensional frequency parameters were displayed. The nondimensional frequency parameter U is defined as follows: '"ua (oh/D
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
479
Fig. 2. Composite layered square plate.
Table 1 Natural frequency parameters for isotropic square plate Mode
1
2
3
4
5
Conforming Q. Non-conforming Q. Non-conforming T. Analytical Solution [19]
3.461 3.456 3.453 3.459
8.361 8.362 8.393 8.358
21.10 21.34 21.34 21.09
27.04 26.81 26.82 27.07
30.52 30.47 30.67 30.56
It is clear from Table 1 that natural frequency parameters, for different modes of vibration, as obtained by means of different finite elements, agree very well with corresponding analytical solutions. 4.2. Composite layered square plate case This case has overall dimensions similar to the previous case shown in Fig. 2, but it has eight layers of a fibrous composite material. Each layer has a thickness of 0.13 mm, and is made of a Hercules type AS/3501-6 Graphite/Epoxy composite with the following properties: E "1.28;10 Pa, E "0.11;10 Pa, k "4.48;10 Pa VY WY VYWY k "k "1.53;10 Pa, o"1.5;10 kg/m, l "0.25 . WYXY XYVY WYVY
480
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486 Table 2 Natural frequency parameters for eight-layer composite square plate Mode
1
2
3
4
5
Conforming-Q Non-conforming-Q Nonconforming-T Analytical Solution [19] Published FEM [20] Experimental [20]
1.812 1.807 1.820 1.813 1.792 1.692
6.521 6.512 6.506 6.553 6.443 6.089
10.44 10.70 10.90 10.48 10.38 10.20
17.11 17.00 16.98 17.29 17.11 15.07
21.23 21.31 21.40 21.49 21.26 19.17
The fibre angles h, in the eight layers are 45, !45, !45, 45, 45, !45, !45, 45°. The non-dimensional frequency parameter for the composite-layered case, is defined as follows:E h VY where D " 12(1!l l ) VYWY WYVY and h represents the total thickness of the plate. The non-dimensional frequency parameters obtained by means of the four-node conforming quadrilateral element, four-node non-conforming quadrilateral element, and three-node non-conforming triangular element, developed in this work, were compared with the analytical solution published by Narita and Leissa [19], and also compared with an experimental work published by Crawley [20], and his finite element results which were based on a moderately thick quadrilateral shallow shell element developed by Lee and Pian [21], as shown in Table 2. The package results have generally proved to be closer to the analytical solution than those based on Lee and Pian’s element. However, some deviation (+8%) between all theoretical results and corresponding experimental results was observed. Crawley [20] mentioned that the experimental results represented an average of the measured natural frequencies of nominally identical samples. Since there was a slight variation in the finished thickness from plate to plate, the measured frequencies were linearly corrected to a reference thickness before averaging. No other correction such as for fibre volume fraction, has been applied. '"ua (oh/D
4.3. Composite layered rectangular plate case This case represents a rectangular plate with aspect ratio a/b"2 : 1, where a is the length in the x-direction and b is the width in the y-direction. The plate has 8 layers of composite material with properties, fibre angles and thickness similar to the previous case but with dimensions: a"162.4 mm, b"76.2 mm. The results obtained were compared with the published experimental and finite element results of Crawley [20], as shown in Table 3. All finite element results of the package were very close to each other, but with a consistent deviation of about 8% from the experimental results. 4.4. Composite layered cylindrical shell cases It is useful to validate the package for the analysis of a curved shell, and a cylindrical shell with published results has been chosen to achieve that. The geometry and boundary conditions of the
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
481
Table 3 Natural frequencies for eight-layer composite rectangular plate Mode
Experimental results [20]
Published FEM [20]
Conf. Q.
Non-conf. Q.
Non-conf. T.
1 2 3 4 5
31.3 185.8 214.0 533.0 653.0
31.9 191.3 228.2 565.3 708.3
31.85 191.16 227.36 562.22 703.78
31.86 192.96 227.20 571.35 703.84
29.50 180.85 224.81 540.12 696.66
Fig. 3. Composite layered cylindrical shell.
cylindrical curved shell are shown in Fig. 3, and the material properties are the same as in the previous two cases. Two different shells were attempted: the first has the fibres lie in angles 0, 0, 30, !30, !30, 30°, 0, 0 and the angles of fibres in the second shell are 0, 45, !45, 90, 90, !45, 45°, 0. The results were compared with published experimental results of Crawley [20] and his finite element results based on the element of Lee and Pian [21] as shown in Tables 4 and 5. It is clear from those tables, that acceptable results have been obtained by using the developed elements. However, the results of the conforming element are more accurate than those of non-conforming elements. This case study represents a shallow shell, where the lack of stiffness values corresponding to h in the equations of the non-conforming elements will lead to ill-conditioning of elemental X rotated matrices and can affect the accuracy of the results.
482
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486 Table 4 Natural frequencies for eight-layer composite cylindrical first shell section h of layers (0, 0, 30, !30, !30, 30, 0, 0°) Mode
Experimental results [20]
Published FEM [20]
Conf. Q.
Non-conf. Q.
Non-conf. T.
1 2 3 4 5
161.0 245.1 555.6 670.0 794.0
165.7 289.6 597.1 718.5 833.3
168.72 295.09 606.06 713.01 816.94
126.34 141.74 452.03 793.62 964.10
141.27 195.37 539.76 982.72 1156.19
Table 5 Natural frequencies for eight-layer composite cylindrical second shell section h of layers (0, 45, !45, 90, 90, 45, !45, 0°) Mode
Experimental results [20]
Published FEM [20]
Conf. Q.
Non-conf. Q.
Non-conf. T.
1 2 3 4 5
177.0 201.8 645.0 754.0 884.8
192.4 236.1 705.8 808.2 980.6
196.53 246.74 720.91 813.96 985.06
190.24 244.77 674.14 1091.27 1376.14
210.59 315.18 773.32 1196.09 1441.25
4.5. Case with variation of fibre angles This case represents a rectangular plate with the length along the x-direction being a"152.4 mm and the length in the y-direction b"76.2 mm. Two plates were considered one with five layers, and the thickness of each layer is 0.13 mm, where the angles of fibres are arranged as h, !h, h, !h, h, and the second is a single layered plate with a thickness equal to the total thickness of the five-layer plate, and the fibres at an angle h. To demonstrate the effect of the fibre angle, the natural frequency analysis has been carried out on plates with different values of h, i.e. h"0, 15, 30, 45, 60, 75, 90°. The first two non-dimensional frequencies as obtained by finite element analysis with the conforming four-node quadrilateral, non-conforming four-node quadrilateral, and non-conforming three-node triangular elements, together with the corresponding analytical solution of Narita and Leissa [19], were plotted against fibre angles, for the five-layer as well as the single- layer plates, as shown in Figs. 4 and 5. Numerical results for the five-layer case showing the first five modes of vibration are demonstrated in Table 6. The following points can be deduced from the given results: (i) Results obtained using different finite elements for the two plates at different fibre angles are in a very close agreement to those obtained by the corresponding analytical solutions.
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
483
Fig. 4. Effect of fibre angle on the first natural frequency for single- and five-layer plates.
Fig. 5. Effect of fibre angle on the second natural frequency for single- and five-layer plates.
(ii) Increasing the fibre angle h, changes the natural frequency, both for the single layer and the five-layer plates, at the natural modes displayed. (iii) The change in the fibre angle h have a greater effect on the results of the five-layer plate than on those of the single layer plate. This may indicate the advantage of using more layers if the fibre angle effect is to be fully exploited. To investigate the thickness effect on the natural frequencies of a multi-layered plate, the five-layer plate was employed with two values of plate thickness being considered, the first is the same as before, but the second thickness was ten times greater. The first four non-dimensional natural
484
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
Table 6 Non-dimensional natural frequencies for five-layer rectangular plate at different fibre angle h Fibres angle (deg) h"0
h"15
h"30
h"45
h"60
h"75
h"90
Mode 1
Conf-Q Nonconf-Q Nonconf-T Analytical
3.511 3.499 3.501 3.513
3.152 3.137 3.145 3.153
2.445 2.423 2.443 2.442
1.703 1.690 1.701 1.700
1.176 1.172 1.177 1.176
0.942 0.942 0.945 0.946
0.894 0.894 0.896 0.897
Mode 2
Conf-Q Nonconf-Q Nonconf-T Analytical
7.057 7.061 7.056 7.068
8.452 8.433 8.455 8.474
10.260 10.230 10.260 10.310
10.070 10.120 10.060 10.080
7.263 7.345 7.367 7.259
5.889 5.905 5.971 5.888
5.388 5.389 5.432 5.390
Mode 3
Conf-Q Nonconf-Q Nonconf-T Analytical
21.910 21.270 21.590 21.970
19.830 20.040 19.680 19.910
15.200 15.840 15.270 15.260
11.130 11.430 11.270 11.190
9.431 9.485 9.485 9.476
6.921 6.938 6.990 6.931
5.601 5.606 5.701 5.602
Mode 4
Conf-Q Nonconf-Q Nonconf-T Analytical
25.230 25.150 26.260 25.840
27.110 26.570 27.480 27.490
31.410 30.920 31.220 31.470
28.980 29.290 28.960 28.980
20.6 21.05 21.2 20.62
16.5 16.6 17.1 16.5
15.7 15.70 16.3 15.7
Mode 5
Conf-Q Nonconf-Q Nonconf-T Analytical
27.760 26.030 27.360 27.840
30.850 30.010 31.27 31.340
36.520 36.190 36.850 37.120
33.870 34.610 34.020 33.950
29.010 29.260 29.190 29.090
21.450 21.540 21.790 21.490
16.930 16.940 17.300 16.950
Fig. 6. Effect of thickness on the first natural frequency for five-layer plate.
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
485
Fig. 7. Effect of thickness on the second natural frequency for five-layer plate.
frequencies for the thin and thick plates were plotted against fibre angles h, as shown in Figs. 6 and 7. It is clear from those figures, that the transverse shear effect tends to reduce the non-dimensional frequency, but this reduction is greater at h values between 15 and 60°, approximately. It is also clear that the effect of the fibre angle on natural frequency may be reduced by increasing the plate thickness.
5. Conclusions It is clear from the previous case studies that conforming and non-conforming elements developed in this work have led to accurate and reliable results. Although, the difference between their corresponding results was insignificant for flat plate cases, the conforming elements have shown more accurate results for the curved shell case. The fact that the conforming element has stiffness and mass contributions corresponding to h , h , and h will make them well conditioned, V W X and provide more accurate results for shallow shell cases than those obtained from non- conforming elements. The conforming elements are also more flexible in use with curved shells and box structures, because they require boundary conditions to be in terms of slope angles with respect to global axes, whilst non-conforming elements require them in terms of local axes, which are not easy for the user to know. Although the developed elements consider transverse shear effects, they do not suffer from shear locking as do Mindlin-type elements, and reduced integration techniques are not required for the elements developed in this work. Transverse shear effects, which are increased by increasing thickness, tend to reduce the non-dimensional frequencies. The fibre angle h can have a significant effect on the stiffness matrix of composite layered plates and shells, and this will result in a corresponding effect on their natural frequencies. This effect can be increased by increasing the number of layers, while keeping the same thickness. The h effect may
486
O. Attia, A. El-Zafrany / International Journal of Mechanical Sciences 41 (1999) 461—486
provide a tool for the designers to tune the natural frequencies away from resonance due to excitation frequency, by simply changing the fibre angles without having to modify the mechanical design.
References [1] Epstein M, Huttelmaier HP. A finite element formulation for multilayered and thick plates. Computers and Structures 1983;16:645. [2] Widera, OEG, Moumene, M. Cutout in laminated plates and shells. Advances in Macro-Mechanics. Composite Material Vessels Components 1984;146:155. [3] Has, TM, Wang JTS. A theory of laminated cylindrical shells consisting of orthotropic laminae. AIAA Journal 1970;8:2141. [4] Witt M, Sobczyk K. Dynamic response of laminated plates to random loading. International Journal of Solids Structures 1980;16:231. [5] Pillasch DW, Majerus JN, Zak AR. Dynamic finite element model for laminated structures. Computers and Structures 1983;16:449. [6] Nelson RB, Lorch DRA. Refined theory for laminated orthotropic plates. Journal of Applied Mechanics 1974;41:177. [7] Reddy JN. A refined nonlinear theory of plates with transverse shear deformation. International Journal of Solids and Structures 1984;20:881. [8] Reddy, JN. A simple higher order theory for laminated composite plates. Journal of Applied Mechanics 1985;51:745. [9] Wang FY. Two dimensional theories deduced from three-dimension theory for a transversely isotropic body—I. Plate problems. International Journal of Solids and Structures 1990;26:445. [10] Chang B, Chen DC, Shabana AA. Effect of the coupling between stretching and bending in the large displacement analysis of plates. International Journal of Numerical Methods in Engineering 1990;30:1233. [11] Wung PM, Reddy JN. A transverse deformation theory of laminated composite plates. Computers and Structures 1991;41:821. [12] Tessler A. An improved plate theory of +1,2,-order for thick composite laminates. International Journal of Solids and Structures 1993;30:981. [13] Cho M, Parmerter R. Finite element for composite plate bending based on efficient higher order theory. AIAA Journal 1994;11:2241. [14] Crisfield MA. Non-linear finite element analysis of solids and structures. New York: Wiley, 1991. [15] El-Zafrany AM, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for triangular elements. International Journal of Numerical Methods in Engineering 1986;23:275. [16] El-Zafrany AM, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for quadrilateral elements. International Journal of Numerical Methods in Engineering 1986;23:1939. [17] Zienkiewicz OC. The finite element method, 3rd ed. New York: McGraw-Hill, 1977. [18] El-Zafrany AM, Fadhil S, Debbih M. An efficient approach for boundary element bending analysis of thin and thick plates. International Journal of Computers and Structures 1995;56:565. [19] Narita Y, Leissa AW. Frequencies and mode shapes of cantilevered laminated composite plates. Journal of Sound and Vibration 1992;54:161. [20] Crawley EF. The natural modes of graphite/epoxy cantilever plates and shells. Journal of Composite Materials 1979;3:195. [21] Lee SW, Pian TH. Improvements of plate and shell finite elements by mixed formulations. AIAA Journal 1978:16.