cw- 7919/ 92/ 0201%o7so3.

aw
@ 1992 Papnon Pms L. td.
Compua~ b Smcrwer Vol . 15. No. 2. ~. 177483, 1982
F%wdi nGrerl Bri l ai n.
A REFINED HIGHER-ORDER Co
PLATE BENDING ELEMENT
TARUN KANT+
Department of Civil Engineering, Indian Institute of Technology, Bombay 400 076. India
and
D. R. J. OWENS and 0. C. ZIENKIEWICZQ
Department of Civil Engineering, University of Wales, University College of Swansea. Swansea SA2 8PP. Wales
(Received 2 December 1980; received for publication 20 lanuory 1981)
Abstract-A general finite element formulation for plate bending problem based on a higher-order displacement
model and a threedimensional state of stress and strain is attempted. The theory incorporates linear and quadratic
variations of transverse normal strain and transverse shearing strains and stresses respectively through the
thickness of the plate. The 9-noded quadrilateral from the family of two dimensional Co continuous isoparametric
elements is then introduced and its performance is evaluated for a wide range of plates under uniformly distributed
load and with different support conditions and ranging from very thick to extremely thin situations. The effect of
full. reduced and selective integration schemes on the final numerical result is examined. The behaviour of this
element with the present formulation is seen to be excellent under all the three integration schemes.
INTRODUCTIOh
Isoparametric plate elements [ I-51 based on Mindlins
theory(61. which obviate most of the problems that beset
elements based on the classical Kirchhoffs plate
theory[7]. have become popular and well established.
However the theory itself has certain limitations, viz.,
the effects of transverse normal strain and transverse
normal stress are not accounted for. the transverse
shearing strains are assumed constant through the
thickness and a fictitious shear correction coefficient is
introduced to account for warping of the cross-section.
Thus the claim that the element is applicable to both thin
and thick plates is not tenable. At the most only
moderately thick plates could be analysed. In this paper
a Co continuity element is developed based on a higher-
order displacement model than hitherto used and a three
dimensional state of stress and strainI8.91. Specifically,
the transverse normal strain and the transverse shearing
strains are now assumed to vary linearly and quadratic-
ally respectively through the thickness of the plate. Thus
the warping of the cross-section is automatically in-
corporated. This element retains the simplicity of
Mindlins fermulation and. at the same time, has no
intrinsic limitations. It is expected that this development
will offer refined and realistic solutions in the case of
really thick situations in general, and composite/layered
systems, stress-distribution near concentrated loads, etc.
in particular. The element which has six degrees of
freedom per node can employ the available family of the
Serendipity and the Lagrange elements[S, IO]. In the
present study a 9-noded Lagrange .quadrilateral is used
which has performed very well on all of the test problems.
Numerical results corresponding to exact, reduced and
*Assistant Professor. presently a Visitor in the Deoartment of
Civil Engineering. University -of Wales, University College
Swansea. Swansea SA2.8PP. Wales. under the Jawaharlal Nehru
Memorial Trust (U.K.1 Scholarship 1979 Award.
*Reader.
#Professor and Head of Department.
OUfx, y, z), 6,(x, y), w*(x, y). etc. will henceforth be written
simply as 1. 6,. a+*. etc.
selective integrations for a wide spectrum of plates-very
thick to extremelythin-are presented.
ELEMENT FORMULATION
The general theoryll, 91 is based on the displacement
model (Fig. 1)
Lltx. y, I)= z&(x, y)+ z3B$tx, y)
u =
L
vtx. y. zl = 20,(X, y) + z??:tx, y)
wtx. !, 21 = wtx, y) + z2w*tx, y)
1
ma
where U. 1 and W define the displacement components
in the coordinate directions x, y and z respectively. Both
the x and y coordinates lie in the reference plane of the
plate which is assumed unstrained. The terms 8, and 6,
are the usual average rotations of the normals to the
reference plane along the x and y directions respectively,
while w defines the transverse displacement at the
reference surface. The terms @, 8: and w* are the
corresponding higher-order terms in the Taylors series
expansion used in the present theory and are also defined
at the reference plane. Thus the generalized displace-
t
z
W
8,
V
)_
c
Y
J
X
Fig. I. Positive set of displacement components.
177
178 T. KANT et ai.
ment of the reference plane S is expressed in terms of
six independent variables,
s = [w, RX, 0,. w*, et, 0:). (2)
The six three dimensional strain components are then
related to the generaiised displacement components 6 by
the relation,
E
1[
=-K +zK*
. I
Q, =zK, +zK;
Y
yxy = zK,, t z'K:,
Yxz = d, +z24$,
yy* = d, + Z'bC
E,=zzw*
(31
where the new curvature terms K,, KY, etc. called the
generalised strain components vector Q in the two
dimensionai finite element context, are related to the
generalised displacement components 6 by the following
matrix relation.
o~oooo-
0$
0 0 0
Od&$O 0 0
~10000
~01000
0000~0
A
00000;;;
oooo~-&
000;30
000;03
oooioo
= L6 (4)
The total potential energy 1~ for the present theory[9)
is given by
rrudA-
j-/z++p.-l(w+;w*)dA
(5)
where pz+ and pL- are the transverse distributed loads on
the positive and negative extreme z planes respectively
and h is the total thickness of the plate. The generalised
stress component vector (I which is the integral of the
physical stress components through the thickness of the
plate (see Ref. (91 and Figs. 2-4) is given by
v = IMP. M,. M,,, Qx, Q,. ML M:, Id:,, Q:.Q:, MAT.
(6)
The generalized stress vector u and the generalised
strain vector E are partitioned as follows:
and
where
oh! = IM,, M,, MrvIT. g,,z = [M:, M:, M:,],
@1,3=Mz,a,, =[Qx, Qy17.u,z=[Q:.Qtl'.
+,t = [Kc. Ky, Kxy]. e,,z =[KC. KF. K:,],
43 - - 2w*. f,I = [&. &I. e.2 = [bT, C#Jt]. (9)
For a linear elastic material the constitutive relation can
be written as
ghl = &,,rb, + &)hZChZ + Qb)h3 zw* 1 IOal
ah2 = &E~, f Ddh2 f & Zw c lob,
ghh) = 4,%,,i - D&h2 + &, ?W* f IOCJ
~~I=QslE,i~D,?E,: ClOdl
ff,: = &%,I + 4&? 1 IOel
where the elasticity matrices 4,. Q,,*. &,, pVz. etc. for a
homogeneous and isotropic plate with Youngs modulus
E. Poissons ratio v and thickness h are expressed in the
following manner:
&z= Eh5
8O(l + v)(l -2~)
?h)h4 =4,8(1
Eh
+V)(l-ZV)
I)h.(= Eh5
8Ofl + u)(l -2~)
Eh
&= 12(1 +v)(l-zV)
(I - VI
(IIf)
( I Ig)
(llhl
Eh IO
DS, =-
i 1
16Otl+v) 0 I
Illil
IL -1 0
0
0 I -2v
,
IL
0 I
h 1
0
-?v
--i-- -1
L 0
I-o
0
--J 0
I-2V
2
(I la1
illb)
I I IC)
(IId)
IlIe,
A refined higher-order Co plate bending element 179
!?_
2
I/
F
I!
Y
2
\Reference
surface
Fig. 2. Positive set of stress components.
/ / ,i (P,*i,) / / / /
Fig. 3. Positive set of stress resultants-forces
dy
CAS Vol IS. No. 2-F
Fig. 4. Positive set of stress resultants-couples
180 T. KANT et al.
It is seen that the strain energy expression of eqn (5)
contains only the first derivatives of the components of
the generalized displacement vector 6 and thus only Co
continuity is required for the shape functions to be used
in the element formulation. If the same shape function is
used to define all the components of the generalized
displacement vector S, then
6 = f NiSi
i-l
(12)
where N, is the shape function associated with node i, &
is the value of 6 corresponding to node i and m is the
number of nodes in the element.
With the generalized displacement vector 6 known at
all points within the element. the generalized strain vec-
tor Q at any point is determined with the aid of eqns (4)
and (12). as follows:
where
Bi = LN,
(14)
and is given explicitly for the present case as follows:
,
0~0000
oo;ooo
o;;ooo
~l0000
~01000
0000~0
00000~
OOOOLL
dy b
oooi30
000;03
000200
V, (IS)
Having obtained the D and 8 matrices as given by
eqns (1 I) and (IS) respectively the element stiffness
matrix 0. can be conveniently expressed in the follow-
ing form [ IO]:
K= BTQB dA =
BTDBIJI d5 d?
(16)
For the present element which has six degrees of
freedom per node as opposed to Mindlins element[I, IO]
which has only three degrees of freedom per node, it is
considered appropriate to economise on the cost of
computation involved in numerical quadrature for
evaluating the integral in eqn (16). It is generally agreed
that the foregoing computation accounts for a very sub-
stantial part of the total computation cost of the com-
plete analysis[lll. To achieve this end, the element
stiffness matrix is partitioned in both the row and column
directions into blocks. Each block corresponds to the
degrees of freedom at a node[12]. The expression of the
stiffness coefficients for a typical block ij could be writ-
ten as follows:
K = .i, 2,
WP w, IIIBiTm,
(17)
where W, and W, are weighting coefficients, n is the
number of numerical quadrature points in each direction
and & and & are the strain-displacement matrices based
on the ith and jth shape functions, respectively. Due to
symmetry of the stiffness matrix K, only the blocks K"
lying on one side of the main diagonal are formed (Fig.
5). Further, the suggestion of Gupta[l3] is incorporated
in the program which stipulates explicit multiplication of
the Bi, D and tj matrices instead of carrying out the full
matrix multiphcation of the triple product. This is
expected to reduce the computing time for element
stiffness formulation by a factor of nine. For the con-
venience of any prospective user, these 6 x 6 submatrices
corresponding to flexure and transverse shear effects are
presented in Appendix I.
The consistent load vector P due to distributed load p
can be written in general as:
P=
I
+I t,
NTpdA =
A IJ
Nrp//l dl dq
(18)
-I -I
which takes the following form suitable for numerical
integration when transformed in the context of the
energy expression (5) for the theory presented in this
paper
I
1
:
p=pi, J, w~w~lzWi p14 (A+ f P*-). (19)
S _
0
I
0
It is noted that eqn (19) is different from the usual ones
because the loads are not assumed to act on the
reference surface in the present development.
I ,
I
$* I Kz3 ;_-_--_________J K2m
1 - I-
I I
I
Symmetrical ,
I
\
I
\ I
\
\
\
cl_
K
mm
m: number of nodes
1
Fig. 5. Nodal partitions of the element stiffness matrix.
A refined higher-order Co plate bending element 181
A computer program was developed based on the
foregoing approach to generate the element stiffness
matrices, the consistent load vector and a few related
processes for the stress compu~tion. This program was
integrated with the modified version of the one described
in detail by Hinton and Owen[l4].
NtlME3UCAl E-
A square plate of side a = 1 and subjected to a trans-
verse distributed load p2 = 1 is analysed with three types
of ~undary conditions. A Poissons ratio (v) of 0.3 is
assumed. The numerical experimentation is carried out
with four of the nine noded Lagrangian quadrilateral
elements in a quarter plate. All the computations were
performed on a CDC 7600 computer in single precision
with 14 significant digits of real rounded arithmetic ac-
curacy. Full. reduced and sefective integration schemes
based on gauss-Legendre product rules, viz. 3 x 3 and
3 x 3f.E). 2 x 2 and ? x2(R) and 3 x 3 and 2 x 2(S), have
been employed for flexure and shear contributions res-
pectively to the element stiffness matrix computation.
The results for the plate analysed with the three boun-
dary conditions are presented in Tables 1-3 for various
values of the aspect ratio u/h. Typical execution time is
0.52 CP seconds per analysis.
CONCLUSIONS
The performance of the 9-noded Lagrangian iso-
parametric element has been studied in conjunction with
the refined higher-order plate theory presented in this
paper-and also elsewhere[8,9]. On the basis of the
results of Tables l-3 it can be positively concluded that
the elements buoyance is excelient in all the situa-
tions studiid in this paper with no apparent preferential
choice for full (exact), reduced and selective integration
schemes for computation of the element stiffness
matrix[l-5,15,16]. This appears to be very significant
and perhaps is a pointer towards the use of the present
relined higher-order theory in preference to Mindlin
theory hi~erto used in most of the analyses. Qualita-
tively this could be due to the better represen~tions of
the cross-sectional deformations and the stress-strain
law. However, its behaviour in context of the few boun-
dary constraints which lead to the development of
mechanisms (or near mechanisms) commonly interpreted
as zero energy modesf5f needs to be seen.
REFERENCES
E. Hinton. A. Razzaque. 0. C. Zienkiewicz and J. D. Davies.
Simple finite element solution for plates of homogeneous,
sandwich and cellular construction. Proc. of theInsfitutian
of Ciuil Engineers, Pan Ii. 59.4365 (1975).
T, J. R. Hughes. R. L. Taylor and W. Kanoknukulchai. A
simple and efficient finite element for plate bending. inr. j.
Num. Me&. Engng. 11, 15294543 (1977).
E. D. L. Pugh E. Hinton and 0. C. Zicnkiewicz. A study of
quadrilateral plate bending elements with reduced in-
tegration. Int. 1. Num. MeGr. Engng. 12. 10594079 (19781.
Table I. Simply supported square plate under uniform loading with w = @, = W* = f?: = t&l along its four edges
a/h 1ntes.
centI- Central+ tWll.r+ flid-edge+
TYW
01s01acemant Bendin Ploment
s
T\,lst ing l%rnent Shear Force
P&?/O Pa
Pa
Pa
E 0.00653 0.0521 - 0.0223 0.311
2 R o.GG054 0.0509 - 0.0216 0.285
S O.OOB5~ 0.0522 - 0.0219 0.285
E 0.00612 0.0499 - 0.0276 0.315
3 Fi 0.00612 0.0488 - 0.0262 0.285
s 0.00611 0.0500 - 0.0272 0.265
E 0.00522 0.0491 - 0.0295 C.322
4 R 0.00523 O.D4RD - 0.0273 0.285
s 0.00523 0.0492 - 0.0292 0.285
E O.OG480 0.0467 - 0.0304 0.330
5 R 0.00482 0.0477 - 0.0267 0.285
S 0.00481 0.0488 - 0.0302 0.285
E 0.00423 0.0480 - 0.0315 0.381
IO R 0.00426 0.0472 - 0.0299 0.285
s 0.00426 0.0484 - 0.0314 0.2&S
10.00427 Exact;
E 0.00398 0.0470 - 0.0317 0.573
5c R G .00408 0.0471 - 0.0302 G.ZBS
5 0.00408 0.0482 - 0.0319 0.285
E 0.00397 0.0469 - 0.0317 0.597
1OC R 0.00407 0.0471 - 0.0302 0.205
5 0.00407 0.0482 - 0.0319 0.264
E 0.05396 0.0468 - 0.0317 0.606
1OOG R 0.00407 0.0473 - 0.0302 C.285
S c .0040? 0.0482 - 0.0319 0.234
E 0.0039E 0.0468 - 0.0317 O.SOE
10000 R C.OO4G7 0.0471 - 0.0302 0.285
s 0.00407 0.0482 - 0.0319 0.2P4
I I*+ 0.004OE 0.0479 - 0.0325 0.338
1
Classic+? thin piate solution [Fief. 71
T. KANT et ai.
Table 2. Simply supported square plate under uniform ioading with w = W* = &s*t along its four edges
a,h In-%.
Central cmtr41+ corner+ Ri d- edBe+
Type
01sp1acsmnt Bendi ng mmmt Twi sti n FI omant
4
Shosr FOTCD
' pa' / 0 *Da* *XI *Pa
E 0.00949 0.0602 - 0.00206 0.324
2 R 0.00944 0.0596 - Q.00302 0.30t
S 0.00943 0.0605 - 0.00095 0.301
E 0.00605 0.0573 - 0.00509 0.336
3 R 0. 00893 0. 0566 - o.oo545 0.310
S 0.00693 0.0577 - 0.00337 0.309
E 0.00592 5.0552 - 0.00760 0.345
4 R 0.00593 0.0548 - 0.00751 0.314
S 0.00592 0.0556 - 0.00536 0.312
E 0.00536 0.0536 - 0.01051 0.353
5 R 0.00541 0. 0534 - 0. 00026 0. 316
s 0. 00540 0. 0546 - 0. 00716 0. 313
E 0. 00446 0. 0503 - 0. 02Oc 0. 397
10 A 0. 00456 5. 0503 - 0. 0152 0. 322
S 0. 30456 Cl . 0521 - 3.0135 0. 316
E o. I 30400 3. 0472 - 0. 0306 0. 561
50 R 0. 00419 0. 0483 - 0.0218 0. 341
S 0. 00420 0. 0511 - a. 0216 0. 339
E 0.00307 0.0469 - 0. 0314 0. 563
100 A 0. 00416 I l . 0482 - 0. 0221 0. 343
s 0. 00416 0. 0511 - 0.0221 0. 341
E 0. 00396 0. 0466 - 0. 0317 0. 591
1000 R 0. 00417 0. 0461 - 0.0223 0. 344
S 0. 004' 17 0. 0511 - 0. 0223 0. 342
E 0.00396 0. 0466 - 0. 0317 3. 591
10000 R 0. 00417 3. 046! - 3.0223 3. 344
5 0. 00417 G. 0511 - 0. 0223 0. 342
+
val ues U. neareat Gauss poi nt.
Table 3. Clamped square plate under uniform loading with w = l+ = 8, = w* = 8: = 6: = o(c) along its four edges
a/h
Intag.
Central C&WC ilid-Edge+ nid-ea*~+
: YW
Ofsplacment Bendi ng mnent Bendi n i %moot
f
Shssr Fores
- pa' / D - 5. 2
*Pa l Pa
E 0. 00607 0. 0335 - 0. 0262 0. 332 . -
2 R 0. 00611 0. 0320 - 0.0100 0. 297
s J . 00609 0. 0338 - 0. 0267 0. 297
E 0. 00349 0. 3263 - 0. 0300 0. 352
3 R 0. 00355 0. 0270 - 0. 0242 0. 306
S 0. 00352 0. 0267 - 0. 0319 0. 306
E 0. 50253 0. 3261 - 0. 0320 0. 374
4 R C. 00261 0. 0252 - 0. 0264 0. 315
5 0. 00256 0. 0267 - 0. 0344 0. 315
E G. 00206 0. 0250 - 0. 0337 0. 396
5 R 0. 00216 0. 0243 - 0. 0277 0. 323
s 0. 00211 0. 0256 - 0.0350 0. 323
6 0. 00136 D*O226 - o. bw 0. 516
10 R 0. 00152 0. 0226 - 0.0299 0. 344
5 0. 00146 0. 0236 - 0. 0387 0. 344
E 0.001w 0. 0207 - 0. 0194 0. 906
so R 0. 00129 0. 0222 - 0. 0308 0. 356
S 0. 00123 0. 0227 - 0. 0309 0. 356
E 0. 00097 0. 0205 _ 0. 0177 0. 956
100 R 0. 00126 0. 0222 - 0. 0308 0. 357
S 0. 00123 0. 0227 - 0. 0390 0. 356
E 0. 00096 0. 0204
1000 R 0. 00126 0. 0222
S 0. 00222 0. 0227
- 0. 0170 0. 975
- 0. 0308 0. 357
- 0. 0390 0. 356
E 0. 00096 0. 0204 - 0. 0170 0. 975
10000 R 0. 00123 0. 0222 - C. 0308 0. 357
S 0. 00122 0. 0227 - 0. 0399 0. 356
: 1
i t
0. 00126 0. 0231 - 0. 0513
+
t+
Cl assi cal thi n pl ate sol uti on [ Ref . 71.
A refined higher-order Co plate bending element 183
4. T. J. R. Hughes, M. Cohen and M. Haroun. Reduced and
selective integration techniques in the finite element analysis
of plates. Nucl. Engng Design 46.203-222 (1978).
5. T. J. R. Hughes and hf. Cohen, The Heterosis finite ele-
ment for plate bending. Compuf. Structures 9. 44.5-450
(1978).
R. D. Mindlin, Influence of rotary inertia and shear on
flexural motions of isotonic elastic plates. 1. ADD/. Mech. 18.
31-38 (1951). ..
S. P. Timoshenko and S. Woinowsky-Krieger. Theory of
Plates and Shells. McGraw-Hill. New York (1959).
E. Reissner. On transverse bending of plates, including the
effects of transverse shear deformation. Int. .I. Solids Struc-
tures 11,569-573 (1975).
T. Kant, The segmentation method in the analysis of elastic
plates with two opposite edges simply supported using a
refined higher order theory. Rep. C/R/366/80. Department of
Civil Engineering, University of Wales. Swansea SA2 8PP,
1980 (to be published).
IO.
II.
12.
13.
14.
15.
16.
0. C. Zienkiewicz. The Finite Elemenf Method. McGraw-Hill,
U.K. (19771.
Y. K. Cheung
Finite E~emenf
and M. F. Yeo, A Practical Introduction to
Analysis. DD. 94%. Pittnan, London (1979).
A. K. Noor and S. J:Har&y. Evaluation of element stiffness
matrices on CDC STAR-100 computer. Compuf. Structures
9. 151-161 (1978).
A. K. Gupta and B. Mohraz. A method of computing
numerically integrated stiffness matrices. Int. 1. Num. Mefh.
Engng. $83-89 (1972).
E. Hinton and D. R. J. Owen, Finite Element Programming.
Academic Press. London (1977).
0. C. Zienkiewicz. R. L. Taylor and J. M. Too, Reduced
integration techniques in general analysis of plates and
shells. Inf. 1. Num. Meth. Engng. 3, 275-290 (1971).
0. C. Zienkiewicz and E. Hinton. Reduced integration, func-
tion smoothing and non-conformity in finite element analysis
(with special reference to thick plates). L Franklin Inst. 302,
443-461 (1976).
tat Flexure contribufion.
41,-2.i
N. 9
0
x,, I-2,
- N,N
0,11-2.1
.?
2
- N;.dN, 0
z
p,
5 Dell-2.) N, N,
5D,W2r)N,3 +&(I-Zv)N,N,
0
+DD:ll-2v)N. 9 $ D:(I-2v) N,N,
0
$D,(I-2) N, 9
0
+D,(I-2u) N, N, 4 Dfil-2~1 N. N,
fb) Transverse shear contribution.
tK, =_Eh/(l +v)(l - 2s). D, =K,h%2. of =3D,h2/20, D:* =5mh2/28