- Email: [email protected]
Consistent and variationally correct finite elements for higher-order laminated plate theory
Composite Structures 29 (1994) 445-456 © 1994 Elsevier Science Limited Printed in Great Britain. All rights reserved 0263-8223/94/$7.0l) ELSEVIER
Consistent and variationally correct finite elements for higher-order laminated plate theory P. Rama Mohan, B. P. Naganarayana Aircraft Design Bureau, Hindusthan Aeronautics Ltd., Structural Sciences Division, National Aeronautical Laboratory, Bangalore - 560 O17, India
&
G. Prathap Jawaharlal Nehru Centre for Advanced Scientific Research, Indian Institute of Science Campus, Bangalore - 560 012, India
Higher-order plate theories and the corresponding finite element formulations have become popular recently, owing to the increasing demand for accuracy in transverse stress predictions in thick and laminated structures. In this paper, several displacement-type finite elements based on the Lo-Christensen-Wu higher-order description of transverse structural deformation are examined. The errors of locking, delayed convergence and stress oscillations are critically studied. The accuracy of the elements is then demonstrated using some wellestablished bench-mark tests. It is also observed with reference to a highprecision three-dimensional finite element solution that the higher-order transverse deformation theories and the corresponding finite element formulations, in general, cannot capture the transverse stress distribution correctly near the structural boundary of thick beams and plates.
1 INTRODUCTION
proper a m o u n t of transverse shear energy. In addition, the transverse stress distribution can be considerably inaccurate in the case of thick and laminated structures. As an i m p r o v e m e n t over the first-order shear deformable theories, Lo et aL 3,4 have presented a plate theory based on an assumed higher-order displacement field. This theory permits a quadratic variation of the transverse shearing strains and linear variation of transverse normal strain through the thickness direction. Eater, Reddy presented a third-order shear deformable plate theory, 5 which again neglects transverse normal strain effects, but assumes a parabolic distribution of transverse shearing strains through the thickness by enforcing zero shear conditions on the lateral surfaces a priori. As a result, the formulation needs C~-continuous description for the displacement fields in the finite element formulation. Finite elements based on C J-continuous displacement fields are not desirable for use in general purpose package applications. Liao et al. ~ p r o p o s e d a hybrid stress model based on the higher-order theory to improve the transverse dis-
A n u m b e r of theories exist in the literature for the analysis of thick and layered-composite structures. T h e classical theory based on the Love-Kirchhoff hypothesis that normals to the middle surface remain normal and unstretched after deformation was the starting point in the development of general plate theories. This theory neglects the transverse shear and transverse normal strain effects in the thickness direction. Today, the inaccuracy of such a hypothesis when applied to thick and layered plates is well known. Reissner I and Mindlin 2 provided first-order shear deformable theories based on the assumption that normals to the middle surface remain straight and unstretched but need not be normal to the deformed mid-surface. Such theories allow for a constant distribution of transverse shear strain through the thickness direction but still neglect the transverse normal strain. One drawback of such theories is that a fictitious shear correction factor is required in order to maintain the 445
446
P. R a m a Mohan, B. P. Naganarayana, G. Prathap
tribution of deformations and stresses using independent fields for the displacements u, v, w and the transverse shear stresses Oxz, Oyz. Again, the stress-based finite element formulations are not convenient for use from the point of view of general purpose finite element packages. The higher-order displacement field proposed by Lo et al. 3,4 requires only C°-continuity of the solution. This permits a great deal of flexibility in the specification of approximation functions over an element and the corresponding finite elements are best suited for general purpose applications. Unfortunately, most of these elements become too stiff when used to model thin structures and some elements become completely rigid as their span-to-thickness ratio becomes very large. This type of behavior is called artificial stiffening or locking.7,8 Much effort has been made to understand clearly and to eliminate the source of the errors of locking and the associated stress oscillations 7,~ especially in the case of finite elements developed using first-order shear deformable theories.9, ,~ Kant et al. ~ have presented a C°-continuous finite element formulation for the higher-order theory proposed by Lo et al. 3"4 A uniform reduced integration strategy is adopted for modeling the shear stiffnesss of the element, as an attempt to eliminate locking in a heuristic way. This paper examines critically the problems of shear locking and stress oscillations in the case of finite elements developed using the general higher-order displacement field proposed by Lo et al. 3'4 The sources of these errors are explicitly located and eliminated by applying the concepts of field consistency7"8 and variational correctness,2,13 in the case of two-noded and threenoded beam elements (B2 and B3) and the observations are then extended to the case of the four-noded and nine-noded plate elements (Q4 and Q9). The distribution of interlaminar shear and normal stresses across the thickness are evaluated at the centroid of the quadratic elements (B3 and Q9) by using equilibrium equations while the constitutive equations are used for recovering the transverse stress distributions in the case of the linear elements (B2 and Q4). Numerical differentiation is used for evaluating the second derivatives in the former case. Thus the one- and two-dimensional quadratic elements are able to give quartic variations of shear stress and cubic variation of the bending stress across the structural thickness.
The transverse stress distribution captured by such finite elements is then compared with a corresponding high-precision three-dimensional finite element solution. It is observed that the present finite elements are very accurate throughout the structural domain except near the boundary, particularly when it is clamped. In ,fact, none of the higher-order transverse deforrnation theories and the corresponding finite elements can recover accurate transverse stress distribution across the thickness near the structural boundary. This is because the transverse stress distribution needs a much higher order of description when compared to the quartic shear stress and cubic bending stress distribution (that can be generally achieved by the higher-order transverse deformation theories and the finite elements~ near the structural boundary. In addition, most of the higher-order theories assume zero shear stress at the free (top and bottom) surfaces of the beam/ plate everywhere. This is not true in the case of the clamped boundary where, in fact. the transverse stress distribution becomes singular near the top and bottom surfaces, which act as re-entrant corners. Thus. enforcement of quartic shear and cubic bending stress distribution and zero shear stress on the free surfaces lead to extraneous oscillations in the transverse stress distribution near the structural boundary. As a result, the higher-order transverse deformation theories and the corresponding finite elements, in general. become unreliable near the structural boundary.
2 FINITE ELEMENT FORMULATION In this paper, the higher order CU-displacement field description proposed by Lo et al. 3"4 is chosen for the finite element formulations. A family of layered-composite finite elements -- two-noded linear and three-noded quadratic beam elements CB2 and B3) and four-noded linear and ninenoded quadratic rectangular plate elements (Q4 and Q9) - is developed for a close examination. 2.1 Field description The displacement field proposed by Lo e t al. 3"4 consists of in-plane displacements expanded as cubic functions in the thickness coordinate. 3 The polynomial expansion for transverse displacement is truncated at one order lower the the expansion for in-plane displacements. It can be seen (see eqns (2e) and (2f) below) that this choice leads to a
Higher-order laminated plate theory consistent definition of interpolation for the transverse shear strains with reference to the thickness coordinate z. Here, the same assumptions are used in the finite element formulations such that the three-dimensional displacement field is defined using the mid-surface degrees of freedom u, 0,..., w* as follows:
The stress-strain relations for the L th lamina in the laminate coordinate axes system can be obtained by the usual transformation rules ~4 and are given by
OxxlL11012 [ 3Q16L[xqL
U(x, y, z)= u(x, y) + zO(x, y) + z2u*(x, y) + z 30*(x,y)
V( x, y, z) = v( x, y) + z qk(x, y) + z2v*( x, y) + z3 q~*(x, y)
(2a)
~ ~ , + z..3.h 6yy = V,v=V,v-Jr-Z~,v-Jr-Z2V,y ~,y,
(2b)
gxy= g,y-+- V,x~-(bl, yq-V,x) q- z( O,y-}- ~,x)
+z2(U,y+V,*)+z3(O,*+¢, *)
(2c)
6~ = W,: = ~p+ 2zw*
(2f)
Each layer in the laminate is assumed to be in a three-dimensional stress state so that the constitutive relationships for a typical lamina L with reference to the fiber-matrix coordinate axis (1-2-3) can be written as 14
111L { [clci2c30IL C12 C22 C23 0
L Q45 Q55 J [ 6xzj
dS.
(3)
The elemental equations can be obtained by minimizing the potential energy function FI with respect to the unknown nodal displacement vector. 2.2
T h e linear t w o - n o d e d b e a m e l e m e n t (B2)
The linear two-noded beam element with seven degrees of freedom per node is shown in Fig. l(a). It is assumed that V=0, U - - U ( x , z ) a n d W= W (x,z) such that the strain field given in eqn (2) reduces to:
.~622~
/
6xx
0
Q16 Q26 Q36 Q66
II = -2 va~jeodV-
1]),y)
+ z2(3~ * + w,,*).
0
OxyJ
(2e)
V,~+ W,y=(~b+ W,y)+ Z(2V* +
a12J
Q13 Q23 Q33 Q36
(2d)
+ z2(30 * + w,*)
022~ :
oz~/
in which {a,x, ayy, o~z, axy, oyz, o~z}and {6xx, 6yy, 6~z, e,y, eyz, 6~z} are the stress and strain components, respectively, with respect to laminate axes and Qij are the three-dimensional elastic constants in the laminate axes of the L th lamina. The potential energy function I1 of a laminated composite finite element of surface area S and volume V subjected to a transverse loading p(x, y) on a surface which is at a distance z = ¢ from the mid-surface can be expressed as follows
exz = U,~+ W,~=(0+ w,~)+z(2u*+ q),x)
6 yz
Q,z Q22 Q23 Q26
0-xz
W(x, y, z) = w(x, y) + Z~p(x, y) + Z 2 w*(x, y).
6x,:= U,, = u,~+zO,,+z:U*,x+Z30, *
]
Oyy} =
{'}L:rQ 4 4 ILI yzl L
(1)
The in-plane and transverse strain components associated with the displacement field in eqn (1) are given by:
447
0
C66
1612J
U,x=U,x+ZO, x+Z2U,*+z3¢, *
6 zz -~- W,z = ~ + 2zw*
(4a) (4b)
6 xz -~- U,~+ W,x=(0+ w,,)+z(2u*+ qJ, x)
0131
C55 [6131
where {01j , 022, 0-33, 0-12, 0-23, 0"13} are the stress and {ell, 622, 633, 612, 623, 613} a r e the strain components referred to the lamina coordinate axes and C 0 are the elastic constants of the L th lamina.
+z2(30"+
(4c)
Eyy= e,y = eyz = 0.
(4d)
The displacement components for the element are interpolated using linear Lagrangian interpolation functions so that, in an isoparametric formulation, we can express the active degrees of
P. Rama Mohan, B. P. Naganarayana, G. Prathap
448
I-
(doF)
I.
.I
~
= (u,O, ta",e',w,P,w'}
(a) 2 - n o d e d
2~--~ -
~,~,w,P,w" }
(b) 3-noded beam etement
bear~ eteMen±
3
4
=
g .I
4
7
3
9~
e
0
0
%
0
=
qJ
2
(do~') =
v,¢~ v~,~ ~,
w,p~W ~")
Wjp, W"~")
(c) 4-noded
~2
pl~te etenen±
(d) 9-noded pro.re
ete~en-t
Fig. 1. The family of higher-order elements: (a) two-noded beam element; (b) three-noded beam element; (c) four-noded beam plate element; (d) nine-noded plate element.
freedom as:
{} a°21 U*
/a21 a22
0*
laB, a32
{:} I::::::]{;l w*
(5)
L b2,622J
where the constants a0~,..., b22 are functions of the nodal degrees of freedom. Substituting eqns (5) in eqn (4c) we can express the transverse shear strain as
gxz={{all +b02/l}5-a12~} + z{ {2a2, + bt2/l} + 2a22~}
(6)
+ z2{{3a3, + b22/1} + 3a32 }. If this is written in terms of Legendre polynomials, ~°,~3 it can be shown that the transverse shear strain energy comprises squares of the coefficients of this expansion. It is then simple to identify the constraints on the transverse shear strain in the thin plate limits as
True all + bo2/l~O
Spurious a12--,0
2a21 + bt2/l--,O
a22~0
3a~ 1+ b22/l-,0
a32-, 0.
i7)
From eqns (7), we can argue that the true constraints give the correct Kirchhoff's condition at the element centroid, i.e. (ex=)~=0 ~ 0. On the other hand the spurious constraints reflect a (spurious) condition equivalent to Oe,,x/Oz~O. Thus the spurious constraints lock the element from bending and lead to highly underestimated displacements and bending moments. It can be shown that the locking error is directly proportional to (1/t) 2 as in the case of first-order elements. 7 Again from eqn (6), it can be observed that these spurious constraints lead to large linear oscillations in the shear stress predictions in the case of the inconsistent element which we identify as B2.0. The errors of locking and stress oscillations can be analytically derived in an a priori sense as was done for the corresponding first-order two-noded beam element, v,8 However, owing to the identical nature of the errors, the relevant derivations are not presented here. The preceding observations are confirmed later by a proper numerical example. From the method of Legendre polynomial expansion, ~°,~-~ the variationally correct way of eliminating the inconsistent terms, having associated with the spurious constraint, from the transverse shear strain field is to simply drop the inconsistent linear Legendre terms in ~ to eqn (6) such that the consistent shear strain field becomes,
r~={alt + boz/l} + z{2a2~ + b12/l} +z2{3a31 + b22/1}. It is demonstrated later, using appropriate numerical examples, that such a consistent element, B2.1, is free of locking and stress oscillations. B2.1 exhibits the variationally correct convergence rate and constant stresses over the element length.
2.3 The quadratic three-noded beam element (B3) The quadratic three-noded beam element with seven degrees of freedom per node is shown in Fig. l(b). The same assumptions that were employed in the case of B2 a r e used here ~also
Higher-order laminated plate theory such that the strain field given in eqns (4) remain valid. The displacement components for the element are interpolated using the quadratic Langrangian interpolation functions so that, in an isoparametric formulation, we can express the active degrees of freedom as:
{} Ia°la°a°]l =
u* O*
al 1 a12 a13 a21 a22 a23 a31 a32 a33
w*
~2
[_b2, b22 b23
(8)
ex== {a~ l + at3~3 + bo2/l} + {at2 + 2bo3/l}
+ Z{ {2a21 + 2a23/3 + bl2/l} + {2a22
+2b,3/l}~} + z2{ {3a3t + a33 + b22/l } + {3a32
+ 2b23/l}~}.
× ~-(a13/3)(1 - 3~ 2)
It is demonstrated later, using appropriate numerical examples, that such a consistent element, B3.1, is free of delayed convergence and stress oscillations.
+ z{ {2a21 + 2a23/3 + b12/l} + {2a22 + 2b13/l} × ~-(2a23/3)(1 - 3~ 2)}
× ~-- a33 ( 1 - 392)}.
inconsistent quadratic element (B3.0) is equivalent to that of the consistent linear element (B2.1). From eqn (9), it is clear that the spurious constraint is associated with the quadratic Legendre term and hence the inconsistent element suffers large quadratic oscillations while predicting the transverse shear stress. These observations are established using an appropriate numerical example later. From the method of Legendre polynomial expansion, the variationally correct way of eliminating the inconsistent terms from the transverse shear strain field is to simply drop the inconsistent quadratic Legendre terms in ~ from eqn (9) such that the consistent shear strain field becomes rxz = {a~, + a13/3 + bo2/l} + {a,e + 2bo3/I}~
~2
Substituting eqns (8) into eqn (4c) we can express the transverse shear strain as
-+- Z2{ {3a3, + a~3 + b22/l } + {3a32 +
449
2b23/1 } (9)
Again, from the method of Legendre polynomial expansion the constraint on the transverse shear strain energy leads to three discretized constraints in the limits of vanishing thickness. It can be observed that the constraints associated with the constant and linear Legendre terms are the true discretized constraints representing the Kirchhofff condition in the limits of vanishing thickness. On the other hand, the third constraint associated with the quadratic Legendre term is a spurious constraint on the linear part of the bending strain, i.e.
Thus, such an inconsistent element, B3.0, can capture only the constant bending moment part over the element length and the linear bending moment part will be highly underestimated. In other words, the inconsistent quadratic element, B3.0, does not lock as its linear counterpart B2.0 but has a lower convergence rate than the expected variationally correct rate. It can be observed that the rate of convergence of the
2.4 The linear four-noded plate element (Q4)
The four-noded plate element configuration is as shown in Fig. 1(c). The element displacement field is described by 11 degrees of freedom: u, 0, u*, 0", e, ~b, v*, ~b*, w, V2, w*. The regular bilinear Lagrangian shape functions are used for displacement field interpolation. It can be shown that the shear strain fields exz and ey= (eqns (2e) and (2f)) can be expressed as
Sx==Ao + Al ~ + A2rl + A3~rl e.,,~= B o + B~ ~ + B2rl+ B.~r I. For an undistorted (i.e. rectangular) element configuration, the terms associated with ~, s~r/in the shear strain field ex., (i.e. constants A~ and A3) are field-inconsistent and lead to spurious constraints in the penalty limits. Similarly, the constants B2 and B3 in the shear strain field % are field-inconsistent and lead to spurious constraints in the penalty limits. These spurious constraints lead to the conditions
ae,.,laz--, O,
Oe,>./Oz~ 0,
Oev/az-~ O.
Thus the inconsistent element, Q4.0, suffers from locking, highly underestimated bending moments
450
P. Rama Mohan, B. P. Naganarayana, G. Prathap
and large oscillations in shear stress predictions. Thus the shear strain fields have to be smoothed to their respective consistent form as Lz = A0 + fi'2 q and ryz =/~0 + B1 ~" Any variationally correct method of field-redistribution 1° can be used for the purpose. The fieldredistribution is equivalent to the use of 1 × 2 and 2 × 1 Gaussian quadrature rules for the shear strain energy components associated with the shear strain fields S~z and eyz, respectively, for a rectangular element. The field-consistent element, Q4.1, is free of the errors of locking and stress oscillations.
2.5 The quadratic nine-noded plate element (Q9) The element configuration for the nine-noded plate element is shown in Fig. l(d). Again the element displacement field description is based on 11 degrees of freedom. It can be shown that the terms associated with the quadratic Legendre terms ( 1 - 3~ 2) are field-inconsistent in exz while those associated with ( 1 - 3 r / 2 ) are field-inconsistent in eyz. The corresponding spurious constraints lead to the conditions of the type
Thus the linear variation in the bending moments is highly underestimated and as a result, the field-inconsistent element, Q9.0, exhibits delayed convergence, constant bending stresses and large quadratic oscillations in shear stress predictions similar to its first-order counterpart. ~3 The shear strain fields are expanded as Legendre polynomials on certain principal lines, and the inconsistent quadratic Legendre terms are dropped to get the consistent element version, Q9.1. The consistent element is free of delayed convergence and gives correct linear variation to the bending and shear stress fields as shown later in numerical examples.
2.6 Stress recovery The transverse stress field recovered from the constitutive relationships in a displacement-type finite element formulation is discontinuous at the interlaminar region of laminated structures. The equilibrium equations are often integrated through the thickness of the structure to give the
desirable continuous transverse stress distribution across the thickness. Such a procedure needs first derivatives of the normal stresses with reference to the in-plane coordinates; this cannot be done directly in the case of linear finite elements. In this paper, the stresses are calculated using the constitutive relationship for the linear elements (B2 and Q4) while equilibrium equations are used for predicting the shear stresses from the normal stresses that are recovered from the corresponding constitutive relationships in the case of the quadratic elements (B3 and Q9). All the stress components are recovered at the element centroids in the case of linear elements while the inplane stresses are sampled at the integration points corresponding to reduced Gaussian quadrature and the transverse stresses are sampled at the element centroids in case of quadratic elements. In the next section, the quadratic elements are used for solving some well-established benchmark tests, since they naturally can produce more accurate stress distribution across the structural thickness when compared to linear elements. These elements can sense quadratic, cubic and quartic variations for the transverse normal~ inplane normal and transverse shear stress components respectively. It is observed that the order of stress distribution mentioned above becomes insufficient near the structural boundary. In addition, the equilibrium equations are generally solved with an additional assumption that the shear stresses vanish on the top and bottom surfaces. This assumption is invalid particularly near the clamped boundary where, in fact. the stress distribution becomes singular near the top and bottom corners, which act as re-entrant corners. in the case of thick structures. Thus imposition of quartic distribution and vanishing transverse shear stress conditions on either the top and bottom surfaces leads to extraneous perturbations in the transverse shear distribution near a boundary as will be demonstrated later. But this behavior is inherent in general higher-order transverse deformation formulation and consequently they cannot recognize the stress singularity near the reentrant corners on the clamped boundary. However, the present element formulations are found to be very accurate away from the structural boundary when compared to a high-precision three-dimensional finite element solution. For comparison of the transverse shear stress distribution obtained from the C ° elements pre-
Higher-order laminated plate theory sented in this paper, a C 1-continuous theory proposed by Krishnamurthy ~5 is considered. This theory can recover constant shear and cubic bending stress distributions across the thickness near the clamped boundary; and quadratic shear and linear bending stress distributions across the thickness away from the clamped edges. It is observed later, that the transverse stress distribution is not accurate anywhere when compared to a corresponding high-precision three-dimensional finite element solution. ~6 However, the transverse shear stress distribution at the clamped edge in the interior region is most reliable when compared to other higher-order transverse deformation theories.
451
X
X
Theory/B3.1/Q9.1 B2.0/Q4.0 B3.O/Q9.0/ E32.1/Q4.1
*
+ x
co
'
'
'
I
'
3
'
'
I
'
'
'
I
'
'
5
4
II
'
I
'
6
'
'
t
'
7
'
I
8
No.EL
(a) * 1 element case + 2 element case
//*
3 NUMERICAL EXAMPLES In this section, we consider first an example to demonstrate the errors of locking and stress oscillation in field-inconsistent elements. Then we consider a problem to establish reliability of higher-order transverse deformation theories and finite element formulations in recovering accurate transverse stress distribution near structural boundary. Finally we consider a few well-established laminated plate problems to demonstrate the efficiency of the new field-consistent ninenoded plate element, Q9.1.
'
/ ~.~-
g--.@
e4-
3
4
5
o90/t5 (b)
6
7
8
Fig. 2. (a) Convergence rate with different consistent and inconsistent elements. (b) Locking error with reference to thickness ratio (l/h).
3.1 Cantilever slender beam under tip shear force A cantilever beam of length l = 10, breadth b = 1 and thickness t = 1-0 subjected to a tip shear force P = 1, is modeled with different consistent and inconsistent elements. Young's modulus of the material is 10 while the Poisson's ratio is zero. The rate of convergence with different elements is depicted in Fig. 2(a). The inconsistent linear elements (B2.0, Q4.0) lock severely as expected. In fact, the locking error is proportional to the square of the aspect ratio of the elements as in the case of their first-order counterparts (see Fig. 2(b)). The consistent linear elements show the variationally correct convergence rate (the discretization error for the present problem modeled with linear elements is given by 1/4n 2 where n is the number of elements used). It can also be noticed that the inconsistent quadratic elements (B3.0, Q9.0) do not lock but show a delayed rate of convergence that is equivalent to that of con-
sistent linear elements. The consistent quadratic elements (B3.1, Q9.1) give an exact solution even with a single-element model. The bending moment predictions from various element versions are plotted in Fig. 3. Again, as in Section 2.2, the inconsistent linear elements lock severely showing highly underestimated bending moments. The consistent linear elements show correct bending moments at the element centroids. On the other hand, the inconsistent quadratic elements cannot depict the linear variation of the bending moment in the penalty limits and as a result the bending moment is correct only at the element centroids. The consistent quadratic elements produce variationally correct linear bending moment variation over the element domain. Figure 4 shows shear force predictions from various elements. It can be observed that the
452
P. Rama Mohan, B. P. Naganarayana, G. Prathap j
£~
*
~o
+ X o
\
2
Theory/B3.1/Q9.1 B3.O/Q9.0 B2.1/Q4.1 B2.0/Q4.0
/ /
\\
/ "\
Fig. 5. \\ \
2
4
6
8
10
Bending moment distribution with different consistent and inconsistent elements.
* + x
Theory/B2.1/B3.1/ Q4.1/Q9,1 B2.0/Q4.0 B3.0/Q9.0
7/J
/
~< ,/
/
//
I 0
2
4
6
8
10
x
Fig. 4.
x
A deep cantilever beam under uniformly distributed shear body force.
\
x
~o. , ,* ' \
f
I ! I
i \
Ot--i
intensity in z-direction
:J \, \
~
body force of unit
/
\
Fig. 3.
Uniformly distributed
// \
0
"
J/"
Shear force distribution with different consistent and inconsistent elements.
inconsistent linear elements show large linear oscillations while the inconsistent quadratic elements show large quadratic oscillations while all the consistent element versions are able to predict the constant shear stress condition over the structural domain correctly. It should be noted that the behavior of different element versions is very similar to that of their first-order counterparts.
3.2 Cantilever deep beam under uniformly distributed shear load A deep cantilever beam of length 1= 10, breadth b= 1 and thickness t = 10 (Fig. 5) is subjected to uniformly distributed body force of unit intensity in the z-direction. The material properties are the same as in the previous example. The structure is
modeled by using consistent nine-noded plate elements (Qg.I) and the eight-noded consistent hexahedral elements (H8) of Ref. 16. A n adaptively graded mesh generated by the three-dimensional finite element package ~6 is used as the basis for both meshes. The finite element mesh is highly refined near the top and bottom corners at the root of the cantilever beam by using mesh gradings in the ratio of 1 : 2 : 4 : 8 : 1 6 : 3 2 : 6 4 : 1 2 8 and 1:1:4:8:8:4:1:1:1:1:4:8:8:4:1:t in the xand z-directions, respectively, in the case of H8. The same mesh grading is retained in the xdirection while 20 layers of equal thickness are used in the z-direction in the case of Q9.1. Stresses are sampled at the centroids of H8 elements while stress recovery is done at the interlaminar points on the centroidal locations of Q9.1 elements. The transverse shear and bending stress distributions across the structural thickness are sampled at the locations x=0-0196, x = 0 . 9 3 6 , x = 1.843 for both the models. The transverse shear stress distributions at different locations from Q9 and H8 models are plotted in Fig. 6 at locations x = 0"0196. x = 0" 196 and x = 1.843. respectively. The stress distributions derived from the classical beam theory and from another independent higher-order shear deformation theory ~5 are also plotted for comparison. As discussed in Section 2.5, Q9.1 shows very accurate quartic shear stress variation when compared to the three-dimensional solution (HS/ away from the clamped edge. However. near the clamped edge Q9.1 shows extraneous quartic oscillations such that the shear stress is not correct at any predetermined point across the thickness. On the other hand, many higher-order theories and the corresponding finite element formulations show quadratic shear stress variations (which are inaccurate when compared to the three-dimensional solution at all locations) across the thickness throughout the beam length. Krishnamurthy's theory provides a transverse shear
Higher-order laminated plate theory ---.... ---
Q9.1 HEXA8 EBT AVK
\,,\ to
N
]
4
/
JJtl _..P
o , , ~,'?-7-,--T--,__,-:~,-q 0
5
10
~z
15
20
Ca) --
Q9.1 HEXA8
-
N
t
I
_b. o
,
,
,
i
,
,
,~-~
~ ,
,
4
2
,
i
6
8
,
,
,
i 10
,
;
,
i 12
,
,
; 1 14
4 53
stress distribution which is constant at the clamped edges and gradually becomes quadratic away from them, as shown in Fig. 6. The accuracy of the shear stress distribution obtained from Krishnamurthy's theory is very good in the interior region at the clamped edge when compared to the three-dimensional solution. However, the stress distribution does not compare well with the three-dimensional solution away from the clamped edge as shown in Fig. 6(c). Again, this cannot recognize the stress singularities at the reentrant corners on the clamped edge. The shear stress distribution at z = 9-91 along the axis of the beam is obtained using the threedimensional finite element mesh and is plotted in Fig. 7. The three-dimensional finite element solution clearly shows the singular nature of the stress near the clamped edge. However, as a rule most of the higher-order finite element formulations cannot depict the stress singularity near the clamped edge, as shown in Fig. 6. All the higher-order solutions can recover the cubic variation of the bending stress across the beam thickness accurately. However, the bending stress becomes singular near the top and bottom surfaces at the clamped edge. Again, the threedimensional solution can capture this nature correctly while all the beam/plate theories fail to do the same. Krishnamurthy's theory is unreliable even at locations away from the clamped edge giving a linear variation to the bending stress, which is actually cubic in nature, over the thickness. However, owing to the similar behavior
CF-X Z
(b) o
o
¢o
----
_
EBT (O-xx) HEXA8( 0;',x) HEXA8 ( d~xz)
o N
----
A8
o
b...... o
, , 1 1 1 , , [ , , , i , i 2
4
, j , , , ] , l , ] , , i 6
8
10
12
\,
j 14
%z (c) Fig. 6. Transverse shear stress distribution for a deep cantilever beam: (a) x = 0"0196; (b) x = 0.1960; (c) x = 1.8432.
5
10
15
20
x Fig. 7.
Bending and shear stress distribution along the axial direction at z = 9"91.
P. Rama Mohan, B. P. Naganarayana, G. Prathap
454
when compared to transverse shear stress distribution, the results are not shown here. The transverse stress distributions recovered from the higher-order transverse deformation elements are perturbed near the free edge also. Since the magnitude of the stresses is infinitesimally small near a free edge, the impact of the inaccuracy is negligible. Again, the results are not shown here.
solutions 17 and Reddy's solutions? From the results, one can notice that the present solutions are in close agreement with Pagano's solutions compared to Reddy's solution, especially the transverse shear stress. This may be attributed to the fact that imposition of any additional constraints in a displacement approach, as in Ref. 5. reduces the flexibility of the model to a certain extent.
3.3 Simply supported 0/90/0 rectangular plate under sinusoidal load
3.4 Simply supported 0/90/0 square plate under sinusoidal load
A simply supported, three-ply (0/90/0) rectangular laminate with length-to-breadth ratio a/b--3 and layers of equal thickness, subjected to sinusoidally distributed transverse load is considered for comparison of maximum deflection and interlaminar stresses. The material properties used are given by
A simply supported, square 0/90/0 laminated plate with a length-to-thickness ratio (a/h) equal to 4 subjected to sinusoidal distributed transverse load on the top surface is considered for validation of interlaminar stresses. Each layer of the laminate is of equal thickness and the material constants with respect to its principal material axes are specified as: E ~ = 2 5 × 1 0 6 psi, Ee2=E33=106 psi, GI2=GI3=0"5×lO ~ psi, G23 = 0"2 × 106 psi, vl2 = v23 = v~ = 0"25. For comparison the following normalized quantities are used in the study: {6x~, d~ = (1/pos) {oxz, aye}; dz=(1/po s2) Oz where s = a/h and P0 is the maximum load intensity at the plate centroid. Figure 8 shows the distribution of 6x:, 6~: and 6. across the thickness and is compared with the solution of Pagano,17 and Liao et al. 6 From Fig. 8 it can be observed that the present solutions are in close agreement with Pagano's elasticity solution, especially the distributions of 6xz and 6~. Here, 6: is sampled at the plate centroid (a/2,a/2) while 6~= and 6y; are sampled at the mid-points of the plate edges y= a/2 and x = a/2, respectively.
E I / E 2 = 2 5 , GI2/E 2 = 0"5, G23/E2 = 0.2,
E2 = E3, G13 = Glz, '1,'12= 1,'23= "V13-----0"25.
The non-dimensionalized deflections vP and interlaminar shear stresses dxz, 6yz, are given by: v~=100 (E2h3/qoa4)w; {6~,d~zt=(h/qo a) {axe, aye} where a is the edge length of the square plate and h is its thickness, and are presented in Table 1. Here, • is sampled at the plate centroid while 6~, and 6y~ are sampled at the midpoints of the plate edges y= b/2 and x = a/2, respectively, at z=0. Results obtained for different thickness ratios (a/h) are compared with Pagano's elasticity
3.5 Simply supported 0/90 square plate under sinusoidal load Table 1. Non-dimensionalized deflections and stresses in a rectangular laminate under sinusoidal load (0/90/0, h~= h/ 3, b/a = 3)
a/h
Source
~
b~=
6v~
4
Q9.I Ref. 17 Ref. 5 Q9.1 Ref. 17 Ref. 5 Q9.1 Ref. 17 Ref. 5 Q9.1 Ref. 17 Ref. 5
2'626 2.820 2-641 0.867 0.919 0.862 0.595 0.610 0.594 0.507 0.508 0.507
0-3718 0.3870 0"2724 0"4220 0.4200 0.2859 0.4314 0-4340 0.2880 0.4344 0.4390 0,2885
0"0317 0"0334 0.0348 0.0146 0.0152 0.0170 0-0116 0.0119 0.0139 0.0106 0.0108 0.0129
10 20 100
Finally, a simply supported square cross-ply, antisymmetric (0/90) laminated plate under sinusoidal transverse load is considered for comparison of deflections and transverse stresses. The set of material properties used is the same as that specified in Section 3.3 The non-dimensionalized deflections w and interlaminar shear stresses 6x., 6yz, with ~ = 10 (E2h3/qo a4) wand {oxz, 6yZ}=( h/qo a) {oxz, oy~}are presented in Table 2. Here, ~ is sampled at the plate centroid while ox~ and oy~ are sampled at the mid-points of the plate edges y= b/2, z = 0-2 h and x= a/2, z = - 0"2 h, respectively. The results for maximum stresses and deflection are corn-
Higher-order laminatedplate theory -
-
.
.
.
.
.
.
-
-
Table 2. Maximum deflection and transverse stresses for a simply supported unsyrnmetric cross-ply (0/90) square plate under sinusoidal load
Q9.1 Pagano Liao HOPE
-
.
e4
¢5
5
a/h
Source
v~
G~
6yz
4
Q9.1 Ref. 17 Ref. 18 Q9.1 Ref. 17 Ref. 18
0"2152
0'3091 0'3127 0.2868 0'3318 0.3310 0.2957
0"3028 0.3188 0.2865 0'3315 -0.2956
10
LI c
¢-1
455
-
-
0'2020 0"1264 -0.1220
/5- - / I 0.2
0.0
0,4
--
~o
0.6
d
(a)
-
(25-
~
-
Q9.1 T.Kant Pagano
-
~.
c-4
o-
Q9.1 PagGno LiGo
--
~/
.....
I ....
/
HOPE
r ' / ~ l ' ' ' ' r
....
Gz
I ' ' Frq
.c N
o d
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
Fig. 9. Simply supported 0/90 square plate under sinusoidal load. Transverse distribution of normalized shear stress dx~.
0.05
0.00
0.10
0.15
0.20
0.25
(b)
Q9.1 Pagano Liao/HOPE
..... ,*
/
)j~//-
/
,:::5
o c~
~,'1
....
i'''t
....
I ....
i ....
I ....
i ....
I
c5 I
¢5 ..// 0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
pared with Pagano's three-dimensional elasticity solutions and are presented for different a/h ratios. T h e distribution of 6, z, through the thickness is shown in Fig. 9 and are c o m p a r e d with three-dimensional elasticity ~7 and Kant's ~8 solutions. T h e present solutions are close to the threedimensional elasticity solutions of Kant. j8 This m a y be attributed to fast convergence of results because of consistent shear strain field definitions which is equivalent to the use of 2 x 3 and 3 x 2 Gaussian rules for the shear energy c o m p o n e n t s associated with the shear strain fields exz and eyz, respectively. It should be noted that a uniformly r e d u c e d integration rule (2 x 2) is applied for the whole shear strain energy in Ref. 18 which, in turn, results in spurious zero energy m o d e s in the case of the n i n e - n o d e d plate element.~3
(c) Fig. 8. Simply supported 0/90/0 square plate under sinusoidal load: (a) transverse distribution of normalized shear stress 6~:; (b) transverse distribution of normalized shear stress 6~; (c) transverse distribution of normalized normal stress 6 z.
4 CONCLUSIONS In this paper, the errors related to the finite elem e n t formulations based on the C°-continuous
456
P. Rama Mohan, B. P. Naganarayana, G. Prathap
higher-order description of the displacement 3 are critically examined. The cause for the errors of locking, delayed convergence and stress oscillations present in a conventional finite element formulation corresponding to such a higher-order transverse deformation theory has been explicitly identified and eliminated for linear and quadratic beam elements. The consistency and correctness rules are then extended to formulate linear and quadratic laminated plate elements which are free from the errors of locking and stress oscillations. It is demonstrated through numerical examples that the errors for such higher-order finite element formulations follow very similar patterns as in the case of the corresponding first-order finite element formulations. It is also observed that the stress distribution across the structural thickness from any higher-order theory is not reliable near the structural boundary. Particularly these theories and the corresponding finite element formulations cannot depict the singularities in transverse stress distribution that can arise due to the boundary conditions. The present consistent higherorder transverse deformation finite elements are shown to be highly accurate in predicting transverse stress distribution in some well-established laminated plate problems. So far, discussions have been confined to rectangular plate elements. It is highly desirable to have higher-order plate elements which can assume an arbitrary quadrilateral planform. Work is now under way to apply the edge-consistency corrections 19 (which were very effective for the first-order plate elements) to the higher-order plate elements as well.
ACKNOWLEDGEMENTS The authors are grateful to Prof. R. Narasimha, Director and Dr B. R. Somashekar, Head, Structural Sciences Division of the National Aeronautical Laboratory, Bangalore, India for their constant encouragement and keen interest in the subject. Thanks are also due to Dr Kota Harinarayan, Director, Aeronautical Development Agency for permission to use the computer facili' ties at his organization.
REFERENCES 1. Reissner. E., The effect of transverse shear deformation on the bending of elastic plates. ASME J. Appl. Mech. 12 1945 A69-A77. 2. Mindlin. R. D.. Influence of rotatory inertia and shear deformation on flexural motions of lsotropic elastic plates. ASME J. Appl. Mech.. 18 ( 1951 31-8. 3. Lo. K. H., Christensen. R. M. & Wu. E. M.. A high-order theory of plate deformation, part i: homogeneous plates. ASMEJ. Appl. Mech.. 44 19771663-8. 4. Lo. K. H.. Christensen. R. M. & Wu, E. M.. A high-order theory of plate deformation, part 2: laminated plates. ASME J. Appl. Mech.. 44 t 1977) 669-76. 5. Red@. J. N.. A simple higher order theory for laminated composite plates. J. Appl. Mech.. 51 ( 1984 745-52. 6. IAao. M. L.. Jing, H. S. & Hwang, M., Improvements on the higher order plate element with partial hybrid stress model. Cornput. Struct.. 42 l) ( 1992 45-51. 7. Prathap, G. & Bhashyam. G R.. Reduced imegration and the shear-flexible beam elements. Int. J. Nurner. Meths. Engng, 18 (1982) 195-210. 8. Prathap, G.. Field-consistency -- toward a science of constrained multi-strain-field finite elemenl formulauons. S A D H A N A . Proc Engng Sci. lndian Acad. Sci.. 9 1986) 345-53. 9. Ramesh Babu. C.. Field-consistency in the finite element formulation of multi-strain-field problems in structural mechanics. PhD thesis. IIT. Madras, 1986. 10. Naganarayana, B. R. Consistency and correctness principles in quadratic displacement type finite elements. PhD thesis. Indian Institute of Science. Bangalore. 199 I. 11. Kant. T.. Owen. D. R. J. & Zienkiewiez. O. C.. A refined higher-order C t' plate bending element. Comlmt. Struct., 15,11982) 177-83. 12. Prathap. G.. A variational basis for least-squares fieldredistribution of strain functions in the finite element formulation of constrained media elasticity. TM ST 8801. N A L Bangaiore. India. 1988. 13. Naganarayana. B. P.. Prathap, G.. Dattaguru. B. & Ramamurthy, T. S.. A field-consistent and variationally correct representation of transverse shear stratus in the 9-noded plate element. (~mgmt. Meths. Appl. Mech. Engng, 97 374 (1992) 355-74. 14. Jones, R. M., Mechanics of Composite Materials. McGraw-Hill, 1975. 15. Krishnamurthy, A. V., Toward a consistent beam theory. AIAA J., 22 (6)(1984) 811-16. 16. Prathap, G. & Naganarayana, [3. P., 3D-FEES -- lheoretical Manual, PD ST 9005, N A L Bangalore, India, 1990. 17. Pagano. N. J., Exact solutions /or rectangular bidirectional composite and sandwich plates. J. Cornp, Mater., 4 (1970) 20-34. 18. Kant, T. & Pandya, B. N., Finite clement stress analysis of unsymmetrically laminated composite plates based on a refined higher-order theory. Proc. Int. Conf. Composite Materials and Structures, 6 - 9 January 1988, Indian Institute of Technology, Madras, India. 19. Prathap, G. & Somashekar, B. R., Field- and edgeconsistency synthesis of a 4-noded quadrilateral plate bending element. Int. J. Numer. Meths. Engng, 26 (1986) 1693-708.
Consistent and variationally correct finite elements for higher-order laminated plate theory P. Rama Mohan, B. P. Naganarayana Aircraft Design Bureau, Hindusthan Aeronautics Ltd., Structural Sciences Division, National Aeronautical Laboratory, Bangalore - 560 O17, India
&
G. Prathap Jawaharlal Nehru Centre for Advanced Scientific Research, Indian Institute of Science Campus, Bangalore - 560 012, India
Higher-order plate theories and the corresponding finite element formulations have become popular recently, owing to the increasing demand for accuracy in transverse stress predictions in thick and laminated structures. In this paper, several displacement-type finite elements based on the Lo-Christensen-Wu higher-order description of transverse structural deformation are examined. The errors of locking, delayed convergence and stress oscillations are critically studied. The accuracy of the elements is then demonstrated using some wellestablished bench-mark tests. It is also observed with reference to a highprecision three-dimensional finite element solution that the higher-order transverse deformation theories and the corresponding finite element formulations, in general, cannot capture the transverse stress distribution correctly near the structural boundary of thick beams and plates.
1 INTRODUCTION
proper a m o u n t of transverse shear energy. In addition, the transverse stress distribution can be considerably inaccurate in the case of thick and laminated structures. As an i m p r o v e m e n t over the first-order shear deformable theories, Lo et aL 3,4 have presented a plate theory based on an assumed higher-order displacement field. This theory permits a quadratic variation of the transverse shearing strains and linear variation of transverse normal strain through the thickness direction. Eater, Reddy presented a third-order shear deformable plate theory, 5 which again neglects transverse normal strain effects, but assumes a parabolic distribution of transverse shearing strains through the thickness by enforcing zero shear conditions on the lateral surfaces a priori. As a result, the formulation needs C~-continuous description for the displacement fields in the finite element formulation. Finite elements based on C J-continuous displacement fields are not desirable for use in general purpose package applications. Liao et al. ~ p r o p o s e d a hybrid stress model based on the higher-order theory to improve the transverse dis-
A n u m b e r of theories exist in the literature for the analysis of thick and layered-composite structures. T h e classical theory based on the Love-Kirchhoff hypothesis that normals to the middle surface remain normal and unstretched after deformation was the starting point in the development of general plate theories. This theory neglects the transverse shear and transverse normal strain effects in the thickness direction. Today, the inaccuracy of such a hypothesis when applied to thick and layered plates is well known. Reissner I and Mindlin 2 provided first-order shear deformable theories based on the assumption that normals to the middle surface remain straight and unstretched but need not be normal to the deformed mid-surface. Such theories allow for a constant distribution of transverse shear strain through the thickness direction but still neglect the transverse normal strain. One drawback of such theories is that a fictitious shear correction factor is required in order to maintain the 445
446
P. R a m a Mohan, B. P. Naganarayana, G. Prathap
tribution of deformations and stresses using independent fields for the displacements u, v, w and the transverse shear stresses Oxz, Oyz. Again, the stress-based finite element formulations are not convenient for use from the point of view of general purpose finite element packages. The higher-order displacement field proposed by Lo et al. 3,4 requires only C°-continuity of the solution. This permits a great deal of flexibility in the specification of approximation functions over an element and the corresponding finite elements are best suited for general purpose applications. Unfortunately, most of these elements become too stiff when used to model thin structures and some elements become completely rigid as their span-to-thickness ratio becomes very large. This type of behavior is called artificial stiffening or locking.7,8 Much effort has been made to understand clearly and to eliminate the source of the errors of locking and the associated stress oscillations 7,~ especially in the case of finite elements developed using first-order shear deformable theories.9, ,~ Kant et al. ~ have presented a C°-continuous finite element formulation for the higher-order theory proposed by Lo et al. 3"4 A uniform reduced integration strategy is adopted for modeling the shear stiffnesss of the element, as an attempt to eliminate locking in a heuristic way. This paper examines critically the problems of shear locking and stress oscillations in the case of finite elements developed using the general higher-order displacement field proposed by Lo et al. 3'4 The sources of these errors are explicitly located and eliminated by applying the concepts of field consistency7"8 and variational correctness,2,13 in the case of two-noded and threenoded beam elements (B2 and B3) and the observations are then extended to the case of the four-noded and nine-noded plate elements (Q4 and Q9). The distribution of interlaminar shear and normal stresses across the thickness are evaluated at the centroid of the quadratic elements (B3 and Q9) by using equilibrium equations while the constitutive equations are used for recovering the transverse stress distributions in the case of the linear elements (B2 and Q4). Numerical differentiation is used for evaluating the second derivatives in the former case. Thus the one- and two-dimensional quadratic elements are able to give quartic variations of shear stress and cubic variation of the bending stress across the structural thickness.
The transverse stress distribution captured by such finite elements is then compared with a corresponding high-precision three-dimensional finite element solution. It is observed that the present finite elements are very accurate throughout the structural domain except near the boundary, particularly when it is clamped. In ,fact, none of the higher-order transverse deforrnation theories and the corresponding finite elements can recover accurate transverse stress distribution across the thickness near the structural boundary. This is because the transverse stress distribution needs a much higher order of description when compared to the quartic shear stress and cubic bending stress distribution (that can be generally achieved by the higher-order transverse deformation theories and the finite elements~ near the structural boundary. In addition, most of the higher-order theories assume zero shear stress at the free (top and bottom) surfaces of the beam/ plate everywhere. This is not true in the case of the clamped boundary where, in fact. the transverse stress distribution becomes singular near the top and bottom surfaces, which act as re-entrant corners. Thus. enforcement of quartic shear and cubic bending stress distribution and zero shear stress on the free surfaces lead to extraneous oscillations in the transverse stress distribution near the structural boundary. As a result, the higher-order transverse deformation theories and the corresponding finite elements, in general. become unreliable near the structural boundary.
2 FINITE ELEMENT FORMULATION In this paper, the higher order CU-displacement field description proposed by Lo et al. 3"4 is chosen for the finite element formulations. A family of layered-composite finite elements -- two-noded linear and three-noded quadratic beam elements CB2 and B3) and four-noded linear and ninenoded quadratic rectangular plate elements (Q4 and Q9) - is developed for a close examination. 2.1 Field description The displacement field proposed by Lo e t al. 3"4 consists of in-plane displacements expanded as cubic functions in the thickness coordinate. 3 The polynomial expansion for transverse displacement is truncated at one order lower the the expansion for in-plane displacements. It can be seen (see eqns (2e) and (2f) below) that this choice leads to a
Higher-order laminated plate theory consistent definition of interpolation for the transverse shear strains with reference to the thickness coordinate z. Here, the same assumptions are used in the finite element formulations such that the three-dimensional displacement field is defined using the mid-surface degrees of freedom u, 0,..., w* as follows:
The stress-strain relations for the L th lamina in the laminate coordinate axes system can be obtained by the usual transformation rules ~4 and are given by
OxxlL11012 [ 3Q16L[xqL
U(x, y, z)= u(x, y) + zO(x, y) + z2u*(x, y) + z 30*(x,y)
V( x, y, z) = v( x, y) + z qk(x, y) + z2v*( x, y) + z3 q~*(x, y)
(2a)
~ ~ , + z..3.h 6yy = V,v=V,v-Jr-Z~,v-Jr-Z2V,y ~,y,
(2b)
gxy= g,y-+- V,x~-(bl, yq-V,x) q- z( O,y-}- ~,x)
+z2(U,y+V,*)+z3(O,*+¢, *)
(2c)
6~ = W,: = ~p+ 2zw*
(2f)
Each layer in the laminate is assumed to be in a three-dimensional stress state so that the constitutive relationships for a typical lamina L with reference to the fiber-matrix coordinate axis (1-2-3) can be written as 14
111L { [clci2c30IL C12 C22 C23 0
L Q45 Q55 J [ 6xzj
dS.
(3)
The elemental equations can be obtained by minimizing the potential energy function FI with respect to the unknown nodal displacement vector. 2.2
T h e linear t w o - n o d e d b e a m e l e m e n t (B2)
The linear two-noded beam element with seven degrees of freedom per node is shown in Fig. l(a). It is assumed that V=0, U - - U ( x , z ) a n d W= W (x,z) such that the strain field given in eqn (2) reduces to:
.~622~
/
6xx
0
Q16 Q26 Q36 Q66
II = -2 va~jeodV-
1]),y)
+ z2(3~ * + w,,*).
0
OxyJ
(2e)
V,~+ W,y=(~b+ W,y)+ Z(2V* +
a12J
Q13 Q23 Q33 Q36
(2d)
+ z2(30 * + w,*)
022~ :
oz~/
in which {a,x, ayy, o~z, axy, oyz, o~z}and {6xx, 6yy, 6~z, e,y, eyz, 6~z} are the stress and strain components, respectively, with respect to laminate axes and Qij are the three-dimensional elastic constants in the laminate axes of the L th lamina. The potential energy function I1 of a laminated composite finite element of surface area S and volume V subjected to a transverse loading p(x, y) on a surface which is at a distance z = ¢ from the mid-surface can be expressed as follows
exz = U,~+ W,~=(0+ w,~)+z(2u*+ q),x)
6 yz
Q,z Q22 Q23 Q26
0-xz
W(x, y, z) = w(x, y) + Z~p(x, y) + Z 2 w*(x, y).
6x,:= U,, = u,~+zO,,+z:U*,x+Z30, *
]
Oyy} =
{'}L:rQ 4 4 ILI yzl L
(1)
The in-plane and transverse strain components associated with the displacement field in eqn (1) are given by:
447
0
C66
1612J
U,x=U,x+ZO, x+Z2U,*+z3¢, *
6 zz -~- W,z = ~ + 2zw*
(4a) (4b)
6 xz -~- U,~+ W,x=(0+ w,,)+z(2u*+ qJ, x)
0131
C55 [6131
where {01j , 022, 0-33, 0-12, 0-23, 0"13} are the stress and {ell, 622, 633, 612, 623, 613} a r e the strain components referred to the lamina coordinate axes and C 0 are the elastic constants of the L th lamina.
+z2(30"+
(4c)
Eyy= e,y = eyz = 0.
(4d)
The displacement components for the element are interpolated using linear Lagrangian interpolation functions so that, in an isoparametric formulation, we can express the active degrees of
P. Rama Mohan, B. P. Naganarayana, G. Prathap
448
I-
(doF)
I.
.I
~
= (u,O, ta",e',w,P,w'}
(a) 2 - n o d e d
2~--~ -
~,~,w,P,w" }
(b) 3-noded beam etement
bear~ eteMen±
3
4
=
g .I
4
7
3
9~
e
0
0
%
0
qJ
2
(do~') =
v,¢~ v~,~ ~,
w,p~W ~")
Wjp, W"~")
(c) 4-noded
~2
pl~te etenen±
(d) 9-noded pro.re
ete~en-t
Fig. 1. The family of higher-order elements: (a) two-noded beam element; (b) three-noded beam element; (c) four-noded beam plate element; (d) nine-noded plate element.
freedom as:
{} a°21 U*
/a21 a22
0*
laB, a32
{:} I::::::]{;l w*
(5)
L b2,622J
where the constants a0~,..., b22 are functions of the nodal degrees of freedom. Substituting eqns (5) in eqn (4c) we can express the transverse shear strain as
gxz={{all +b02/l}5-a12~} + z{ {2a2, + bt2/l} + 2a22~}
(6)
+ z2{{3a3, + b22/1} + 3a32 }. If this is written in terms of Legendre polynomials, ~°,~3 it can be shown that the transverse shear strain energy comprises squares of the coefficients of this expansion. It is then simple to identify the constraints on the transverse shear strain in the thin plate limits as
True all + bo2/l~O
Spurious a12--,0
2a21 + bt2/l--,O
a22~0
3a~ 1+ b22/l-,0
a32-, 0.
i7)
From eqns (7), we can argue that the true constraints give the correct Kirchhoff's condition at the element centroid, i.e. (ex=)~=0 ~ 0. On the other hand the spurious constraints reflect a (spurious) condition equivalent to Oe,,x/Oz~O. Thus the spurious constraints lock the element from bending and lead to highly underestimated displacements and bending moments. It can be shown that the locking error is directly proportional to (1/t) 2 as in the case of first-order elements. 7 Again from eqn (6), it can be observed that these spurious constraints lead to large linear oscillations in the shear stress predictions in the case of the inconsistent element which we identify as B2.0. The errors of locking and stress oscillations can be analytically derived in an a priori sense as was done for the corresponding first-order two-noded beam element, v,8 However, owing to the identical nature of the errors, the relevant derivations are not presented here. The preceding observations are confirmed later by a proper numerical example. From the method of Legendre polynomial expansion, ~°,~-~ the variationally correct way of eliminating the inconsistent terms, having associated with the spurious constraint, from the transverse shear strain field is to simply drop the inconsistent linear Legendre terms in ~ to eqn (6) such that the consistent shear strain field becomes,
r~={alt + boz/l} + z{2a2~ + b12/l} +z2{3a31 + b22/1}. It is demonstrated later, using appropriate numerical examples, that such a consistent element, B2.1, is free of locking and stress oscillations. B2.1 exhibits the variationally correct convergence rate and constant stresses over the element length.
2.3 The quadratic three-noded beam element (B3) The quadratic three-noded beam element with seven degrees of freedom per node is shown in Fig. l(b). The same assumptions that were employed in the case of B2 a r e used here ~also
Higher-order laminated plate theory such that the strain field given in eqns (4) remain valid. The displacement components for the element are interpolated using the quadratic Langrangian interpolation functions so that, in an isoparametric formulation, we can express the active degrees of freedom as:
{} Ia°la°a°]l =
u* O*
al 1 a12 a13 a21 a22 a23 a31 a32 a33
w*
~2
[_b2, b22 b23
(8)
ex== {a~ l + at3~3 + bo2/l} + {at2 + 2bo3/l}
+ Z{ {2a21 + 2a23/3 + bl2/l} + {2a22
+2b,3/l}~} + z2{ {3a3t + a33 + b22/l } + {3a32
+ 2b23/l}~}.
× ~-(a13/3)(1 - 3~ 2)
It is demonstrated later, using appropriate numerical examples, that such a consistent element, B3.1, is free of delayed convergence and stress oscillations.
+ z{ {2a21 + 2a23/3 + b12/l} + {2a22 + 2b13/l} × ~-(2a23/3)(1 - 3~ 2)}
× ~-- a33 ( 1 - 392)}.
inconsistent quadratic element (B3.0) is equivalent to that of the consistent linear element (B2.1). From eqn (9), it is clear that the spurious constraint is associated with the quadratic Legendre term and hence the inconsistent element suffers large quadratic oscillations while predicting the transverse shear stress. These observations are established using an appropriate numerical example later. From the method of Legendre polynomial expansion, the variationally correct way of eliminating the inconsistent terms from the transverse shear strain field is to simply drop the inconsistent quadratic Legendre terms in ~ from eqn (9) such that the consistent shear strain field becomes rxz = {a~, + a13/3 + bo2/l} + {a,e + 2bo3/I}~
~2
Substituting eqns (8) into eqn (4c) we can express the transverse shear strain as
-+- Z2{ {3a3, + a~3 + b22/l } + {3a32 +
449
2b23/1 } (9)
Again, from the method of Legendre polynomial expansion the constraint on the transverse shear strain energy leads to three discretized constraints in the limits of vanishing thickness. It can be observed that the constraints associated with the constant and linear Legendre terms are the true discretized constraints representing the Kirchhofff condition in the limits of vanishing thickness. On the other hand, the third constraint associated with the quadratic Legendre term is a spurious constraint on the linear part of the bending strain, i.e.
Thus, such an inconsistent element, B3.0, can capture only the constant bending moment part over the element length and the linear bending moment part will be highly underestimated. In other words, the inconsistent quadratic element, B3.0, does not lock as its linear counterpart B2.0 but has a lower convergence rate than the expected variationally correct rate. It can be observed that the rate of convergence of the
2.4 The linear four-noded plate element (Q4)
The four-noded plate element configuration is as shown in Fig. 1(c). The element displacement field is described by 11 degrees of freedom: u, 0, u*, 0", e, ~b, v*, ~b*, w, V2, w*. The regular bilinear Lagrangian shape functions are used for displacement field interpolation. It can be shown that the shear strain fields exz and ey= (eqns (2e) and (2f)) can be expressed as
Sx==Ao + Al ~ + A2rl + A3~rl e.,,~= B o + B~ ~ + B2rl+ B.~r I. For an undistorted (i.e. rectangular) element configuration, the terms associated with ~, s~r/in the shear strain field ex., (i.e. constants A~ and A3) are field-inconsistent and lead to spurious constraints in the penalty limits. Similarly, the constants B2 and B3 in the shear strain field % are field-inconsistent and lead to spurious constraints in the penalty limits. These spurious constraints lead to the conditions
ae,.,laz--, O,
Oe,>./Oz~ 0,
Oev/az-~ O.
Thus the inconsistent element, Q4.0, suffers from locking, highly underestimated bending moments
450
P. Rama Mohan, B. P. Naganarayana, G. Prathap
and large oscillations in shear stress predictions. Thus the shear strain fields have to be smoothed to their respective consistent form as Lz = A0 + fi'2 q and ryz =/~0 + B1 ~" Any variationally correct method of field-redistribution 1° can be used for the purpose. The fieldredistribution is equivalent to the use of 1 × 2 and 2 × 1 Gaussian quadrature rules for the shear strain energy components associated with the shear strain fields S~z and eyz, respectively, for a rectangular element. The field-consistent element, Q4.1, is free of the errors of locking and stress oscillations.
2.5 The quadratic nine-noded plate element (Q9) The element configuration for the nine-noded plate element is shown in Fig. l(d). Again the element displacement field description is based on 11 degrees of freedom. It can be shown that the terms associated with the quadratic Legendre terms ( 1 - 3~ 2) are field-inconsistent in exz while those associated with ( 1 - 3 r / 2 ) are field-inconsistent in eyz. The corresponding spurious constraints lead to the conditions of the type
Thus the linear variation in the bending moments is highly underestimated and as a result, the field-inconsistent element, Q9.0, exhibits delayed convergence, constant bending stresses and large quadratic oscillations in shear stress predictions similar to its first-order counterpart. ~3 The shear strain fields are expanded as Legendre polynomials on certain principal lines, and the inconsistent quadratic Legendre terms are dropped to get the consistent element version, Q9.1. The consistent element is free of delayed convergence and gives correct linear variation to the bending and shear stress fields as shown later in numerical examples.
2.6 Stress recovery The transverse stress field recovered from the constitutive relationships in a displacement-type finite element formulation is discontinuous at the interlaminar region of laminated structures. The equilibrium equations are often integrated through the thickness of the structure to give the
desirable continuous transverse stress distribution across the thickness. Such a procedure needs first derivatives of the normal stresses with reference to the in-plane coordinates; this cannot be done directly in the case of linear finite elements. In this paper, the stresses are calculated using the constitutive relationship for the linear elements (B2 and Q4) while equilibrium equations are used for predicting the shear stresses from the normal stresses that are recovered from the corresponding constitutive relationships in the case of the quadratic elements (B3 and Q9). All the stress components are recovered at the element centroids in the case of linear elements while the inplane stresses are sampled at the integration points corresponding to reduced Gaussian quadrature and the transverse stresses are sampled at the element centroids in case of quadratic elements. In the next section, the quadratic elements are used for solving some well-established benchmark tests, since they naturally can produce more accurate stress distribution across the structural thickness when compared to linear elements. These elements can sense quadratic, cubic and quartic variations for the transverse normal~ inplane normal and transverse shear stress components respectively. It is observed that the order of stress distribution mentioned above becomes insufficient near the structural boundary. In addition, the equilibrium equations are generally solved with an additional assumption that the shear stresses vanish on the top and bottom surfaces. This assumption is invalid particularly near the clamped boundary where, in fact. the stress distribution becomes singular near the top and bottom corners, which act as re-entrant corners. in the case of thick structures. Thus imposition of quartic distribution and vanishing transverse shear stress conditions on either the top and bottom surfaces leads to extraneous perturbations in the transverse shear distribution near a boundary as will be demonstrated later. But this behavior is inherent in general higher-order transverse deformation formulation and consequently they cannot recognize the stress singularity near the reentrant corners on the clamped boundary. However, the present element formulations are found to be very accurate away from the structural boundary when compared to a high-precision three-dimensional finite element solution. For comparison of the transverse shear stress distribution obtained from the C ° elements pre-
Higher-order laminated plate theory sented in this paper, a C 1-continuous theory proposed by Krishnamurthy ~5 is considered. This theory can recover constant shear and cubic bending stress distributions across the thickness near the clamped boundary; and quadratic shear and linear bending stress distributions across the thickness away from the clamped edges. It is observed later, that the transverse stress distribution is not accurate anywhere when compared to a corresponding high-precision three-dimensional finite element solution. ~6 However, the transverse shear stress distribution at the clamped edge in the interior region is most reliable when compared to other higher-order transverse deformation theories.
451
X
X
Theory/B3.1/Q9.1 B2.0/Q4.0 B3.O/Q9.0/ E32.1/Q4.1
*
+ x
co
'
'
'
I
'
3
'
'
I
'
'
'
I
'
'
5
4
II
'
I
'
6
'
'
t
'
7
'
I
8
No.EL
(a) * 1 element case + 2 element case
//*
3 NUMERICAL EXAMPLES In this section, we consider first an example to demonstrate the errors of locking and stress oscillation in field-inconsistent elements. Then we consider a problem to establish reliability of higher-order transverse deformation theories and finite element formulations in recovering accurate transverse stress distribution near structural boundary. Finally we consider a few well-established laminated plate problems to demonstrate the efficiency of the new field-consistent ninenoded plate element, Q9.1.
'
/ ~.~-
g--.@
e4-
3
4
5
o90/t5 (b)
6
7
8
Fig. 2. (a) Convergence rate with different consistent and inconsistent elements. (b) Locking error with reference to thickness ratio (l/h).
3.1 Cantilever slender beam under tip shear force A cantilever beam of length l = 10, breadth b = 1 and thickness t = 1-0 subjected to a tip shear force P = 1, is modeled with different consistent and inconsistent elements. Young's modulus of the material is 10 while the Poisson's ratio is zero. The rate of convergence with different elements is depicted in Fig. 2(a). The inconsistent linear elements (B2.0, Q4.0) lock severely as expected. In fact, the locking error is proportional to the square of the aspect ratio of the elements as in the case of their first-order counterparts (see Fig. 2(b)). The consistent linear elements show the variationally correct convergence rate (the discretization error for the present problem modeled with linear elements is given by 1/4n 2 where n is the number of elements used). It can also be noticed that the inconsistent quadratic elements (B3.0, Q9.0) do not lock but show a delayed rate of convergence that is equivalent to that of con-
sistent linear elements. The consistent quadratic elements (B3.1, Q9.1) give an exact solution even with a single-element model. The bending moment predictions from various element versions are plotted in Fig. 3. Again, as in Section 2.2, the inconsistent linear elements lock severely showing highly underestimated bending moments. The consistent linear elements show correct bending moments at the element centroids. On the other hand, the inconsistent quadratic elements cannot depict the linear variation of the bending moment in the penalty limits and as a result the bending moment is correct only at the element centroids. The consistent quadratic elements produce variationally correct linear bending moment variation over the element domain. Figure 4 shows shear force predictions from various elements. It can be observed that the
452
P. Rama Mohan, B. P. Naganarayana, G. Prathap j
£~
*
~o
+ X o
\
2
Theory/B3.1/Q9.1 B3.O/Q9.0 B2.1/Q4.1 B2.0/Q4.0
/ /
\\
/ "\
Fig. 5. \\ \
2
4
6
8
10
Bending moment distribution with different consistent and inconsistent elements.
* + x
Theory/B2.1/B3.1/ Q4.1/Q9,1 B2.0/Q4.0 B3.0/Q9.0
7/J
/
~< ,/
/
//
I 0
2
4
6
8
10
x
Fig. 4.
x
A deep cantilever beam under uniformly distributed shear body force.
\
x
~o. , ,* ' \
f
I ! I
i \
Ot--i
intensity in z-direction
:J \, \
~
body force of unit
/
\
Fig. 3.
Uniformly distributed
// \
0
"
J/"
Shear force distribution with different consistent and inconsistent elements.
inconsistent linear elements show large linear oscillations while the inconsistent quadratic elements show large quadratic oscillations while all the consistent element versions are able to predict the constant shear stress condition over the structural domain correctly. It should be noted that the behavior of different element versions is very similar to that of their first-order counterparts.
3.2 Cantilever deep beam under uniformly distributed shear load A deep cantilever beam of length 1= 10, breadth b= 1 and thickness t = 10 (Fig. 5) is subjected to uniformly distributed body force of unit intensity in the z-direction. The material properties are the same as in the previous example. The structure is
modeled by using consistent nine-noded plate elements (Qg.I) and the eight-noded consistent hexahedral elements (H8) of Ref. 16. A n adaptively graded mesh generated by the three-dimensional finite element package ~6 is used as the basis for both meshes. The finite element mesh is highly refined near the top and bottom corners at the root of the cantilever beam by using mesh gradings in the ratio of 1 : 2 : 4 : 8 : 1 6 : 3 2 : 6 4 : 1 2 8 and 1:1:4:8:8:4:1:1:1:1:4:8:8:4:1:t in the xand z-directions, respectively, in the case of H8. The same mesh grading is retained in the xdirection while 20 layers of equal thickness are used in the z-direction in the case of Q9.1. Stresses are sampled at the centroids of H8 elements while stress recovery is done at the interlaminar points on the centroidal locations of Q9.1 elements. The transverse shear and bending stress distributions across the structural thickness are sampled at the locations x=0-0196, x = 0 . 9 3 6 , x = 1.843 for both the models. The transverse shear stress distributions at different locations from Q9 and H8 models are plotted in Fig. 6 at locations x = 0"0196. x = 0" 196 and x = 1.843. respectively. The stress distributions derived from the classical beam theory and from another independent higher-order shear deformation theory ~5 are also plotted for comparison. As discussed in Section 2.5, Q9.1 shows very accurate quartic shear stress variation when compared to the three-dimensional solution (HS/ away from the clamped edge. However. near the clamped edge Q9.1 shows extraneous quartic oscillations such that the shear stress is not correct at any predetermined point across the thickness. On the other hand, many higher-order theories and the corresponding finite element formulations show quadratic shear stress variations (which are inaccurate when compared to the three-dimensional solution at all locations) across the thickness throughout the beam length. Krishnamurthy's theory provides a transverse shear
Higher-order laminated plate theory ---.... ---
Q9.1 HEXA8 EBT AVK
\,,\ to
N
]
4
/
JJtl _..P
o , , ~,'?-7-,--T--,__,-:~,-q 0
5
10
~z
15
20
Ca) --
Q9.1 HEXA8
-
N
t
I
_b. o
,
,
,
i
,
,
,~-~
~ ,
,
4
2
,
i
6
8
,
,
,
i 10
,
;
,
i 12
,
,
; 1 14
4 53
stress distribution which is constant at the clamped edges and gradually becomes quadratic away from them, as shown in Fig. 6. The accuracy of the shear stress distribution obtained from Krishnamurthy's theory is very good in the interior region at the clamped edge when compared to the three-dimensional solution. However, the stress distribution does not compare well with the three-dimensional solution away from the clamped edge as shown in Fig. 6(c). Again, this cannot recognize the stress singularities at the reentrant corners on the clamped edge. The shear stress distribution at z = 9-91 along the axis of the beam is obtained using the threedimensional finite element mesh and is plotted in Fig. 7. The three-dimensional finite element solution clearly shows the singular nature of the stress near the clamped edge. However, as a rule most of the higher-order finite element formulations cannot depict the stress singularity near the clamped edge, as shown in Fig. 6. All the higher-order solutions can recover the cubic variation of the bending stress across the beam thickness accurately. However, the bending stress becomes singular near the top and bottom surfaces at the clamped edge. Again, the threedimensional solution can capture this nature correctly while all the beam/plate theories fail to do the same. Krishnamurthy's theory is unreliable even at locations away from the clamped edge giving a linear variation to the bending stress, which is actually cubic in nature, over the thickness. However, owing to the similar behavior
CF-X Z
(b) o
o
¢o
----
_
EBT (O-xx) HEXA8( 0;',x) HEXA8 ( d~xz)
o N
----
A8
o
b...... o
, , 1 1 1 , , [ , , , i , i 2
4
, j , , , ] , l , ] , , i 6
8
10
12
\,
j 14
%z (c) Fig. 6. Transverse shear stress distribution for a deep cantilever beam: (a) x = 0"0196; (b) x = 0.1960; (c) x = 1.8432.
5
10
15
20
x Fig. 7.
Bending and shear stress distribution along the axial direction at z = 9"91.
P. Rama Mohan, B. P. Naganarayana, G. Prathap
454
when compared to transverse shear stress distribution, the results are not shown here. The transverse stress distributions recovered from the higher-order transverse deformation elements are perturbed near the free edge also. Since the magnitude of the stresses is infinitesimally small near a free edge, the impact of the inaccuracy is negligible. Again, the results are not shown here.
solutions 17 and Reddy's solutions? From the results, one can notice that the present solutions are in close agreement with Pagano's solutions compared to Reddy's solution, especially the transverse shear stress. This may be attributed to the fact that imposition of any additional constraints in a displacement approach, as in Ref. 5. reduces the flexibility of the model to a certain extent.
3.3 Simply supported 0/90/0 rectangular plate under sinusoidal load
3.4 Simply supported 0/90/0 square plate under sinusoidal load
A simply supported, three-ply (0/90/0) rectangular laminate with length-to-breadth ratio a/b--3 and layers of equal thickness, subjected to sinusoidally distributed transverse load is considered for comparison of maximum deflection and interlaminar stresses. The material properties used are given by
A simply supported, square 0/90/0 laminated plate with a length-to-thickness ratio (a/h) equal to 4 subjected to sinusoidal distributed transverse load on the top surface is considered for validation of interlaminar stresses. Each layer of the laminate is of equal thickness and the material constants with respect to its principal material axes are specified as: E ~ = 2 5 × 1 0 6 psi, Ee2=E33=106 psi, GI2=GI3=0"5×lO ~ psi, G23 = 0"2 × 106 psi, vl2 = v23 = v~ = 0"25. For comparison the following normalized quantities are used in the study: {6x~, d~ = (1/pos) {oxz, aye}; dz=(1/po s2) Oz where s = a/h and P0 is the maximum load intensity at the plate centroid. Figure 8 shows the distribution of 6x:, 6~: and 6. across the thickness and is compared with the solution of Pagano,17 and Liao et al. 6 From Fig. 8 it can be observed that the present solutions are in close agreement with Pagano's elasticity solution, especially the distributions of 6xz and 6~. Here, 6: is sampled at the plate centroid (a/2,a/2) while 6~= and 6y; are sampled at the mid-points of the plate edges y= a/2 and x = a/2, respectively.
E I / E 2 = 2 5 , GI2/E 2 = 0"5, G23/E2 = 0.2,
E2 = E3, G13 = Glz, '1,'12= 1,'23= "V13-----0"25.
The non-dimensionalized deflections vP and interlaminar shear stresses dxz, 6yz, are given by: v~=100 (E2h3/qoa4)w; {6~,d~zt=(h/qo a) {axe, aye} where a is the edge length of the square plate and h is its thickness, and are presented in Table 1. Here, • is sampled at the plate centroid while 6~, and 6y~ are sampled at the midpoints of the plate edges y= b/2 and x = a/2, respectively, at z=0. Results obtained for different thickness ratios (a/h) are compared with Pagano's elasticity
3.5 Simply supported 0/90 square plate under sinusoidal load Table 1. Non-dimensionalized deflections and stresses in a rectangular laminate under sinusoidal load (0/90/0, h~= h/ 3, b/a = 3)
a/h
Source
~
b~=
6v~
4
Q9.I Ref. 17 Ref. 5 Q9.1 Ref. 17 Ref. 5 Q9.1 Ref. 17 Ref. 5 Q9.1 Ref. 17 Ref. 5
2'626 2.820 2-641 0.867 0.919 0.862 0.595 0.610 0.594 0.507 0.508 0.507
0-3718 0.3870 0"2724 0"4220 0.4200 0.2859 0.4314 0-4340 0.2880 0.4344 0.4390 0,2885
0"0317 0"0334 0.0348 0.0146 0.0152 0.0170 0-0116 0.0119 0.0139 0.0106 0.0108 0.0129
10 20 100
Finally, a simply supported square cross-ply, antisymmetric (0/90) laminated plate under sinusoidal transverse load is considered for comparison of deflections and transverse stresses. The set of material properties used is the same as that specified in Section 3.3 The non-dimensionalized deflections w and interlaminar shear stresses 6x., 6yz, with ~ = 10 (E2h3/qo a4) wand {oxz, 6yZ}=( h/qo a) {oxz, oy~}are presented in Table 2. Here, ~ is sampled at the plate centroid while ox~ and oy~ are sampled at the mid-points of the plate edges y= b/2, z = 0-2 h and x= a/2, z = - 0"2 h, respectively. The results for maximum stresses and deflection are corn-
Higher-order laminatedplate theory -
-
.
.
.
.
.
.
-
-
Table 2. Maximum deflection and transverse stresses for a simply supported unsyrnmetric cross-ply (0/90) square plate under sinusoidal load
Q9.1 Pagano Liao HOPE
-
.
e4
¢5
5
a/h
Source
v~
G~
6yz
4
Q9.1 Ref. 17 Ref. 18 Q9.1 Ref. 17 Ref. 18
0"2152
0'3091 0'3127 0.2868 0'3318 0.3310 0.2957
0"3028 0.3188 0.2865 0'3315 -0.2956
10
LI c
¢-1
455
-
-
0'2020 0"1264 -0.1220
/5- - / I 0.2
0.0
0,4
--
~o
0.6
d
(a)
-
(25-
~
-
Q9.1 T.Kant Pagano
-
~.
c-4
o-
Q9.1 PagGno LiGo
--
~/
.....
I ....
/
HOPE
r ' / ~ l ' ' ' ' r
....
Gz
I ' ' Frq
.c N
o d
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
Fig. 9. Simply supported 0/90 square plate under sinusoidal load. Transverse distribution of normalized shear stress dx~.
0.05
0.00
0.10
0.15
0.20
0.25
(b)
Q9.1 Pagano Liao/HOPE
..... ,*
/
)j~//-
/
,:::5
o c~
~,'1
....
i'''t
....
I ....
i ....
I ....
i ....
I
c5 I
¢5 ..// 0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
pared with Pagano's three-dimensional elasticity solutions and are presented for different a/h ratios. T h e distribution of 6, z, through the thickness is shown in Fig. 9 and are c o m p a r e d with three-dimensional elasticity ~7 and Kant's ~8 solutions. T h e present solutions are close to the threedimensional elasticity solutions of Kant. j8 This m a y be attributed to fast convergence of results because of consistent shear strain field definitions which is equivalent to the use of 2 x 3 and 3 x 2 Gaussian rules for the shear energy c o m p o n e n t s associated with the shear strain fields exz and eyz, respectively. It should be noted that a uniformly r e d u c e d integration rule (2 x 2) is applied for the whole shear strain energy in Ref. 18 which, in turn, results in spurious zero energy m o d e s in the case of the n i n e - n o d e d plate element.~3
(c) Fig. 8. Simply supported 0/90/0 square plate under sinusoidal load: (a) transverse distribution of normalized shear stress 6~:; (b) transverse distribution of normalized shear stress 6~; (c) transverse distribution of normalized normal stress 6 z.
4 CONCLUSIONS In this paper, the errors related to the finite elem e n t formulations based on the C°-continuous
456
P. Rama Mohan, B. P. Naganarayana, G. Prathap
higher-order description of the displacement 3 are critically examined. The cause for the errors of locking, delayed convergence and stress oscillations present in a conventional finite element formulation corresponding to such a higher-order transverse deformation theory has been explicitly identified and eliminated for linear and quadratic beam elements. The consistency and correctness rules are then extended to formulate linear and quadratic laminated plate elements which are free from the errors of locking and stress oscillations. It is demonstrated through numerical examples that the errors for such higher-order finite element formulations follow very similar patterns as in the case of the corresponding first-order finite element formulations. It is also observed that the stress distribution across the structural thickness from any higher-order theory is not reliable near the structural boundary. Particularly these theories and the corresponding finite element formulations cannot depict the singularities in transverse stress distribution that can arise due to the boundary conditions. The present consistent higherorder transverse deformation finite elements are shown to be highly accurate in predicting transverse stress distribution in some well-established laminated plate problems. So far, discussions have been confined to rectangular plate elements. It is highly desirable to have higher-order plate elements which can assume an arbitrary quadrilateral planform. Work is now under way to apply the edge-consistency corrections 19 (which were very effective for the first-order plate elements) to the higher-order plate elements as well.
ACKNOWLEDGEMENTS The authors are grateful to Prof. R. Narasimha, Director and Dr B. R. Somashekar, Head, Structural Sciences Division of the National Aeronautical Laboratory, Bangalore, India for their constant encouragement and keen interest in the subject. Thanks are also due to Dr Kota Harinarayan, Director, Aeronautical Development Agency for permission to use the computer facili' ties at his organization.
REFERENCES 1. Reissner. E., The effect of transverse shear deformation on the bending of elastic plates. ASME J. Appl. Mech. 12 1945 A69-A77. 2. Mindlin. R. D.. Influence of rotatory inertia and shear deformation on flexural motions of lsotropic elastic plates. ASME J. Appl. Mech.. 18 ( 1951 31-8. 3. Lo. K. H., Christensen. R. M. & Wu. E. M.. A high-order theory of plate deformation, part i: homogeneous plates. ASMEJ. Appl. Mech.. 44 19771663-8. 4. Lo. K. H.. Christensen. R. M. & Wu, E. M.. A high-order theory of plate deformation, part 2: laminated plates. ASME J. Appl. Mech.. 44 t 1977) 669-76. 5. Red@. J. N.. A simple higher order theory for laminated composite plates. J. Appl. Mech.. 51 ( 1984 745-52. 6. IAao. M. L.. Jing, H. S. & Hwang, M., Improvements on the higher order plate element with partial hybrid stress model. Cornput. Struct.. 42 l) ( 1992 45-51. 7. Prathap, G. & Bhashyam. G R.. Reduced imegration and the shear-flexible beam elements. Int. J. Nurner. Meths. Engng, 18 (1982) 195-210. 8. Prathap, G.. Field-consistency -- toward a science of constrained multi-strain-field finite elemenl formulauons. S A D H A N A . Proc Engng Sci. lndian Acad. Sci.. 9 1986) 345-53. 9. Ramesh Babu. C.. Field-consistency in the finite element formulation of multi-strain-field problems in structural mechanics. PhD thesis. IIT. Madras, 1986. 10. Naganarayana, B. R. Consistency and correctness principles in quadratic displacement type finite elements. PhD thesis. Indian Institute of Science. Bangalore. 199 I. 11. Kant. T.. Owen. D. R. J. & Zienkiewiez. O. C.. A refined higher-order C t' plate bending element. Comlmt. Struct., 15,11982) 177-83. 12. Prathap. G.. A variational basis for least-squares fieldredistribution of strain functions in the finite element formulation of constrained media elasticity. TM ST 8801. N A L Bangaiore. India. 1988. 13. Naganarayana. B. P.. Prathap, G.. Dattaguru. B. & Ramamurthy, T. S.. A field-consistent and variationally correct representation of transverse shear stratus in the 9-noded plate element. (~mgmt. Meths. Appl. Mech. Engng, 97 374 (1992) 355-74. 14. Jones, R. M., Mechanics of Composite Materials. McGraw-Hill, 1975. 15. Krishnamurthy, A. V., Toward a consistent beam theory. AIAA J., 22 (6)(1984) 811-16. 16. Prathap, G. & Naganarayana, [3. P., 3D-FEES -- lheoretical Manual, PD ST 9005, N A L Bangalore, India, 1990. 17. Pagano. N. J., Exact solutions /or rectangular bidirectional composite and sandwich plates. J. Cornp, Mater., 4 (1970) 20-34. 18. Kant, T. & Pandya, B. N., Finite clement stress analysis of unsymmetrically laminated composite plates based on a refined higher-order theory. Proc. Int. Conf. Composite Materials and Structures, 6 - 9 January 1988, Indian Institute of Technology, Madras, India. 19. Prathap, G. & Somashekar, B. R., Field- and edgeconsistency synthesis of a 4-noded quadrilateral plate bending element. Int. J. Numer. Meths. Engng, 26 (1986) 1693-708.