- Email: [email protected]
A general quadrilateral multilayered plate element with continuous interlaminar stresses
Compurers&Structures Vol.41, No. I,pp. 91-105, 1993 0
PrintedinGreat Britain.
aw-7949/93 s6.00+0.00 1993 Rrgamonl'ms Ltd
A GENERAL QUADRILATERAL MULTILAYERED PLATE ELEMENT WITH CONTINUOUS INTERLAMINAR STRESSES-f M. DI SCIUVA Department of Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruxxi 24, 10129 Turin, Italy (Received
16 December
1991)
Abstract-The purpose of this work is the formulation of a general quadrilateral multilayered anisotropic plate element. The plate element is formulated on the basis of a refined third-order shear deformation plate theory recently proposed by the author. The plate model makes use of a displacement field that fulfils a priori the geometric and stress continuity conditions at the interfaces between the layers, and requires only five generalized displacements to describe the kinematics of the laminate deformation. To test the performance of the developed plate element, both closed-form and finite element-based solutions are given and compared with results from three-dimensional elasticity and other approximate two-dimensional plate models.
linear or nonlinear distribution of the in-plane displacements or the transverse shearing stresses in the thickness direction. These approaches consider all layers as one equivalent single anisotropic layer. CPT, based on neglecting transverse shear strains and transverse normal strains and stress in the laminate [7], and the first-order shear deformation theory (FSDT), in which the transverse shear strains are assumed to be constant through the thickness [8] and the transverse normal strain and stress are neglected, are well-known examples of plate models with assumed global linear distribution. The so-called higher-order shear deformation theories (HSDT) are examples of plate models with assumed global nonlinear distribution [g-16]. This assumption generally results in incompatible shearing stresses between every two adjacent layers. In fact, from the theory of elasticity we know that the displacements and stresses at the interface k between the kth and (k + 1)th bonded layers must satisfy the following contact conditions$
1. INTRODUCTION The use of fibre reinforced laminates in aerospace, automotive, shipbuilding and other industries has increased tremendously in the past several years. This
is largely due to the high strength-to-weight ratios of composites as well as their ability to be tailored to meet design requirements of strength and stiffness. Coinciding with these new applications is the interest in the accurate prediction of the detailed response and failure characteristics of laminated anisotropic plates. The anisotropy and nonhomogeneity and larger ratio of longitudinal to transverse moduli of these new materials demand improvement in the classical plate theory (CPT), plate theory based on the well-known Kirchhoff assumptions. As a result, modelling approaches which take into account transverse shear deformation have been the topic of serious research in the last two decades. It suffices here to mention the survey papers of Bert [1], Noor and Burton [2,3], Kapania and Raciti [4, 51, and Reddy [6]. Some of these modelling approaches are extensions of similar approaches used for isotropic plates. Belonging to this category are the so-called smeared laminate models, that is, those models based on an assumed global (i.e. one for the whole laminate)
geometric continuity
conditions:
ri,(z;) = ti,(z$) stress continuity &,(z;)
t This paper represents an updated and completed version of the one presented at the European Conference on New Advances in Computational Mechanics, Giens (Var)-France (April, 1991). $ In this paper, if not otherwise stated, the Einsteinian summation convention over repeated subscripts will be adopted with Latin indices ranging from 1 to 3 and Greek indices from 1 to 2.
conditions:
= e&:
);
G(ZL ) = &,(zk+ ).
In these relations, zii (i = 1-3) are the displacement components along an orthogonal frame xi on the reference surface Q, xj = z being the normal to R. The minus and plus refer to the values of the attached functions for z = zk - 0 and z = z, + 0, 91
92
M. DI SCIUVA
respectively, where z, is the value of the z-coordinate at the interface k between the layers k and k + 1 (Fig. 1). From the geometric continuity conditions follow the strain continuity conditions &,a@;) = &,a@: ). Continuity conditions cannot be imposed on the remaining stresses and strains; in general, these are discontinuous at the interfaces. If the transverse shearing strains are assumed to be continuous functions of the thickness coordinate z, it is not difficult to realize that the contact conditions on the transverse shearing stresses cannot be satisfied if the layers have different mechanical properties. In the literature there are a number of approaches in which the contact conditions are satisfied. These models are based on an assumed piecewise linear or nonlinear displacement or stress distribution in the thickness direction and are known as discrete-layer models [17-291. Most of these models consider the single layer as a plate and impose the contact conditions as constraint conditions on the transverse stresses at the interfaces. The major drawback is that the number of governing equations increases with the number of layers in the laminate. The approach proposed by the author [21-251 overcomes this drawback as it makes use of an assumed displacement field that fulfils a priori the geometric and stress continuity conditions at the interfaces between the layers, and requires only five generalized displacements to describe the kinematics of the laminate deformation, as in FSDT. In its original version, a piecewise linear variation of the in-plane displacement (the so-called zig-zag model [21-231) has been assumed. This model may be considered the counterpart of the FSDT in the class of the discrete-layers models, and as with the FSDT, it is unable to fulfil the conditions of zero transverse shear stresses on the bounding surfaces of the plate. To overcome this restriction, the author [24,25] proposed an improvement of the zig-zag model (the so-called refined higher-order shear deformation theory, RHSDT model), by assuming a displacement field which allows a nonlinear variation of the inplane displacements through the laminate thickness and fulfils a priori the geometric and stress continuity conditions at the interfaces and the static condition of zero transverse shear stresses on the top and bottom surfaces, at least for symmetric laminates. Moreover, only five generalized displacements are required to describe the kinematics of the laminate, as in the FSDT and in the zig-zag model (in the following called the RFSDT). In the last two decades, there has been a considerable effort towards the development of shear deformable, multilayered anisotropic plate elements. This
development follows the same path as in the analytical modelling approaches. In effect, as a basis for the formulation of plate elements of this type, many investigators have employed smeared laminate models [30-40], while others make use of discretelayer models [35,36,41-511. In the first approach, the multilayered anisotropic plate element is reduced to an equivalent single anisotropic layer element. Many of the finite element formulations belonging to the second group make use at least of one solid or plate element per layer in the thickness direction. Then, this approach requires numerous unknowns and, as a consequence, these formulations become nontractable in the solution of practical engineering problems. The rectangular plate element formulated by the author [35,36] on the basis of the RFSDT model, zig-zag model approach, overcomes this drawback because the number of nodal parameters is independent of the number of layers. Based on the RHSDT model [24,25], in this paper we formulate a general quadrilateral multilayered anisotropic plate element. To test the performance of the developed plate element, both closed-form and finite elementbased solutions are given and compared with results from three-dimensional elasticity and other approximate two-dimensional plate models. The paper is structured as follows: based on the assumed piecewise cubic displacement field and on the von Karman nonlinear strain-displacement relations, in Sec. 2 we derive the nonlinear equations of motion and the static stability equations of multilayered anisotropic plates by making use of the dynamic version of the principle of virtual work. In Sec. 3 we formulate the quadrilateral plate element. To give an insight on the accuracy and reliability of the developed finite element, Sec. 4 contains analytical and FEM results on bending under sinusoidal transverse loading, undamped frequencies and buckling loads under uniform axial compression. Whenever possible, a comparison with results from three-dimensional elasticity and other approximate two-dimensional models is also given. 2. ANALYTICAL
FORMULATION
2.1. Kinematics Let us consider a multilayered plate of total thickness h and reference surface R composed of N orthotropic layers of uniform thickness perfectly bonded together. The principal axes of elasticity of each layer can be arbitrarily oriented with respect to the laminate orthogonal reference frame xi, where x, is an orthogonal frame on the reference surface Q (here chosen to be the midsurface), and xj = z is the normal to R (see Fig. 1). The material properties and the thickness of each layer may be entirely different. The discrete-layer plate model considered here is based on the following piecewise cubic
A general quadrilateral anisotropic plate element
93
Fig. 1. Geometrical laminate configuration. approximation of the in-plane displacement the plate thickness
across
r&(xa, 2; t) = u.(xs; t) - 2 &w(xfl; t) +f,&)g&;
t)
r&(xfl, 2; t) = w(xg; t) or, in matrix notation ii = Z,d, where I’= [u”, r& r&l, d: = [d; dj, d;] d; = b,
4,
d;=[w
S;],
@:,=[a,~
determined; z, are the coordinates of the N - 1 interfaces; Yk is the Heaviside unit function (it has the value 0 for z < z,, and 1 for z 3 zk). The tracers 6, and 6, identify the contributions through (dF= 1; CRT (6, = 8r = U.jk = 0), FSDT-modified 6, = aDlP= 0), HSDT (6, = 8r = 1; a@ = 0), RFSDT (a,= 1; 6,= 0) RHSDT (S,= 6, = 1) plate models. The notation a,( .) represents the partial derivative a(.)/aa. Note that the displacement components li, are continuous functions of the thickness coordinate z for all values of the parameters aplr . Thus, the geometric continuity conditions are automatically satisfied and we must only satisfy the stress continuity conditions on u,,. By satisfying these stress continuity conditions, it may be shown [24,25] that the following expressions for a,, hold?
a2w] k-l
k-l
%lk
al
d;=k,
0
0
--z
[ 001
0 0
%lq
q-1
1 0 0 --z s2”,=01
x ullq 1 +AIU, 1 q=l
tal)k+
=Allk
0
.flI A* fil122. 0
0I
k-l
auk --&k
caJ)k
+
taf)k
+
caf)k
+
x q=I
hq
k-l
N-l
u22k =
A22k
k=l
with
~22,
k-l
The functions fur(z) specify the distribution of the transverse shear strains along the thickness. In this study we make use of the following expression [24,25]
Am(z) =W@) + c %k(Z-%)Yll
k-l
1 q=l
1 q=I
U22q
1
+&k
1 1
c q”
k-l
+
&k
+
A2,k
c q-1
Ilk
x q-1
&%d?::k+,-Qt:,Q::,+, _, gk
Here u.(xP; t) and w(x#; t) are the displacements in the x, and z directions, respectively, of a point belonging to R; g,(xs; t) are the transverse shear rotations for z = 0; Sf is the Kronecker delta (St = 1 if a = /?; St = 0 if a #/I); a,, are parameters to be
A
I2k
=
Q::kQ::k+,
-
Q::kQ::k+,
Bk
A2x=Q::kQ:k+,
-Q::kQf:k+,_, 9k
t Obviously, the expressions for u~,,~are related to the assumed expressions for the strains. Those given here hold for the strain-displacement relations given in Sec. 2.2 and for the lamina constitutive relations given in this section.
U21q
k-l
where A
%q
st,=(Qt:k+,)‘-Ql:+,Q::k+,.
UI~,
94
M. DI SCIWA
Thus, for a specified laminate, uegk are known constants only depending on the transverse shear mechanical properties of the constituent layers. In the above relations, the notation (afit represents a,j valued at z = z,. Moreover, Q$ and Qz& are the so-called reduced stiffnesses of the k th orthotropic layer in the plate coordinate system (x,, z). They relate the stress components d to the strain components P (lamina constitutive relations)
respectively. Notice that we have used c for the linear portion of the strain. We can write i2=E+iiNL with 2ENL= A W’ 8 IV
C = %“,dc,
where -zU
or, in matrix notation
F
0
0
dF
1
and
a=Qg, where
I Qft
Q=
Q::
0
Qt:
0 0 0 Q::
Q:: Q:: Q:: 0 Qi: Q:: Q!: 0 0 0 0 Qt: 0 0 0 Q:: Q::
Q is the (5 x 5) matrix of the transformed stiffnesses of the generic layer. 2.2. Strain-displacement
e; = A;=
reduced
0
a2
0
0
a,
0
a,
a,w 0 0
azw
1
1.
a,w
0
0
a,w
0
0
Moreover. let
e;000 0e; I1 Qf= 0e; 00 d, = 53:ds
relations
In deriving the static stability equations, we make use of the von Karmin strain-displacement relations (moderately large deflections and small strains)
[
a,
with
i
and
where
where d, is the vector of the generalized displacements d: = [d;
w
d;].
and 26a8=a,u,+a,u,,
k,,=
-asBw
where L,~ and k,, are the (membrane) extensional strains and the curvatures of the reference surface,
2.3. Governing equations
In the present study we are interested in deriving the linear equations of motion and the static stability equations. Both the problems can be approached by making use of the principle of virtual work 6# = 6L’+
aLin,
95
A general quadrilateral anisotropic plate element where 6@ is the virtual variation of the strain energy; 6L’ is the virtual variation of the work done by the applied external forces; 6L” is the virtual variation of the work done by the inertial forces.
where &I+, =
JR
SG2=
JR
2.3.1. Strain energy. Substituting the expressions for the von Kbrmin strains in the expression of the virtual variation of the strain energy 6@ =
(a-‘~%) dR, s Cl
where B is the stress state corresponding to a generic time of motion and 60 the strain state related to the virtual displacements Slii, we obtain
S@ =
Jn
W’6d,dQ+
d:Xo Sdc di2 =
(97;d,)‘Xo(ka;
ad,) dR
8 :N 68, dQ.
In terms of force and moment stress resultants the explicit expression for the virtual variation of the strain energy of the multilayered anisotropic plate reads
Jn
@;N&‘,dfl
with
+
&N&ia,4 6w+ (ass,,- T,) k,l dfi
+
Jr
w = (%2”:d) being the vector of the force and moment resultants for unit length
stress
~‘=[NI,;N,~;N~,;N,,;M,,;M,,;M~,;M~~; S,,;
Jfl
[Nm6u, - Mg, da,
W +
S, Sg,
+ (v, + a,M,,, + N,,,) 6w] 6r, where V, = n, a,M,,
and
Sn; S,,; 8,; T,; Tzl = IN’;M’;S’;T’l,
where ( . )n, =
(N,,, Kb, S,,) = K&)(1>
z&r)>
T, = and
hi2 (...)= J’ J (...)&=;
Here, and in the following, t stands for transpose and
4
( . )2n -
n2(.
)*n
on r, r denoting the boundary of R, and n, the direction cosines of the outward unit normal n to r with respect to the a-axis and t the tangent unit vector. Let S,,(s) be the 2.3.2. Work of external forces. portion of the lateral cylindrical surface of the sth layer on which the external loadings per unit area A(s) are prescribed; let /jZ =J%(x,, -h/2) +pZ(x,, +/z/2) be the transverse loading per unit area. Under these conditions, it is found that
(...)dz.
-h/2
s=l
2%- I
Moreover, by taking into account the constitutive equations, we obtain W = .%,,d, + +f-, A,,@,.,
where we have introduced the following resultants and moment resultants of the external forces
where XO =
W’:Q~‘,),
x, = (S:Q>.
Then 2.3.3. Work of internaI forces. By using the approximate expressions for Sri, (as well as for &), an
M. Dr !%XJVA
96
approximate expression for the virtual variation the work of inertial forces is obtained as
of
where
equal to zero, the following nonlinear motion result
@g,)
h=(s:P~“) and 53: is a matrix of differential that
operators
such
d, = B:d,. We have
a,S,,-
T,=m&,-m~asti
equations of
+m$&.
Boundary conditions : Let ru be the portion of r on which are prescribed the values of the generalized coordinates and let r, be that portion of r on which are prescribed external tractions. Thus the set of boundary conditions reads Prescribed on r.
Natural on t’,
u, = ii w= ;
N,=Nm v,+~,M,,+N,,=
+m,ii,-m,d,G M”=iV” s,=S,.
w,n = C’,” g.=&
v=+a,li;i,,
+nsm@+
2.4. Static stability equations 0000010 0100001 In expanded form and after performing grations by parts, the result is
[(mot&- m, a,$ + mf&)
6L” = -
some inte-
6u,
sR
where we have introduced the following inertia resultants (i =O, 1,2) (m,; m?; m$) = (P(z’; z%~;_L_t$,s)>. In the previous relations, p denotes the material mass density, and the overdot indicates differentiation with respect to time. Equations of motion: By substituting the expressions for SQ, 6L’, and 6L” into the principle of virtual work, rearranging the various contributions and setting the coefficients of the virtual variations of the generalized coordinates in the integral over R
The static stability equations of the laminated plates subjected to combined in-plane normal and shear loads may be formulated by using the Euler method of the adjacent equilibrium configurations [52]. In this approach, we assume the existence of at least one additional, distinct equilibrium configuration (the adjacent equilibrium configuration) in the very close neighbourhood of the original equilibrium configuration (the equilibrium configuration just prior to the occurrence of the buckling). If such an adjacent equilibrium configuration exists, the body may change suddenly from the original equilibrium configuration (here characterized by ui, w”, etc.) to the other under the stimulus of infinitesimal disturbance u,, w, etc. Although these terms have been used to denote the displacement field in the previous derivation of nonlinear equations of motion, it is clear that in the present stability equations they assume a different significance. More precisely, in the nonlinear equations of motion u,, w are the displacements measured from the reference stress-free state, while here they are measured from original equilibrium configuration. Using the expressions for the nonlinear strains as given by von Karrnbn, on the assumption that the external forces and moments in both equilibrium configuration are identical, the principle of virtual work of the adjacent equilibrium configuration reads
s
((a-“’ W’ + d” 6’))
n
dR = 0,
A general quadrilateral anisotropic plate element where the prime and double prime represent the quantities which contain terms that are linear or quadratic in the derivatives of the infinitesimal disturbance displacements P’=~+~(A,&,+A,8~),
97
and
2Z”=A,8,.
Noting that
Let us assume a linearized form for the relations between the incremental stresses and strains, i.e. d’ = QG’, then
A,,&, = A,@@ yields C=i+A,@,. By taking the virtual variations with respect to the disturbance quantities, yields
W’=A,68,.
where X0 and X, have been defined previously, and X2 = (Q). In explicit form
Substituting the expressions for the strains in the expression of the principle of virtual work for the adjacent state, we obtain
+ where &Do=
(,‘I
Si) dR
SW’= -
+
(tf”.2Z’~) Sd, dR
sn =
s n
8”6d EdR
sr
n, $,(N:swo) dT
s s
(N:pa,pw +d,N~fld,w)Sw
di-2
0
sn =
(N& i&w0 + d,N;, 8,~‘) SWd&2
6@’ = -
n, as(N;, w) dT.
=
The expression for S@Ois formally the same as that given in Sec. 2.3.1. 3. THE FINITE ELEMENT FORMULATION
This section deals with the displacement finite element formulation of the multilayered plate models summarized in the previous section. Let the region of the plate be divided into a number of quadrilateral elements with straight sides (Fig. 2).
rl’ (-1.+1)
@
* 5
(-1.4)
0
Fig. 2. Mapping of the geometry.
@
(+1.-l)
M. DI
98
.%JlJVA
Table 1. Centre deflection I? of a square plate, laminate ORT3LU, simply-supported on all edges under transverse s&so&l loading. I+,-~, = 0.43 1
Table 2. Centre deflection $ of a rectangular plate (a2 = 3a,), laminate ORT3LU, simply-supported on all edges under transverse sinusoidal loading. J,, = 0.503
AT
FSDT
HSDT
RFSDT
RHSDT
AT
FSDT
HSDT
RFSDT
RHSDT
Ref. [53]
4 10 20 50 100
1.776 0.669 0.492 0.441 0.434
1.922 0.713 0.504 0.443 0.434
1.891 0.723 0.508 0.444 0.434
1.972 0.754 0.517 0.445 0.435
4 10 20 50 100
2.363 0.803 0.578 0.515 0.506
2.641 0.862 0.594 0.518 0.507
2.717 0.881 0.599 0.519 0.507
2.757 0.919 0.610 0.520 0.508
2.820 0.919 0.610 0.520 0.508
The four-node quadrilateral plate element formulated on the basis of the RHSDT has 10 DOF at each nodal point, thus giving a total of 40 DOF. In the following, this plate element will be named RHQ40. The nodal parameters are the five generalized displacements ul, w, g,; the two total rotations 19,w; the two curvatures &,w and the twist a,w of the reference surface a. Notice that following the scheme of Sec. 2, it is possible to derive from the RHQ40 plate element the corresponding plate elements formulated on the CPT (here named CQ40), on the FSDT-modified (here named FQ40), on the RFSDT (here named RFQ40), and on the HSDT (here named HQ40). 3.1. Shupe functions To make the computation of the involved integrals uniform, a topological transformation (mapping) of the geometry element from the physical plane x, , x2 to the transformed plane 5, r) (Fig. 2) has been accomplished using the Lagrange linear interpolation functions
principle) for the plate mode1 considered here (CPT, FSDT-modified, HSDT, RFSDT, and RHSDT) requires C’-approximation of the transverse deflection w, and Co-approximation of the in-plane displacement w, and the shear rotations g,. Thus, Hermite cubic or higher-order interpolation functions should be used for approximating w, and linear or higherorder Lagrange interpolation functions for approximating u, and g,. In formulating the plate element RHQ40, the Lagrangian linear interpolation function Li has been used for approximating d,, and d8; Hermite sixth-order polynomials Hi for approximating w. Thus, an isoparametric formulation has been used for d, and d,, and a subparametric formulation for w. Note that the element is able to model a constant strain state and a rigid body displacement field in the membrane and transverse behaviour. Moreover, the element is conforming because of the continuity of the displacements and the slopes across adjacent elements. Let A: = [A: A;,. A;],
where .Xi = a (1 + tit) (1 + qiq), x,.~being the coordinates of the nodal point i in the physical plane and &, vi those in the transformed plane. This transformation maps a generic quadrilateral region with straight sides in the physical plane into a square region of boundary [- 1 < 5,~ < 11 in the natural plane. The finite element formulations based on the displacement approach (i.e. on total potential energy Table 3. Maximum
with A:= [ul
41, A; = k’, i&l,
Al. = [w’
alwl
T 0.437 TO.734 TO.620 TO.845 -0.755 +0.801
a12d
azzwq
be the vector of the physical degrees of freedom surrounding element e. Each of the 10 subvectors has
stresses in a square plate, laminate ORT3LU, transverse sinusoidal loading
FSDT HSDT RFSDT RHSDT Ref. [53]
a2Wl a,,wl
simply-supported
TO.477 TO.503 TO.661 TO.598 -0.556 +0.534
TO.037 TO.050 + 0.045 TO.051 -0.0511 +0.0505
+0.256
0.217
+0.381 +0.369 +0.364 +0.357 +0.357
+0.111 +0.117 +0.123 -0.124 +0.1228
+0.395
0.082
FSDT HSDT RFSDT RHSDT Ref. [53]
TO.513 TO.568 TO.549 TO.595 -0.590 f0.590
TO.254 TO.269 T 0.289 TO.292 -0.288 +0.285
70.025 T 0.028 TO.027 7 0.029 -2.89 +2.89
CPT
TO.539
TO.180
TO.021
0.037 0.286 0.262 0.241
on all edges under
0.197 0.209 0.254 0.230
A general quadrilateral anisotropic plate element
99
Table 4. Maximum stresses in a rectangular plate (a2 = h,), laminate ORT3LU, simply-supported on all edges under transverse sinusoidal loading AT
4
10
elements sponding
Model FSDT HSDT RFSDT RHSDT Ref. [53]
011 TO.613 T I .036 TO.925 f 1.205 -1.10 +1.14
512
u22
F 0.093 TO.103 TO.114 TO.108 -0.119 +0.109
T 0.020 TO.026 F 0.026 T 0.028 - 0.0269 f0.0281
013
93
0.436 0.383 0.366 0.329
0.028 0.030 0.032 0.030
+0.351
+0.0334 +0.014 +0.015 +0.015 +0.015 +0.0152 f0.0108
FSDT HSDT RFSDT RHSDT Ref. [53]
TO.621 TO.692 T 0.675 TO.730 -0.725 +0.726
TO.037 F 0.040 T-o.041 TO.042 -0.0435 +0.0418
TO.010 TO.012 TO.012 TO.012 -0.0123 +0.0120
+ 0.439 +0.430 + 0.427 +0.419 + 0.420
CPT
f 0.623
TO.0252
T 0.0083
+0.440
which are the nodal values variable; for example II; = [uf
u:
u;
of the corre-
and
“y;=X*=
uf].
where the submatrices 1x4
N-0 [ 0
.N
1
are row matrices
of dimension
he the appropriate matrices of the chosen shape functions in terms of the natural coordinates 5, r~such that the following formal relations hold for the element e =A’&, = [Np], with i = l-4. The explicit expressions for the elements of these matrices are given in the Appendix. 3.2. Element elastic st$ness matrix The element elastic stiffness matrix K, is obtained from the linear part 69’ of the strain energy 6@ (see Sec. 2.3.1). Taking into account that in the present finite element formulation
Then
where yields
Table 5. Frequency parameter of the fundamental bending mode cG of a square plate, laminate ORT3LU, simplysupported on all edges AT
4 10 20 30 50 100
CF’T
FSDT
HSDT
RFSDT
RHSDT
14.501 15.104 15.197 15.214 15.233 15.227
7.413 12.163 14.230 14.757 15.053 15.183
7.116 11.790 14.060 14.670 15.019 15.175
7.163 11.700 14.012 14.645 15.009 15.172
7.027 Il.462 13.888 14.579 14.983 15.165
Table 6. Asymptotic buckling stress parameter K, of a square plate, laminate ORT3LU, simply-supported on all edges under uniaxial compression fl,, AT
4 10 20 30 50 100
CPT
FSDT
HSDT
RFSDT
RHSDT
23.495 23.495 23.495 23.495 23.495 23.495
5.706 15.138 20.588 22.101 22.972 23.362
5.272 14.221 20.099 21.843 22.870 23.362
5.358 14.008 19.961 21.767 22.839 23.327
5.138 13.436 19.608 21.572 22.760 23.306
M. DI SCIWA
loo
Table 7. Maximum normalized deflection iir for a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges, under transverse sinusoidal load. Numerical integration scheme with 3 x 3 and 5 x 5 Gauss points. Convergence analysis
NE
3x3
1 2 3 4 5 Anal.
3x3
5x5
3x3
5x5
0.690 0.709 0.711 0.712 0.712
0.675 0.705 0.710 0.711 0.711
0.726 0.750 0.752 0.753 0.754
0.708 0.745 0.750 0.752 0.753
5X5
0.700 0.720 0.722 0.723 0.723
RHQ40
HQ40
RFQ40
0.684 0.716 0.720 0.722 0.722
0.713
0.723
0.754
Table 8. Fundamental frequency parameter CT,for a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges. Numerical integration scheme with 3 x 3 and 5 x 5 Gauss points. Convergence analysis HQ40
RFQ40
RHQ40
NE
3x3
5X5
3x3
5X5
3x3
5X5
1 2 3 4 5 Anal.
11.848 11.695 11.680 11.675 11.673
11.982 11.727 11.694 11.683 11.678
11.990 11.839 11.825 11.820 11.818
12.116 11.870 11.838 11.827 11.823
11.644 11.475 11.458 11.453 11.450
11.793 11.511 11.474 11.462 11.456
11.700
11.463
11.790
where
From the previous positions, it is readily found that
B, = %r”9;J(r is the strain matrix relating the strain components of the DOF in the physical plane. In the above expression, J is the Jacobian of the transformation from the natural to the physical plane given previously, and the matrix JT relates the partial derivatives of the shape functions with respect to the natural coordinates to the partial derivatives of the shape functions with respect to the physical coordinates
where
is the matrix which relates the slopes reference surface to the nodal degrees the transverse displacement w. Upon the expression for 8, into the virtual
of the element of freedom for substitution of variation of the
strain energy 6@“, yields
where the definitions follow from
for the submatrices
J,,
etc.
and “0, , ‘VI2have the same form as 8,) 0, by replacing 1 with e, and 2 with 1. The same holds for “PS. Then I
I
ss-I
-I
K,=
B;X-,B,J
dt dq.
3.3. Element geometric st@ness matrix As is well known, the element geometric stiffness matrix Kge is obtained from the expression of the virtual variation of the strain energy due to the nonlinear terms in the incremental strains S@: (Sec. 2.4).
/.-.,2-wj
Fig. 3. Mesh with NE = 5.
*I -
A general quadrilateral anisotropic plate element
101
Table 9. Centre deflection I? of a square plate, laminate ORT3LU, simply-supported on all edges under transverse sinusoidal loading. Effect of the thickness ratio on the discretization error AT
FQ40
(Anal.)
HQ40
(Anal.)
RFQ40
(Anal.)
RHQ40
(Anal.)
4 10 20 100
1.566 0.630 0.482 0.433
(1.568) (0.63 1) (0.521) (0.433)
1.917 0.712 0.504 0.434
(1.922) (0.713) (0.504) (0.434)
1.887 0.723 0.507 0.434
(1.891) (0.723) (0.508) (0.434)
1.966 0.754 0.517 0.435
(1.972) (0.754) (0.517) (0.435)
from which follows the expression geometric stiffness matrix
In order to show the accuracy and reliability of the RHSDT model and to assess the convergence characteristics of the corresponding finite element formulation (the RHQ40 finite element), some numerical investigation has been performed. In the following, analytical and FEM results for the following three sample problems are given. 1. Bending under sinusoidal transverse loading
where
(%AQ’N”@N,)J
K,=
dr d+
The consistent element mass matrix is readily obtained from the expression for 6L” by taking into account that
/!?A,
where m and II are the number of halfwaves in the x, and x,-directions, respectively; a, and aI, plate dimen-
sions in the x, - and x2-directions, respectively; and, A, maximum intensity of the distributed loading. 2. Undamped frequencies.
00 f,u= 00I9-1 001 d, = duCf%.K)A,,
3. Buckling under uniform in-plane axial compression and shear of rectangular and square laminated plates. For purposes of comparison, the evaluation
I
i
a,=:,
pz =pO sin u,x, sin j$xz,
3.4. Consistent element mass matrix
where
4. NUMERICAL RESULTS
for the element
and “9: is the same as $’ with 1 replaced by < and 2 replaced by r,t in the derivatives. Moreover, f-’
denotes the inverse of the Jacobian matrix. The result is
has been conducted on sample problems for which exact elasticity solutions (ES) are available. In addition, numerical results from CRT, FSDT with shear correction factor x = 5/6, RFSDT, and HSDT, are also given. Further numerical results showing the accuracy of the proposed approach can be found in [25]. 4.1. Analytical results
Explicit expressions for the elements of the elastic stiffness, geometric stiffness and mass matrices are not given here. It is sufficient to note that the formulation upon which the expression for the element mass matrix rests also includes the in-plane and rotatory inertia terms.
The numerical results here given refer to a threelayered, symmetric cross-ply plate simply supported on all edges for which the exact elasticity solution has been given [53]. The layers have equal thickness with fibres in the outer layers oriented in the x,-direction (0 degree) and those in the inner layer oriented in the x,-direction (90 degrees). The mechanical properties of each layer are RE = EL/ET= 25; RG = GLT/ Er = 0.5; RT = G&E, = 0.2; vLr = 0.25; here L and T are the directions parallel and normal to the fibres, respectively, and vLr is the Poisson’s ratio measuring transverse strain under uniaxial stress parallel to the fibres. In the following, this laminate
Table 10. Centre deflection ii~of a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges under transverse sinusoidal loading. Comparison between RQ4 and RFQ40 elements
Table 11. Fundamental frequency parameter (I, for a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges. Comparison of the convergence characteristics of the RQ4 and RFQ40 elements
where M,=
NE R:$O Anal.
C’~,“.NYd%WdfwWi’4Jdt dv
1
2
3
4
NE
1
2
3
4
0.606 0.700
0.722 0.720 0.723
0.723 0.722
0.723 0.723
RQ4 RFQ40 Anal.
12.938 11.982
11.836 11.727 11.700
11.734 11.694
11.697 11.683
M. DI
102
Table 12. Transverse deflection and buckling loads for a four layers, symmetric (45/-45), angle-ply square plate (AT = lo), simply-supported on all edges (NE = 6 for the whole plate)
FQM RHQ40 Galerkin
I?
N,,,.
8.208 8.439 8.090
18.535 18.110 19.230
SCIUVA
Table 13. Frequency parameter of the fundamental bending
mode (I, for a two layer, antisymmetric (O/90) cross-ply square plate (AT = 10). simply-supported on two opposite edges and clamned on the other two sides NE
R, ,“, 23.133 22.129 23
-60.808 - 54.704 -59
Ref. 1541
will be called ORT3LU. As an example of the accuracy of the proposed RHSDT model, Tables 1 and 2 give results of the maximum deflection
3x3
1 2 3 4 5 Anal. Ref. 1551
of the fundamental load parameter
5x5
29.53854 17.41517 16.48013 16.17408 16.02959 15.709
30.95164 18.44556 17.03638 16.54957 16.31439 15.709
bending mode and the buckling
E,h’ 3 = lOOA,~ Peal
as a function of the plate thickness ratio AT = q/h, for the following simple supported boundary conditions XI
x,=O,a,:
4.2. FEM results
N*,=M,,=S,,=u*=w=g~=0
0, a, :
=
N22=M22=S22=u,=w=g,=0.
For the same plate, Tables normalized maximum stresses
3 and 4 give the
%3 _ *23 =
DoAT
Notice the high accuracy of the proposed approach. Also notice that the results from RFSDT, HSDT, and RHSDT are obtained without introducing any arbitrary shear correction factor.
’
Tables 5 and 6 give results for the undamped quency
fre-
In order to show the convergence characteristics of the RHQ40 finite element, Tables 7 and 8 summarize the FEM results for the centre deflection 5 and the undamped frequency parameter @.Due to symmetry, only one quarter of the plate has been modelled. In the tables, NE represents the number of the elements along the side of one-quarter of the full plate (see Fig. 3). For both sample problems, results are given which are obtained using the Gauss quadrature formula with 3 x 3 and 5 x 5 integrating points. Note that the RFSDT, HSDT, and RHSDT finite elements are free of the locking effect, as shown by the results given in Table 9. Tables 10 and 11 give a comparison with the plate element RQ4 developed and tested by the author [35,36]. The RQ4 plate element is a rectangular anisotropic multilayered plate element formulated on the RFSDT plate model with eight degrees-of-freedom at each nodal point, thus giving a total of 32 DOF. The nodal parameters are the same as those for the RFQ40 without i&w and a,,,~. Moreover, an exact integration has been performed in the evaluation of the stiffness, mass and external loads matrices. The numerical results for the element RFQ40 have been obtained with a numerical integration scheme with 3 x 3 Gauss points for the bending and with 5 x 5 Gauss points for the frequencies.
Table 14. Centre deflection 0 under transverse sinusoidal loading and fundamental frequency parameter 6 of a square plate, laminate ORT3LU with o/h = 10, simply-supported on all edges. Effect of the element planform distortion
Ial
0.250
0.275
0.300
0.325
0.350
0.375
Anal.
0.748 -0.75
0.725 -3.84
0.667 -11.6
0.586 -22.2
0.530 -29.1
0.754
% error
0.750 -0.56
ti % error
11.48 0.11
11.49 0.21
11.65 1.67
11.98 4.53
12.35 7.72
12.39 8.11
11.46
XI,
9
A general quadrilateral anisotropic plate element
103
5. CONCLUDING REMARKS Another group of FEM results refer to a comparison with results obtained by other researchers. The In this paper a general quadrilateral plate element results quoted in Table 12 refer to the buckling loads with 40 DOF has been formulated on the basis of a of a four-layer, symmetric (45/-45), angle-ply plate discrete-layer approach that accounts for piecewise simply-supported on all edges under uniform uniaxial compression and shear. The mechanical properties of cubic distribution of the in-plane displacement across each layer are RE = 14; RG = 0.533; RT = 0.323; the plate thickness and allows continuous transverse shear stresses. As remarked, from this element it is v‘T.= 0.30;vTr= 0.55. possible to derive plate elements of the same planThe results used for comparison have been obform based on CPT, FSDT-modified, RFSDT, and tained using a Galerkin-type approach with HSDT. The numerical investigation shows that: R = T = 8, R and T being the number of base functions in x, - and x,-directions used in the approxi1. All the plate models accounting for the transmate solution. The FEM results have been obtained verse shear effects are adequate in predicting the with NE = 6 for the whole plate because of the lack global behaviour of the plate (transverse deflection, of symmetry for the quarter of plate. The Galerkin values for the shear buckling load are given in [54]. natural frequencies, buckling loads, for example). 2. Among the higher-order shear deformation In Table 13 we give the results on the fundamental plate models, the discrete-layer model here named the frequency of a two-layered, antisymmetric cross-ply RHSDT-model appears to be the most accurate in (O/90), square plate (mechanical properties of each predicting the local behaviour of the laminate (distrilayer: RE = 40;RG =0.6;RT= 0.5; NI = 0.25) with bution of the displacements and stresses across the AT = 10, simply supported on the sides x2 = 0, az and thickness, for example). clamped on the edges xL = 0, a,, modelled by the 3. From a FEM point of view, the RHQ40 plate HSDT. The analytical results have been quoted from element appears to be the most computationally cost Reddy and Khdeir [55]. saving with respect to other discrete-layer based plate The numerical results give some insight into the elements because the element number of DOF is effect of the two degrees of freedom (the nodal values independent from the number of layers in the lamiof the curvatures of the reference surface) and the nate. numerical integration. It appears that for coarse 4. The RHQ40 plate element, and the related plate meshes (NE = 1) the accuracy of the RFQ40 is elements CQ40, FQ40, RFQ40, and HQ40 also apgreater than that of the RQ4 element: this is due to pear to work well when used to model an irregular the two additional degrees of freedom, whose benregion (distored planform). eficial effects are greater for coarse meshes where the kinematics of the single finite elements is very important in defining the kinematics of the whole plate. For REFERENCES fine meshes, it appears to be the numerical integration scheme that plays the most important role. 1. C. W. Bert, A critical evaluation of new plate theories The last results, Table 14, give some insight into the applied to laminated composites. Composite Srruct. 2, 329-347 (1984). effect of the distortion of the element planform. In 2. A. K. Noor and W. S. Burton, Assessment of shear this investigation only one quarter of plate has been deformation theories for multilayered composite plates. modelled, where NE = 2 is used. The distortion of the Appl. Me&. Rev. 41, 1-12 (1989). element has been accomplished by moving the point 3. A. K. Noor and W. S. Burton, Assessment of compuP along the diagonal in the direction of the plate tational models for multilayered anisotropic plates. Composite Struct. 14, 233-265 (1990). centre (see Fig. 4). 4. R. K. Kapania and S. Raciti, Recent advances in
analysis of laminated beams and plates. Part I: Shear effects and buckling. AIAA Jnl27, 923-934 (1989). 5. R. K. Kapania and S. Raciti, Recent advances in analysis of laminated beams and plates. Part II: Vibrations and wave propagation. AIAA Jnl27, 935-946
CL CL
I
x Ia1 = 0.375,
L-4
1t 2
+-,/2----c(
IP
x
2,
la,
=
0.375.
Scripta, Washington, DC (1975). 8. J. M. Whitney and N. J. Pagano, Shear deformation in heterogeneous
anisotropic plates. J. Appl. Mech. 37,
1031-1036 (1970). e
x1
Fig. 4. Mesh with the maximum distortion of the element planform x,, /a, = 0.375. CAS 17,1--H
(1989). 6. J. N. Reddy, A review of refined theories of laminated composite plates. Shock and Vibr. Digest 22, 3-13 (1990). 7. R. M. Jones, Mechanics of Composite Materials.
9.
V. Giavotto, Sulla meccanica di pannelli anisotropi ed eterogenei. Istituto Lombard0 (Memorie SC. Mat.) XXV, 437-480
(1969).
10. J. M. Whitney and C. T. Sun, A higher order theory for extensional motion of laminated composites. J. Sound Vibr. 30, 85-97 (1973).
104
M. DI
11. K. H. Lo, R. M. Christensen and E. M. Wu, A higher-order theory of plate deformation. Part II: Laminated plates. J. Appl. Mech. 44, 669676 (1977). 12. M. Levinson, An accurate, simple theory of the statics and dynamics of elastic plates. Mech. Res. Commun. 7, 343-350 (1980). 13. M. V. V. Murphy, An improved transverse shear deformation theory for laminated anisotropic plates. NASA TP-1903 (1981). 14. J. N. Reddy, A simple higher-order theory for laminated composite plates. J. Appl. Mech. 51, 745-752 (1984). 15. A. V. Krishna and S. Vellaichamy, On higher order shear deformation theory of laminated composite panels. Composite Struct. 8, 247-270 (1987). 16. L. Librescu and A. A. Khdeir. Analvsis of svmmetric cross-ply laminated elastic plates using a higher-order theory: Part I-stress and displacement. Composite Srrucr. 9, 189-213 (1988). 17. C. T. Sun and J. M. Whitney, On theories for the dynamic response of laminated plates. AIAA Jnl 11, 178-183 (1973). 18. S. Srinivas, A refined analysis of composite laminates. J. Sound Vibr. 30, 495-507 (1973). 19. S. T. Mau, A refined laminated plate theory. J. Appl. Mech. 40, 606-607 (1973). 20. P. Seide, An improved approximate theory for the bending of laminated plates. Mech. Today 5, 451-466 (1980). 21. M. Di Sciuva, A refined transverse shear deformation theory for multilayered anisotropic plates. Atti Accadet&a Scienze T&o 118,. 279-295 (1984). 22. M. Di Sciuva. Bendina, vibration and buckling. of simply supported thick multilayered orthotropic p&es: an evaluation of a new displacement model. J. Sound Vibr. 105, 425442 (1986). 23. M. Di Sciuva, An improved shear deformation theory for moderately thick multilayered anisotropic shells and plates. J. Appl. Mech. 54, 589-596 (1987). 24. M. Di Sciuva, Further refinement in the transverse shear deformation theory for multilayered composite plates. Atti Accademia Scienze Torino, in press. 25. M. Di Sciuva, Multilayered anisotropic plate models with continuous interlaminar stresses. Composite Struct., in press. 26. H. Murakami, Laminated composite plate theory with improved in-plane responses. J. Appl. Mech. 53, 661-666 (1986). 27. A. Toledano and H. Murakami, A higher-order laminated plate theory with improved in-plane responses. Int. J. Solids Struct. 23, 111-131 (1987). 28. J. N. Reddy, A generalization of two-dimensional theories of laminated composite plates. Commun. Appl. Numer. Meth. 3, 173-180 (1987). 29. J. N. Reddy, On refined computational models of composite laminates. Int. J. Numer. Meth. Engng 27, 361-382 (1989). 30. C. W. Pryor and R. M. Barker, A finite-element analysis including transverse shear effects for applications to laminated plates. AZAA Jnl9, 912-917 (1971). 31. P. Mantegazza and M. Borri, Elementi finiti per l’analisi di pannelli anisotropi. I’Aerotecnica Missili e Spazio 53, 181-191 (1974). 32. A. K. Noor and M. D. Mathers, Shear-flexible finiteelement models of laminated composite plates and shells. NASA TN D-8844 (1975). 33. J. N. Reddy, Free vibration of antisymmetric, angleply laminated plates including transverse shear deformation by finite element method. J. Sound Vibr. 66,565 (1979). 34. P. G. Bergan and X. Wang, Quadrilateral plate bending elements with shear deformations. Comput. Struct. 19, 25-34 (1984).
&WA
35. M. Di Sciuva, Development of an anisotropic, multilayered, shear-deformable rectangular plate element. Comput. Struct. 21, 789-796 (1985). 36. M. Di Sciuva, Evaluation of some multilayered, sheardeformable plate elements. AIAAIASMEIASCEIAHS 26th Structures, Structural Dynamics and Materials Confirence, Orlando, Florida (1985). Also Comput. Struct. 24, 845-854 (1986). 31. J. J. Engblom and 0. 0. Ochoa, Finite element formulation including interlaminar stress calculations. Comput. Struct. 23, 241-249 (1986). 38. A. T. Chen and T. Y. Yang, Static, dynamic and buckling formulation of a symmetrically laminated plate finite element for a microcomputer. J. Composite Mater. 21, 44453 (1987). 39. S. P. Wanthal and T. Y. Yang, Formulation of three simple triangular plane stress anisotropic finite elements for a microcomputer. J. Composite Mater. 22, 106&1079 (1988). 40. P. R. Heyliger and J. N. Reddy, A higher order beam finite element for bending and vibration problems. J. Sound Vibr. 126, 309-326 (1988). 41. A. S. Mawenya and J. D. Davies, Finite element bending analysis of multilayer plates. Int. J. Numer. Mel. Engng 8, 215-225 (1974). 42. S. C. Panda and R. Natarajan, Finite element analysis of laminated composite plates. Int. J. Numer. Meth. Engng 14, 69-79 (1979). 43. M. Epstein and H. P. Huttelmayer, A finite element formulation for multilayered and thick plates. Comput. Struct. 16, &45650 (1983). 44. 0.0. Ochoa, J. J. Engblom and R. Tucker, A study of the effects of kinematics and material characteristics of the fundamental frequency of composites plates. J. Sound Vibr. 101, 141-148 (1985). 45. R. L. Hinrichsen and A. N. Palazotto, Nonlinear finite element analysis of thick composite plates using cubic spline functions. AIAA Jnl 24, 18361842 (1986). 46. R. A. Chauduri and P. Seide, Triangular finite element for analysis of thick laminated plates. Int. J. Numer. Meth. Engng 24, 1203-1224 (1987). 47. D. R. J. Owen and Z. H. Li, Refined analysis of laminated plates by finite element displacement methods--I. Fundamentals of static analysis. Comma. Struct. 26, 907-914 (1987). 48. D. R. J. Owen and Z. H. Li. Refined analvsis of laminated plates by finite element displacement methods-II. Vibration and stability. Comput. Struct. 26, 915-923 (1987). 49. J. N. Reddy, E. J. Barber0 and J. L. Teply, A plate bending element based on a generalized laminate plate theory. Int. J. Numer. Mel. Engng 28, 2275-2292 (1989). 50. E. J. Barbero, 3. N. Reddy and J. L. Teply, An accurate determination of stresses in thick laminates using a generalized plate theory. Inr. J. Numer. Meth. Engng 29, l-14 (1990). 51. K. M. Rao and H.-R. Meyer-Piening, Analysis of thick laminated anisotropic composite plates by the finite element method. Composite Struct. 15, 185-213 (1990). 52. K. Washizu, Variational Methods in Elasticity and Plasticity. Pergamon Press, Oxford (1968). 53. N. J. Pagano, Exact solutions for rectangular bidirectional composites and sandwich plates. J. Composite Mater. 4, 20-34 (1970). 54. J. M. Whitney, The effect of shear deformation on the bending and buckling of anisotropic laminated plates. Composite Struct. 4, 109-121 (1987). 55. J. N. Reddy and A. A. Khdeir, Buckling and vibration of laminated composite plates using various plate theories. AIAA Jnl 27, 1808-1817 (1990).
105
A general quadrilateral anisotropic plate element APPENDIX
Interpolation
functions for u,, u,, g,, and g,
Ni=f(l+&C)(I+~i~), Interpolation
i=1,2,3,4.
functions for w K’=H,(W,(rl), NY’= 4(5)Hb1),
K=H,(I;)H,(tt) W = H, tTW&)
NYC= H&)H,(‘/), W = H&)Wf),
W = WS)H,ttl) W = H,(S)H,(rl)
W” = H,(WMrtX N;” = H&W&),
N,WI= H,(S)H&) N$ = H,WH,(I)
W~=~&)H,(~), W = H&)H;l(tl),
W’=H&W,(v) N4”((= WWMlj)
H,(r)=&(-3rS+tOr3-15r+8) Hz(r) = h(3rs-
10r3 + 15r + 8)
H,(r)=&{-3rS+r4+10r3-6r*-7r+S) Hd(r) = i(-3rs
- r4 -+ IOr + 6r2 - 7r - 5)
H,(r)=fg(-r5+P+2r3--2r2-r
+ 1)
He(r) = h(r’ f r4 - 2r3 - 2r2 + r + 1).
PrintedinGreat Britain.
aw-7949/93 s6.00+0.00 1993 Rrgamonl'ms Ltd
A GENERAL QUADRILATERAL MULTILAYERED PLATE ELEMENT WITH CONTINUOUS INTERLAMINAR STRESSES-f M. DI SCIUVA Department of Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruxxi 24, 10129 Turin, Italy (Received
16 December
1991)
Abstract-The purpose of this work is the formulation of a general quadrilateral multilayered anisotropic plate element. The plate element is formulated on the basis of a refined third-order shear deformation plate theory recently proposed by the author. The plate model makes use of a displacement field that fulfils a priori the geometric and stress continuity conditions at the interfaces between the layers, and requires only five generalized displacements to describe the kinematics of the laminate deformation. To test the performance of the developed plate element, both closed-form and finite element-based solutions are given and compared with results from three-dimensional elasticity and other approximate two-dimensional plate models.
linear or nonlinear distribution of the in-plane displacements or the transverse shearing stresses in the thickness direction. These approaches consider all layers as one equivalent single anisotropic layer. CPT, based on neglecting transverse shear strains and transverse normal strains and stress in the laminate [7], and the first-order shear deformation theory (FSDT), in which the transverse shear strains are assumed to be constant through the thickness [8] and the transverse normal strain and stress are neglected, are well-known examples of plate models with assumed global linear distribution. The so-called higher-order shear deformation theories (HSDT) are examples of plate models with assumed global nonlinear distribution [g-16]. This assumption generally results in incompatible shearing stresses between every two adjacent layers. In fact, from the theory of elasticity we know that the displacements and stresses at the interface k between the kth and (k + 1)th bonded layers must satisfy the following contact conditions$
1. INTRODUCTION The use of fibre reinforced laminates in aerospace, automotive, shipbuilding and other industries has increased tremendously in the past several years. This
is largely due to the high strength-to-weight ratios of composites as well as their ability to be tailored to meet design requirements of strength and stiffness. Coinciding with these new applications is the interest in the accurate prediction of the detailed response and failure characteristics of laminated anisotropic plates. The anisotropy and nonhomogeneity and larger ratio of longitudinal to transverse moduli of these new materials demand improvement in the classical plate theory (CPT), plate theory based on the well-known Kirchhoff assumptions. As a result, modelling approaches which take into account transverse shear deformation have been the topic of serious research in the last two decades. It suffices here to mention the survey papers of Bert [1], Noor and Burton [2,3], Kapania and Raciti [4, 51, and Reddy [6]. Some of these modelling approaches are extensions of similar approaches used for isotropic plates. Belonging to this category are the so-called smeared laminate models, that is, those models based on an assumed global (i.e. one for the whole laminate)
geometric continuity
conditions:
ri,(z;) = ti,(z$) stress continuity &,(z;)
t This paper represents an updated and completed version of the one presented at the European Conference on New Advances in Computational Mechanics, Giens (Var)-France (April, 1991). $ In this paper, if not otherwise stated, the Einsteinian summation convention over repeated subscripts will be adopted with Latin indices ranging from 1 to 3 and Greek indices from 1 to 2.
conditions:
= e&:
);
G(ZL ) = &,(zk+ ).
In these relations, zii (i = 1-3) are the displacement components along an orthogonal frame xi on the reference surface Q, xj = z being the normal to R. The minus and plus refer to the values of the attached functions for z = zk - 0 and z = z, + 0, 91
92
M. DI SCIUVA
respectively, where z, is the value of the z-coordinate at the interface k between the layers k and k + 1 (Fig. 1). From the geometric continuity conditions follow the strain continuity conditions &,a@;) = &,a@: ). Continuity conditions cannot be imposed on the remaining stresses and strains; in general, these are discontinuous at the interfaces. If the transverse shearing strains are assumed to be continuous functions of the thickness coordinate z, it is not difficult to realize that the contact conditions on the transverse shearing stresses cannot be satisfied if the layers have different mechanical properties. In the literature there are a number of approaches in which the contact conditions are satisfied. These models are based on an assumed piecewise linear or nonlinear displacement or stress distribution in the thickness direction and are known as discrete-layer models [17-291. Most of these models consider the single layer as a plate and impose the contact conditions as constraint conditions on the transverse stresses at the interfaces. The major drawback is that the number of governing equations increases with the number of layers in the laminate. The approach proposed by the author [21-251 overcomes this drawback as it makes use of an assumed displacement field that fulfils a priori the geometric and stress continuity conditions at the interfaces between the layers, and requires only five generalized displacements to describe the kinematics of the laminate deformation, as in FSDT. In its original version, a piecewise linear variation of the in-plane displacement (the so-called zig-zag model [21-231) has been assumed. This model may be considered the counterpart of the FSDT in the class of the discrete-layers models, and as with the FSDT, it is unable to fulfil the conditions of zero transverse shear stresses on the bounding surfaces of the plate. To overcome this restriction, the author [24,25] proposed an improvement of the zig-zag model (the so-called refined higher-order shear deformation theory, RHSDT model), by assuming a displacement field which allows a nonlinear variation of the inplane displacements through the laminate thickness and fulfils a priori the geometric and stress continuity conditions at the interfaces and the static condition of zero transverse shear stresses on the top and bottom surfaces, at least for symmetric laminates. Moreover, only five generalized displacements are required to describe the kinematics of the laminate, as in the FSDT and in the zig-zag model (in the following called the RFSDT). In the last two decades, there has been a considerable effort towards the development of shear deformable, multilayered anisotropic plate elements. This
development follows the same path as in the analytical modelling approaches. In effect, as a basis for the formulation of plate elements of this type, many investigators have employed smeared laminate models [30-40], while others make use of discretelayer models [35,36,41-511. In the first approach, the multilayered anisotropic plate element is reduced to an equivalent single anisotropic layer element. Many of the finite element formulations belonging to the second group make use at least of one solid or plate element per layer in the thickness direction. Then, this approach requires numerous unknowns and, as a consequence, these formulations become nontractable in the solution of practical engineering problems. The rectangular plate element formulated by the author [35,36] on the basis of the RFSDT model, zig-zag model approach, overcomes this drawback because the number of nodal parameters is independent of the number of layers. Based on the RHSDT model [24,25], in this paper we formulate a general quadrilateral multilayered anisotropic plate element. To test the performance of the developed plate element, both closed-form and finite elementbased solutions are given and compared with results from three-dimensional elasticity and other approximate two-dimensional plate models. The paper is structured as follows: based on the assumed piecewise cubic displacement field and on the von Karman nonlinear strain-displacement relations, in Sec. 2 we derive the nonlinear equations of motion and the static stability equations of multilayered anisotropic plates by making use of the dynamic version of the principle of virtual work. In Sec. 3 we formulate the quadrilateral plate element. To give an insight on the accuracy and reliability of the developed finite element, Sec. 4 contains analytical and FEM results on bending under sinusoidal transverse loading, undamped frequencies and buckling loads under uniform axial compression. Whenever possible, a comparison with results from three-dimensional elasticity and other approximate two-dimensional models is also given. 2. ANALYTICAL
FORMULATION
2.1. Kinematics Let us consider a multilayered plate of total thickness h and reference surface R composed of N orthotropic layers of uniform thickness perfectly bonded together. The principal axes of elasticity of each layer can be arbitrarily oriented with respect to the laminate orthogonal reference frame xi, where x, is an orthogonal frame on the reference surface Q (here chosen to be the midsurface), and xj = z is the normal to R (see Fig. 1). The material properties and the thickness of each layer may be entirely different. The discrete-layer plate model considered here is based on the following piecewise cubic
A general quadrilateral anisotropic plate element
93
Fig. 1. Geometrical laminate configuration. approximation of the in-plane displacement the plate thickness
across
r&(xa, 2; t) = u.(xs; t) - 2 &w(xfl; t) +f,&)g&;
t)
r&(xfl, 2; t) = w(xg; t) or, in matrix notation ii = Z,d, where I’= [u”, r& r&l, d: = [d; dj, d;] d; = b,
4,
d;=[w
S;],
@:,=[a,~
determined; z, are the coordinates of the N - 1 interfaces; Yk is the Heaviside unit function (it has the value 0 for z < z,, and 1 for z 3 zk). The tracers 6, and 6, identify the contributions through (dF= 1; CRT (6, = 8r = U.jk = 0), FSDT-modified 6, = aDlP= 0), HSDT (6, = 8r = 1; a@ = 0), RFSDT (a,= 1; 6,= 0) RHSDT (S,= 6, = 1) plate models. The notation a,( .) represents the partial derivative a(.)/aa. Note that the displacement components li, are continuous functions of the thickness coordinate z for all values of the parameters aplr . Thus, the geometric continuity conditions are automatically satisfied and we must only satisfy the stress continuity conditions on u,,. By satisfying these stress continuity conditions, it may be shown [24,25] that the following expressions for a,, hold?
a2w] k-l
k-l
%lk
al
d;=k,
0
0
--z
[ 001
0 0
%lq
q-1
1 0 0 --z s2”,=01
x ullq 1 +AIU, 1 q=l
tal)k+
=Allk
0
.flI A* fil122. 0
0I
k-l
auk --&k
caJ)k
+
taf)k
+
caf)k
+
x q=I
hq
k-l
N-l
u22k =
A22k
k=l
with
~22,
k-l
The functions fur(z) specify the distribution of the transverse shear strains along the thickness. In this study we make use of the following expression [24,25]
Am(z) =W@) + c %k(Z-%)Yll
k-l
1 q=l
1 q=I
U22q
1
+&k
1 1
c q”
k-l
+
&k
+
A2,k
c q-1
Ilk
x q-1
&%d?::k+,-Qt:,Q::,+, _, gk
Here u.(xP; t) and w(x#; t) are the displacements in the x, and z directions, respectively, of a point belonging to R; g,(xs; t) are the transverse shear rotations for z = 0; Sf is the Kronecker delta (St = 1 if a = /?; St = 0 if a #/I); a,, are parameters to be
A
I2k
=
Q::kQ::k+,
-
Q::kQ::k+,
Bk
A2x=Q::kQ:k+,
-Q::kQf:k+,_, 9k
t Obviously, the expressions for u~,,~are related to the assumed expressions for the strains. Those given here hold for the strain-displacement relations given in Sec. 2.2 and for the lamina constitutive relations given in this section.
U21q
k-l
where A
%q
st,=(Qt:k+,)‘-Ql:+,Q::k+,.
UI~,
94
M. DI SCIWA
Thus, for a specified laminate, uegk are known constants only depending on the transverse shear mechanical properties of the constituent layers. In the above relations, the notation (afit represents a,j valued at z = z,. Moreover, Q$ and Qz& are the so-called reduced stiffnesses of the k th orthotropic layer in the plate coordinate system (x,, z). They relate the stress components d to the strain components P (lamina constitutive relations)
respectively. Notice that we have used c for the linear portion of the strain. We can write i2=E+iiNL with 2ENL= A W’ 8 IV
C = %“,dc,
where -zU
or, in matrix notation
F
0
0
dF
1
and
a=Qg, where
I Qft
Q=
Q::
0
Qt:
0 0 0 Q::
Q:: Q:: Q:: 0 Qi: Q:: Q!: 0 0 0 0 Qt: 0 0 0 Q:: Q::
Q is the (5 x 5) matrix of the transformed stiffnesses of the generic layer. 2.2. Strain-displacement
e; = A;=
reduced
0
a2
0
0
a,
0
a,
a,w 0 0
azw
1
1.
a,w
0
0
a,w
0
0
Moreover. let
e;000 0e; I1 Qf= 0e; 00 d, = 53:ds
relations
In deriving the static stability equations, we make use of the von Karmin strain-displacement relations (moderately large deflections and small strains)
[
a,
with
i
and
where
where d, is the vector of the generalized displacements d: = [d;
w
d;].
and 26a8=a,u,+a,u,,
k,,=
-asBw
where L,~ and k,, are the (membrane) extensional strains and the curvatures of the reference surface,
2.3. Governing equations
In the present study we are interested in deriving the linear equations of motion and the static stability equations. Both the problems can be approached by making use of the principle of virtual work 6# = 6L’+
aLin,
95
A general quadrilateral anisotropic plate element where 6@ is the virtual variation of the strain energy; 6L’ is the virtual variation of the work done by the applied external forces; 6L” is the virtual variation of the work done by the inertial forces.
where &I+, =
JR
SG2=
JR
2.3.1. Strain energy. Substituting the expressions for the von Kbrmin strains in the expression of the virtual variation of the strain energy 6@ =
(a-‘~%) dR, s Cl
where B is the stress state corresponding to a generic time of motion and 60 the strain state related to the virtual displacements Slii, we obtain
S@ =
Jn
W’6d,dQ+
d:Xo Sdc di2 =
(97;d,)‘Xo(ka;
ad,) dR
8 :N 68, dQ.
In terms of force and moment stress resultants the explicit expression for the virtual variation of the strain energy of the multilayered anisotropic plate reads
Jn
@;N&‘,dfl
with
+
&N&ia,4 6w+ (ass,,- T,) k,l dfi
+
Jr
w = (%2”:d) being the vector of the force and moment resultants for unit length
stress
~‘=[NI,;N,~;N~,;N,,;M,,;M,,;M~,;M~~; S,,;
Jfl
[Nm6u, - Mg, da,
W +
S, Sg,
+ (v, + a,M,,, + N,,,) 6w] 6r, where V, = n, a,M,,
and
Sn; S,,; 8,; T,; Tzl = IN’;M’;S’;T’l,
where ( . )n, =
(N,,, Kb, S,,) = K&)(1>
z&r)>
T, =
hi2 (...)= J’ J (...)&=;
Here, and in the following, t stands for transpose and
4
( . )2n -
n2(.
)*n
on r, r denoting the boundary of R, and n, the direction cosines of the outward unit normal n to r with respect to the a-axis and t the tangent unit vector. Let S,,(s) be the 2.3.2. Work of external forces. portion of the lateral cylindrical surface of the sth layer on which the external loadings per unit area A(s) are prescribed; let /jZ =J%(x,, -h/2) +pZ(x,, +/z/2) be the transverse loading per unit area. Under these conditions, it is found that
(...)dz.
-h/2
s=l
2%- I
Moreover, by taking into account the constitutive equations, we obtain W = .%,,d, + +f-, A,,@,.,
where we have introduced the following resultants and moment resultants of the external forces
where XO =
W’:Q~‘,),
x, = (S:Q>.
Then 2.3.3. Work of internaI forces. By using the approximate expressions for Sri, (as well as for &), an
M. Dr !%XJVA
96
approximate expression for the virtual variation the work of inertial forces is obtained as
of
where
equal to zero, the following nonlinear motion result
@g,)
h=(s:P~“) and 53: is a matrix of differential that
operators
such
d, = B:d,. We have
a,S,,-
T,=m&,-m~asti
equations of
+m$&.
Boundary conditions : Let ru be the portion of r on which are prescribed the values of the generalized coordinates and let r, be that portion of r on which are prescribed external tractions. Thus the set of boundary conditions reads Prescribed on r.
Natural on t’,
u, = ii w= ;
N,=Nm v,+~,M,,+N,,=
+m,ii,-m,d,G M”=iV” s,=S,.
w,n = C’,” g.=&
v=+a,li;i,,
+nsm@+
2.4. Static stability equations 0000010 0100001 In expanded form and after performing grations by parts, the result is
[(mot&- m, a,$ + mf&)
6L” = -
some inte-
6u,
sR
where we have introduced the following inertia resultants (i =O, 1,2) (m,; m?; m$) = (P(z’; z%~;_L_t$,s)>. In the previous relations, p denotes the material mass density, and the overdot indicates differentiation with respect to time. Equations of motion: By substituting the expressions for SQ, 6L’, and 6L” into the principle of virtual work, rearranging the various contributions and setting the coefficients of the virtual variations of the generalized coordinates in the integral over R
The static stability equations of the laminated plates subjected to combined in-plane normal and shear loads may be formulated by using the Euler method of the adjacent equilibrium configurations [52]. In this approach, we assume the existence of at least one additional, distinct equilibrium configuration (the adjacent equilibrium configuration) in the very close neighbourhood of the original equilibrium configuration (the equilibrium configuration just prior to the occurrence of the buckling). If such an adjacent equilibrium configuration exists, the body may change suddenly from the original equilibrium configuration (here characterized by ui, w”, etc.) to the other under the stimulus of infinitesimal disturbance u,, w, etc. Although these terms have been used to denote the displacement field in the previous derivation of nonlinear equations of motion, it is clear that in the present stability equations they assume a different significance. More precisely, in the nonlinear equations of motion u,, w are the displacements measured from the reference stress-free state, while here they are measured from original equilibrium configuration. Using the expressions for the nonlinear strains as given by von Karrnbn, on the assumption that the external forces and moments in both equilibrium configuration are identical, the principle of virtual work of the adjacent equilibrium configuration reads
s
((a-“’ W’ + d” 6’))
n
dR = 0,
A general quadrilateral anisotropic plate element where the prime and double prime represent the quantities which contain terms that are linear or quadratic in the derivatives of the infinitesimal disturbance displacements P’=~+~(A,&,+A,8~),
97
and
2Z”=A,8,.
Noting that
Let us assume a linearized form for the relations between the incremental stresses and strains, i.e. d’ = QG’, then
A,,&, = A,@@ yields C=i+A,@,. By taking the virtual variations with respect to the disturbance quantities, yields
W’=A,68,.
where X0 and X, have been defined previously, and X2 = (Q). In explicit form
Substituting the expressions for the strains in the expression of the principle of virtual work for the adjacent state, we obtain
+ where &Do=
(,‘I
Si) dR
SW’= -
+
(tf”.2Z’~) Sd, dR
sn =
s n
8”6d EdR
sr
n, $,(N:swo) dT
s s
(N:pa,pw +d,N~fld,w)Sw
di-2
0
sn =
(N& i&w0 + d,N;, 8,~‘) SWd&2
6@’ = -
n, as(N;, w) dT.
=
The expression for S@Ois formally the same as that given in Sec. 2.3.1. 3. THE FINITE ELEMENT FORMULATION
This section deals with the displacement finite element formulation of the multilayered plate models summarized in the previous section. Let the region of the plate be divided into a number of quadrilateral elements with straight sides (Fig. 2).
rl’ (-1.+1)
@
* 5
(-1.4)
0
Fig. 2. Mapping of the geometry.
@
(+1.-l)
M. DI
98
.%JlJVA
Table 1. Centre deflection I? of a square plate, laminate ORT3LU, simply-supported on all edges under transverse s&so&l loading. I+,-~, = 0.43 1
Table 2. Centre deflection $ of a rectangular plate (a2 = 3a,), laminate ORT3LU, simply-supported on all edges under transverse sinusoidal loading. J,, = 0.503
AT
FSDT
HSDT
RFSDT
RHSDT
AT
FSDT
HSDT
RFSDT
RHSDT
Ref. [53]
4 10 20 50 100
1.776 0.669 0.492 0.441 0.434
1.922 0.713 0.504 0.443 0.434
1.891 0.723 0.508 0.444 0.434
1.972 0.754 0.517 0.445 0.435
4 10 20 50 100
2.363 0.803 0.578 0.515 0.506
2.641 0.862 0.594 0.518 0.507
2.717 0.881 0.599 0.519 0.507
2.757 0.919 0.610 0.520 0.508
2.820 0.919 0.610 0.520 0.508
The four-node quadrilateral plate element formulated on the basis of the RHSDT has 10 DOF at each nodal point, thus giving a total of 40 DOF. In the following, this plate element will be named RHQ40. The nodal parameters are the five generalized displacements ul, w, g,; the two total rotations 19,w; the two curvatures &,w and the twist a,w of the reference surface a. Notice that following the scheme of Sec. 2, it is possible to derive from the RHQ40 plate element the corresponding plate elements formulated on the CPT (here named CQ40), on the FSDT-modified (here named FQ40), on the RFSDT (here named RFQ40), and on the HSDT (here named HQ40). 3.1. Shupe functions To make the computation of the involved integrals uniform, a topological transformation (mapping) of the geometry element from the physical plane x, , x2 to the transformed plane 5, r) (Fig. 2) has been accomplished using the Lagrange linear interpolation functions
principle) for the plate mode1 considered here (CPT, FSDT-modified, HSDT, RFSDT, and RHSDT) requires C’-approximation of the transverse deflection w, and Co-approximation of the in-plane displacement w, and the shear rotations g,. Thus, Hermite cubic or higher-order interpolation functions should be used for approximating w, and linear or higherorder Lagrange interpolation functions for approximating u, and g,. In formulating the plate element RHQ40, the Lagrangian linear interpolation function Li has been used for approximating d,, and d8; Hermite sixth-order polynomials Hi for approximating w. Thus, an isoparametric formulation has been used for d, and d,, and a subparametric formulation for w. Note that the element is able to model a constant strain state and a rigid body displacement field in the membrane and transverse behaviour. Moreover, the element is conforming because of the continuity of the displacements and the slopes across adjacent elements. Let A: = [A: A;,. A;],
where .Xi = a (1 + tit) (1 + qiq), x,.~being the coordinates of the nodal point i in the physical plane and &, vi those in the transformed plane. This transformation maps a generic quadrilateral region with straight sides in the physical plane into a square region of boundary [- 1 < 5,~ < 11 in the natural plane. The finite element formulations based on the displacement approach (i.e. on total potential energy Table 3. Maximum
with A:= [ul
41, A; = k’, i&l,
Al. = [w’
alwl
T 0.437 TO.734 TO.620 TO.845 -0.755 +0.801
a12d
azzwq
be the vector of the physical degrees of freedom surrounding element e. Each of the 10 subvectors has
stresses in a square plate, laminate ORT3LU, transverse sinusoidal loading
FSDT HSDT RFSDT RHSDT Ref. [53]
a2Wl a,,wl
simply-supported
TO.477 TO.503 TO.661 TO.598 -0.556 +0.534
TO.037 TO.050 + 0.045 TO.051 -0.0511 +0.0505
+0.256
0.217
+0.381 +0.369 +0.364 +0.357 +0.357
+0.111 +0.117 +0.123 -0.124 +0.1228
+0.395
0.082
FSDT HSDT RFSDT RHSDT Ref. [53]
TO.513 TO.568 TO.549 TO.595 -0.590 f0.590
TO.254 TO.269 T 0.289 TO.292 -0.288 +0.285
70.025 T 0.028 TO.027 7 0.029 -2.89 +2.89
CPT
TO.539
TO.180
TO.021
0.037 0.286 0.262 0.241
on all edges under
0.197 0.209 0.254 0.230
A general quadrilateral anisotropic plate element
99
Table 4. Maximum stresses in a rectangular plate (a2 = h,), laminate ORT3LU, simply-supported on all edges under transverse sinusoidal loading AT
4
10
elements sponding
Model FSDT HSDT RFSDT RHSDT Ref. [53]
011 TO.613 T I .036 TO.925 f 1.205 -1.10 +1.14
512
u22
F 0.093 TO.103 TO.114 TO.108 -0.119 +0.109
T 0.020 TO.026 F 0.026 T 0.028 - 0.0269 f0.0281
013
93
0.436 0.383 0.366 0.329
0.028 0.030 0.032 0.030
+0.351
+0.0334 +0.014 +0.015 +0.015 +0.015 +0.0152 f0.0108
FSDT HSDT RFSDT RHSDT Ref. [53]
TO.621 TO.692 T 0.675 TO.730 -0.725 +0.726
TO.037 F 0.040 T-o.041 TO.042 -0.0435 +0.0418
TO.010 TO.012 TO.012 TO.012 -0.0123 +0.0120
+ 0.439 +0.430 + 0.427 +0.419 + 0.420
CPT
f 0.623
TO.0252
T 0.0083
+0.440
which are the nodal values variable; for example II; = [uf
u:
u;
of the corre-
and
“y;=X*=
uf].
where the submatrices 1x4
N-0 [ 0
.N
1
are row matrices
of dimension
he the appropriate matrices of the chosen shape functions in terms of the natural coordinates 5, r~such that the following formal relations hold for the element e =A’&, = [Np], with i = l-4. The explicit expressions for the elements of these matrices are given in the Appendix. 3.2. Element elastic st$ness matrix The element elastic stiffness matrix K, is obtained from the linear part 69’ of the strain energy 6@ (see Sec. 2.3.1). Taking into account that in the present finite element formulation
Then
where yields
Table 5. Frequency parameter of the fundamental bending mode cG of a square plate, laminate ORT3LU, simplysupported on all edges AT
4 10 20 30 50 100
CF’T
FSDT
HSDT
RFSDT
RHSDT
14.501 15.104 15.197 15.214 15.233 15.227
7.413 12.163 14.230 14.757 15.053 15.183
7.116 11.790 14.060 14.670 15.019 15.175
7.163 11.700 14.012 14.645 15.009 15.172
7.027 Il.462 13.888 14.579 14.983 15.165
Table 6. Asymptotic buckling stress parameter K, of a square plate, laminate ORT3LU, simply-supported on all edges under uniaxial compression fl,, AT
4 10 20 30 50 100
CPT
FSDT
HSDT
RFSDT
RHSDT
23.495 23.495 23.495 23.495 23.495 23.495
5.706 15.138 20.588 22.101 22.972 23.362
5.272 14.221 20.099 21.843 22.870 23.362
5.358 14.008 19.961 21.767 22.839 23.327
5.138 13.436 19.608 21.572 22.760 23.306
M. DI SCIWA
loo
Table 7. Maximum normalized deflection iir for a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges, under transverse sinusoidal load. Numerical integration scheme with 3 x 3 and 5 x 5 Gauss points. Convergence analysis
NE
3x3
1 2 3 4 5 Anal.
3x3
5x5
3x3
5x5
0.690 0.709 0.711 0.712 0.712
0.675 0.705 0.710 0.711 0.711
0.726 0.750 0.752 0.753 0.754
0.708 0.745 0.750 0.752 0.753
5X5
0.700 0.720 0.722 0.723 0.723
RHQ40
HQ40
RFQ40
0.684 0.716 0.720 0.722 0.722
0.713
0.723
0.754
Table 8. Fundamental frequency parameter CT,for a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges. Numerical integration scheme with 3 x 3 and 5 x 5 Gauss points. Convergence analysis HQ40
RFQ40
RHQ40
NE
3x3
5X5
3x3
5X5
3x3
5X5
1 2 3 4 5 Anal.
11.848 11.695 11.680 11.675 11.673
11.982 11.727 11.694 11.683 11.678
11.990 11.839 11.825 11.820 11.818
12.116 11.870 11.838 11.827 11.823
11.644 11.475 11.458 11.453 11.450
11.793 11.511 11.474 11.462 11.456
11.700
11.463
11.790
where
From the previous positions, it is readily found that
B, = %r”9;J(r is the strain matrix relating the strain components of the DOF in the physical plane. In the above expression, J is the Jacobian of the transformation from the natural to the physical plane given previously, and the matrix JT relates the partial derivatives of the shape functions with respect to the natural coordinates to the partial derivatives of the shape functions with respect to the physical coordinates
where
is the matrix which relates the slopes reference surface to the nodal degrees the transverse displacement w. Upon the expression for 8, into the virtual
of the element of freedom for substitution of variation of the
strain energy 6@“, yields
where the definitions follow from
for the submatrices
J,,
etc.
and “0, , ‘VI2have the same form as 8,) 0, by replacing 1 with e, and 2 with 1. The same holds for “PS. Then I
I
ss-I
-I
K,=
B;X-,B,J
dt dq.
3.3. Element geometric st@ness matrix As is well known, the element geometric stiffness matrix Kge is obtained from the expression of the virtual variation of the strain energy due to the nonlinear terms in the incremental strains S@: (Sec. 2.4).
/.-.,2-wj
Fig. 3. Mesh with NE = 5.
*I -
A general quadrilateral anisotropic plate element
101
Table 9. Centre deflection I? of a square plate, laminate ORT3LU, simply-supported on all edges under transverse sinusoidal loading. Effect of the thickness ratio on the discretization error AT
FQ40
(Anal.)
HQ40
(Anal.)
RFQ40
(Anal.)
RHQ40
(Anal.)
4 10 20 100
1.566 0.630 0.482 0.433
(1.568) (0.63 1) (0.521) (0.433)
1.917 0.712 0.504 0.434
(1.922) (0.713) (0.504) (0.434)
1.887 0.723 0.507 0.434
(1.891) (0.723) (0.508) (0.434)
1.966 0.754 0.517 0.435
(1.972) (0.754) (0.517) (0.435)
from which follows the expression geometric stiffness matrix
In order to show the accuracy and reliability of the RHSDT model and to assess the convergence characteristics of the corresponding finite element formulation (the RHQ40 finite element), some numerical investigation has been performed. In the following, analytical and FEM results for the following three sample problems are given. 1. Bending under sinusoidal transverse loading
where
(%AQ’N”@N,)J
K,=
dr d+
The consistent element mass matrix is readily obtained from the expression for 6L” by taking into account that
/!?A,
where m and II are the number of halfwaves in the x, and x,-directions, respectively; a, and aI, plate dimen-
sions in the x, - and x2-directions, respectively; and, A, maximum intensity of the distributed loading. 2. Undamped frequencies.
00 f,u= 00I9-1 001 d, = duCf%.K)A,,
3. Buckling under uniform in-plane axial compression and shear of rectangular and square laminated plates. For purposes of comparison, the evaluation
I
i
a,=:,
pz =pO sin u,x, sin j$xz,
3.4. Consistent element mass matrix
where
4. NUMERICAL RESULTS
for the element
and “9: is the same as $’ with 1 replaced by < and 2 replaced by r,t in the derivatives. Moreover, f-’
denotes the inverse of the Jacobian matrix. The result is
has been conducted on sample problems for which exact elasticity solutions (ES) are available. In addition, numerical results from CRT, FSDT with shear correction factor x = 5/6, RFSDT, and HSDT, are also given. Further numerical results showing the accuracy of the proposed approach can be found in [25]. 4.1. Analytical results
Explicit expressions for the elements of the elastic stiffness, geometric stiffness and mass matrices are not given here. It is sufficient to note that the formulation upon which the expression for the element mass matrix rests also includes the in-plane and rotatory inertia terms.
The numerical results here given refer to a threelayered, symmetric cross-ply plate simply supported on all edges for which the exact elasticity solution has been given [53]. The layers have equal thickness with fibres in the outer layers oriented in the x,-direction (0 degree) and those in the inner layer oriented in the x,-direction (90 degrees). The mechanical properties of each layer are RE = EL/ET= 25; RG = GLT/ Er = 0.5; RT = G&E, = 0.2; vLr = 0.25; here L and T are the directions parallel and normal to the fibres, respectively, and vLr is the Poisson’s ratio measuring transverse strain under uniaxial stress parallel to the fibres. In the following, this laminate
Table 10. Centre deflection ii~of a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges under transverse sinusoidal loading. Comparison between RQ4 and RFQ40 elements
Table 11. Fundamental frequency parameter (I, for a square plate, laminate ORT3LU with AT = 10, simply-supported on all edges. Comparison of the convergence characteristics of the RQ4 and RFQ40 elements
where M,=
NE R:$O Anal.
C’~,“.NYd%WdfwWi’4Jdt dv
1
2
3
4
NE
1
2
3
4
0.606 0.700
0.722 0.720 0.723
0.723 0.722
0.723 0.723
RQ4 RFQ40 Anal.
12.938 11.982
11.836 11.727 11.700
11.734 11.694
11.697 11.683
M. DI
102
Table 12. Transverse deflection and buckling loads for a four layers, symmetric (45/-45), angle-ply square plate (AT = lo), simply-supported on all edges (NE = 6 for the whole plate)
FQM RHQ40 Galerkin
I?
N,,,.
8.208 8.439 8.090
18.535 18.110 19.230
SCIUVA
Table 13. Frequency parameter of the fundamental bending
mode (I, for a two layer, antisymmetric (O/90) cross-ply square plate (AT = 10). simply-supported on two opposite edges and clamned on the other two sides NE
R, ,“, 23.133 22.129 23
-60.808 - 54.704 -59
Ref. 1541
will be called ORT3LU. As an example of the accuracy of the proposed RHSDT model, Tables 1 and 2 give results of the maximum deflection
3x3
1 2 3 4 5 Anal. Ref. 1551
of the fundamental load parameter
5x5
29.53854 17.41517 16.48013 16.17408 16.02959 15.709
30.95164 18.44556 17.03638 16.54957 16.31439 15.709
bending mode and the buckling
E,h’ 3 = lOOA,~ Peal
as a function of the plate thickness ratio AT = q/h, for the following simple supported boundary conditions XI
x,=O,a,:
4.2. FEM results
N*,=M,,=S,,=u*=w=g~=0
0, a, :
=
N22=M22=S22=u,=w=g,=0.
For the same plate, Tables normalized maximum stresses
3 and 4 give the
%3 _ *23 =
DoAT
Notice the high accuracy of the proposed approach. Also notice that the results from RFSDT, HSDT, and RHSDT are obtained without introducing any arbitrary shear correction factor.
’
Tables 5 and 6 give results for the undamped quency
fre-
In order to show the convergence characteristics of the RHQ40 finite element, Tables 7 and 8 summarize the FEM results for the centre deflection 5 and the undamped frequency parameter @.Due to symmetry, only one quarter of the plate has been modelled. In the tables, NE represents the number of the elements along the side of one-quarter of the full plate (see Fig. 3). For both sample problems, results are given which are obtained using the Gauss quadrature formula with 3 x 3 and 5 x 5 integrating points. Note that the RFSDT, HSDT, and RHSDT finite elements are free of the locking effect, as shown by the results given in Table 9. Tables 10 and 11 give a comparison with the plate element RQ4 developed and tested by the author [35,36]. The RQ4 plate element is a rectangular anisotropic multilayered plate element formulated on the RFSDT plate model with eight degrees-of-freedom at each nodal point, thus giving a total of 32 DOF. The nodal parameters are the same as those for the RFQ40 without i&w and a,,,~. Moreover, an exact integration has been performed in the evaluation of the stiffness, mass and external loads matrices. The numerical results for the element RFQ40 have been obtained with a numerical integration scheme with 3 x 3 Gauss points for the bending and with 5 x 5 Gauss points for the frequencies.
Table 14. Centre deflection 0 under transverse sinusoidal loading and fundamental frequency parameter 6 of a square plate, laminate ORT3LU with o/h = 10, simply-supported on all edges. Effect of the element planform distortion
Ial
0.250
0.275
0.300
0.325
0.350
0.375
Anal.
0.748 -0.75
0.725 -3.84
0.667 -11.6
0.586 -22.2
0.530 -29.1
0.754
% error
0.750 -0.56
ti % error
11.48 0.11
11.49 0.21
11.65 1.67
11.98 4.53
12.35 7.72
12.39 8.11
11.46
XI,
9
A general quadrilateral anisotropic plate element
103
5. CONCLUDING REMARKS Another group of FEM results refer to a comparison with results obtained by other researchers. The In this paper a general quadrilateral plate element results quoted in Table 12 refer to the buckling loads with 40 DOF has been formulated on the basis of a of a four-layer, symmetric (45/-45), angle-ply plate discrete-layer approach that accounts for piecewise simply-supported on all edges under uniform uniaxial compression and shear. The mechanical properties of cubic distribution of the in-plane displacement across each layer are RE = 14; RG = 0.533; RT = 0.323; the plate thickness and allows continuous transverse shear stresses. As remarked, from this element it is v‘T.= 0.30;vTr= 0.55. possible to derive plate elements of the same planThe results used for comparison have been obform based on CPT, FSDT-modified, RFSDT, and tained using a Galerkin-type approach with HSDT. The numerical investigation shows that: R = T = 8, R and T being the number of base functions in x, - and x,-directions used in the approxi1. All the plate models accounting for the transmate solution. The FEM results have been obtained verse shear effects are adequate in predicting the with NE = 6 for the whole plate because of the lack global behaviour of the plate (transverse deflection, of symmetry for the quarter of plate. The Galerkin values for the shear buckling load are given in [54]. natural frequencies, buckling loads, for example). 2. Among the higher-order shear deformation In Table 13 we give the results on the fundamental plate models, the discrete-layer model here named the frequency of a two-layered, antisymmetric cross-ply RHSDT-model appears to be the most accurate in (O/90), square plate (mechanical properties of each predicting the local behaviour of the laminate (distrilayer: RE = 40;RG =0.6;RT= 0.5; NI = 0.25) with bution of the displacements and stresses across the AT = 10, simply supported on the sides x2 = 0, az and thickness, for example). clamped on the edges xL = 0, a,, modelled by the 3. From a FEM point of view, the RHQ40 plate HSDT. The analytical results have been quoted from element appears to be the most computationally cost Reddy and Khdeir [55]. saving with respect to other discrete-layer based plate The numerical results give some insight into the elements because the element number of DOF is effect of the two degrees of freedom (the nodal values independent from the number of layers in the lamiof the curvatures of the reference surface) and the nate. numerical integration. It appears that for coarse 4. The RHQ40 plate element, and the related plate meshes (NE = 1) the accuracy of the RFQ40 is elements CQ40, FQ40, RFQ40, and HQ40 also apgreater than that of the RQ4 element: this is due to pear to work well when used to model an irregular the two additional degrees of freedom, whose benregion (distored planform). eficial effects are greater for coarse meshes where the kinematics of the single finite elements is very important in defining the kinematics of the whole plate. For REFERENCES fine meshes, it appears to be the numerical integration scheme that plays the most important role. 1. C. W. Bert, A critical evaluation of new plate theories The last results, Table 14, give some insight into the applied to laminated composites. Composite Srruct. 2, 329-347 (1984). effect of the distortion of the element planform. In 2. A. K. Noor and W. S. Burton, Assessment of shear this investigation only one quarter of plate has been deformation theories for multilayered composite plates. modelled, where NE = 2 is used. The distortion of the Appl. Me&. Rev. 41, 1-12 (1989). element has been accomplished by moving the point 3. A. K. Noor and W. S. Burton, Assessment of compuP along the diagonal in the direction of the plate tational models for multilayered anisotropic plates. Composite Struct. 14, 233-265 (1990). centre (see Fig. 4). 4. R. K. Kapania and S. Raciti, Recent advances in
analysis of laminated beams and plates. Part I: Shear effects and buckling. AIAA Jnl27, 923-934 (1989). 5. R. K. Kapania and S. Raciti, Recent advances in analysis of laminated beams and plates. Part II: Vibrations and wave propagation. AIAA Jnl27, 935-946
CL CL
I
x Ia1 = 0.375,
L-4
1t 2
+-,/2----c(
IP
x
2,
la,
=
0.375.
Scripta, Washington, DC (1975). 8. J. M. Whitney and N. J. Pagano, Shear deformation in heterogeneous
anisotropic plates. J. Appl. Mech. 37,
1031-1036 (1970). e
x1
Fig. 4. Mesh with the maximum distortion of the element planform x,, /a, = 0.375. CAS 17,1--H
(1989). 6. J. N. Reddy, A review of refined theories of laminated composite plates. Shock and Vibr. Digest 22, 3-13 (1990). 7. R. M. Jones, Mechanics of Composite Materials.
9.
V. Giavotto, Sulla meccanica di pannelli anisotropi ed eterogenei. Istituto Lombard0 (Memorie SC. Mat.) XXV, 437-480
(1969).
10. J. M. Whitney and C. T. Sun, A higher order theory for extensional motion of laminated composites. J. Sound Vibr. 30, 85-97 (1973).
104
M. DI
11. K. H. Lo, R. M. Christensen and E. M. Wu, A higher-order theory of plate deformation. Part II: Laminated plates. J. Appl. Mech. 44, 669676 (1977). 12. M. Levinson, An accurate, simple theory of the statics and dynamics of elastic plates. Mech. Res. Commun. 7, 343-350 (1980). 13. M. V. V. Murphy, An improved transverse shear deformation theory for laminated anisotropic plates. NASA TP-1903 (1981). 14. J. N. Reddy, A simple higher-order theory for laminated composite plates. J. Appl. Mech. 51, 745-752 (1984). 15. A. V. Krishna and S. Vellaichamy, On higher order shear deformation theory of laminated composite panels. Composite Struct. 8, 247-270 (1987). 16. L. Librescu and A. A. Khdeir. Analvsis of svmmetric cross-ply laminated elastic plates using a higher-order theory: Part I-stress and displacement. Composite Srrucr. 9, 189-213 (1988). 17. C. T. Sun and J. M. Whitney, On theories for the dynamic response of laminated plates. AIAA Jnl 11, 178-183 (1973). 18. S. Srinivas, A refined analysis of composite laminates. J. Sound Vibr. 30, 495-507 (1973). 19. S. T. Mau, A refined laminated plate theory. J. Appl. Mech. 40, 606-607 (1973). 20. P. Seide, An improved approximate theory for the bending of laminated plates. Mech. Today 5, 451-466 (1980). 21. M. Di Sciuva, A refined transverse shear deformation theory for multilayered anisotropic plates. Atti Accadet&a Scienze T&o 118,. 279-295 (1984). 22. M. Di Sciuva. Bendina, vibration and buckling. of simply supported thick multilayered orthotropic p&es: an evaluation of a new displacement model. J. Sound Vibr. 105, 425442 (1986). 23. M. Di Sciuva, An improved shear deformation theory for moderately thick multilayered anisotropic shells and plates. J. Appl. Mech. 54, 589-596 (1987). 24. M. Di Sciuva, Further refinement in the transverse shear deformation theory for multilayered composite plates. Atti Accademia Scienze Torino, in press. 25. M. Di Sciuva, Multilayered anisotropic plate models with continuous interlaminar stresses. Composite Struct., in press. 26. H. Murakami, Laminated composite plate theory with improved in-plane responses. J. Appl. Mech. 53, 661-666 (1986). 27. A. Toledano and H. Murakami, A higher-order laminated plate theory with improved in-plane responses. Int. J. Solids Struct. 23, 111-131 (1987). 28. J. N. Reddy, A generalization of two-dimensional theories of laminated composite plates. Commun. Appl. Numer. Meth. 3, 173-180 (1987). 29. J. N. Reddy, On refined computational models of composite laminates. Int. J. Numer. Meth. Engng 27, 361-382 (1989). 30. C. W. Pryor and R. M. Barker, A finite-element analysis including transverse shear effects for applications to laminated plates. AZAA Jnl9, 912-917 (1971). 31. P. Mantegazza and M. Borri, Elementi finiti per l’analisi di pannelli anisotropi. I’Aerotecnica Missili e Spazio 53, 181-191 (1974). 32. A. K. Noor and M. D. Mathers, Shear-flexible finiteelement models of laminated composite plates and shells. NASA TN D-8844 (1975). 33. J. N. Reddy, Free vibration of antisymmetric, angleply laminated plates including transverse shear deformation by finite element method. J. Sound Vibr. 66,565 (1979). 34. P. G. Bergan and X. Wang, Quadrilateral plate bending elements with shear deformations. Comput. Struct. 19, 25-34 (1984).
&WA
35. M. Di Sciuva, Development of an anisotropic, multilayered, shear-deformable rectangular plate element. Comput. Struct. 21, 789-796 (1985). 36. M. Di Sciuva, Evaluation of some multilayered, sheardeformable plate elements. AIAAIASMEIASCEIAHS 26th Structures, Structural Dynamics and Materials Confirence, Orlando, Florida (1985). Also Comput. Struct. 24, 845-854 (1986). 31. J. J. Engblom and 0. 0. Ochoa, Finite element formulation including interlaminar stress calculations. Comput. Struct. 23, 241-249 (1986). 38. A. T. Chen and T. Y. Yang, Static, dynamic and buckling formulation of a symmetrically laminated plate finite element for a microcomputer. J. Composite Mater. 21, 44453 (1987). 39. S. P. Wanthal and T. Y. Yang, Formulation of three simple triangular plane stress anisotropic finite elements for a microcomputer. J. Composite Mater. 22, 106&1079 (1988). 40. P. R. Heyliger and J. N. Reddy, A higher order beam finite element for bending and vibration problems. J. Sound Vibr. 126, 309-326 (1988). 41. A. S. Mawenya and J. D. Davies, Finite element bending analysis of multilayer plates. Int. J. Numer. Mel. Engng 8, 215-225 (1974). 42. S. C. Panda and R. Natarajan, Finite element analysis of laminated composite plates. Int. J. Numer. Meth. Engng 14, 69-79 (1979). 43. M. Epstein and H. P. Huttelmayer, A finite element formulation for multilayered and thick plates. Comput. Struct. 16, &45650 (1983). 44. 0.0. Ochoa, J. J. Engblom and R. Tucker, A study of the effects of kinematics and material characteristics of the fundamental frequency of composites plates. J. Sound Vibr. 101, 141-148 (1985). 45. R. L. Hinrichsen and A. N. Palazotto, Nonlinear finite element analysis of thick composite plates using cubic spline functions. AIAA Jnl 24, 18361842 (1986). 46. R. A. Chauduri and P. Seide, Triangular finite element for analysis of thick laminated plates. Int. J. Numer. Meth. Engng 24, 1203-1224 (1987). 47. D. R. J. Owen and Z. H. Li, Refined analysis of laminated plates by finite element displacement methods--I. Fundamentals of static analysis. Comma. Struct. 26, 907-914 (1987). 48. D. R. J. Owen and Z. H. Li. Refined analvsis of laminated plates by finite element displacement methods-II. Vibration and stability. Comput. Struct. 26, 915-923 (1987). 49. J. N. Reddy, E. J. Barber0 and J. L. Teply, A plate bending element based on a generalized laminate plate theory. Int. J. Numer. Mel. Engng 28, 2275-2292 (1989). 50. E. J. Barbero, 3. N. Reddy and J. L. Teply, An accurate determination of stresses in thick laminates using a generalized plate theory. Inr. J. Numer. Meth. Engng 29, l-14 (1990). 51. K. M. Rao and H.-R. Meyer-Piening, Analysis of thick laminated anisotropic composite plates by the finite element method. Composite Struct. 15, 185-213 (1990). 52. K. Washizu, Variational Methods in Elasticity and Plasticity. Pergamon Press, Oxford (1968). 53. N. J. Pagano, Exact solutions for rectangular bidirectional composites and sandwich plates. J. Composite Mater. 4, 20-34 (1970). 54. J. M. Whitney, The effect of shear deformation on the bending and buckling of anisotropic laminated plates. Composite Struct. 4, 109-121 (1987). 55. J. N. Reddy and A. A. Khdeir, Buckling and vibration of laminated composite plates using various plate theories. AIAA Jnl 27, 1808-1817 (1990).
105
A general quadrilateral anisotropic plate element APPENDIX
Interpolation
functions for u,, u,, g,, and g,
Ni=f(l+&C)(I+~i~), Interpolation
i=1,2,3,4.
functions for w K’=H,(W,(rl), NY’= 4(5)Hb1),
K=H,(I;)H,(tt) W = H, tTW&)
NYC= H&)H,(‘/), W = H&)Wf),
W = WS)H,ttl) W = H,(S)H,(rl)
W” = H,(WMrtX N;” = H&W&),
N,WI= H,(S)H&) N$ = H,WH,(I)
W~=~&)H,(~), W = H&)H;l(tl),
W’=H&W,(v) N4”((= WWMlj)
H,(r)=&(-3rS+tOr3-15r+8) Hz(r) = h(3rs-
10r3 + 15r + 8)
H,(r)=&{-3rS+r4+10r3-6r*-7r+S) Hd(r) = i(-3rs
- r4 -+ IOr + 6r2 - 7r - 5)
H,(r)=fg(-r5+P+2r3--2r2-r
+ 1)
He(r) = h(r’ f r4 - 2r3 - 2r2 + r + 1).