- Email: [email protected]
Finite element analysis of geometrically nonlinear plate behaviour using a mindlin formulation
Compufm & Sfnrcfums Vol. I I, pp. 20345 PergamonPress Ltd.. 1980. Printed in Great Britain
FINITE ELEMENT ANALYSIS OF GEOMETRICALLY NONLINEAR PLATE BEHAVIOUR USING A MINDLIN FORMULATION A. PICAt Dipartimento di Scienze dell’informazione, Universita’ di Pisa, Italy and R. D.
WODDt
and E. HINTONt
Department of Civil Engineering, University College of Swansea, Swansea, Wales (Received I5 December 1978; received for publication 31 January 1979)
Ah&rue-The paper presents a finite element analysis of the geometrically nonlinear behaviour of plates using a Mindlin formulation with the assumption of small rotations. A comparison of the performance of Linear, Serendipity, Lagrangian and Heterosis elements is given for square, skew, circular and elliptical plates subjected to distributed and point loading. All results are presented numerically. Generally the higher order elements give good results although no particular element emerges as being “best”.
INlRODUCTlON
Linear elastic analysis of plates using the finite element displacement method may be formulated using either Kirchhoff [I] or Mindlin[Z] plate theory. In the former theory shear deformation is ignored, that is, normals to the midsurface remain normal after deformation and in a finite element context, shape functions have to ensure C’ continuity or otherwise a non-conforming formulation must be adopted. Alternatively in Mindlin plate theory, which allows for shear deformation, normals remain straight, but not necessarily normal to the midsurface after deformation. More significantly shape functions only require c” continuity. This latter requirement makes Mindlin plate elements particularly easy to formulate with the result that they have been successfully applied to many linear problems[3,4]. Such formulations readily accomodate thick or thin plates, variable thickness, curved boundaries and laminated materials. The initial confidence expressed in these elements was later subject to doubt on two accounts. Firstly the Serendipity elements were found to exhibit an overstiff “locking” behaviour in very thin plate situations. Secondly Lagrangian elements with (the necessary) selective integration schemes contained spurious zero energy modes not associated with rigid body movement. These latent problems, in an otherwise elegant formulation, led to extensive investigation[5,6] of Serendipity and Langrangian elements, with reduced and selective integration, in the static and vibration analysis context. Generally it emerged[5] that the 9 noded Lagrangian element with selective integration of the shear term was best, although even this occasionally failed, albeit under exceptional circumstances due to the propagation on excitation of a spurious zero energy mode. This remaining problem appears to have been resolved with the emergence of the combined SerendipitylLagrangian “Heterosis” element [7] which eliminates the offending zero energy mode.
Xecturer.
Geometrically nonlinear (GNL) analysis of plates has naturally followed the same development with the majority of applications being based on the Von Karman non-linear plate equations [8- 101.These admit a coupling between bending and membrane behaviour with the retention of the Kirchhoff normality constraint. Recently, however, the Mindlin plate formulation has appeared in the GNL context [ 111. The aim of this paper is to provide a GNL Mindlin formulation under the assumption of small rotations and to present results for Linear, Serendipity, Lagrangian and “Heterosis” elements. Analyses are carried out for clamped and simply supported square plates, skew plates and clamped circular and elliptical plates. Where appropriate, answers are compared to “exact” thin plate solutions. Because of the general absence of numerical results in the literature, no comparison has been possible with other linite element solutions with the exception of the results obtained using the thin plate programs developed by the second authorH21. Consequently this paper takes the opportunity to present all the results numerically to aid future studies. Numerical integration and nodal stress calculation are discussed.
PLATE DEFIXMATIONS
In Mindlin plate theory the displacements u = [u,u,wl’ at (x,y,z) are expressed as functions of the midplane (z = 0) translations il, r? and rt~and independent normal rotations 6, and f3,,see Fig. l(a), as, u = W,y)-Z&(X,Y)
(la)
U= 8 (X,Y)- 2% (X.Y)
(lb)
w = $(x,y)
(lc)
0, and 6,. are the rotations of the normals, to the undeformed midplane, in the xz and yz planes respectively. These normals are not necessarily normal to the midplane after deformation and consequently shear deformation is permitted, see Fig. l(b). 203
204
A. PICA et al. and finally the nonlinear component of inplane strain 1 a+ * 2%( ) ‘_ lak* 4 -* : -- ( ) 1 2 ay a* a+ -.. ax ay
(44
Equations (4) represent the generalised strains. [ai
Mad-plane
I b]
deformattons
Deformotton of In x-z plane
cross-sectIon
STRESSES
Fig. I. Deformations.
The Piola-Kirchhoff stress vector u associated with the Green’s strains Q of eqn (3) is,
STRAINS
For a Mindlin plate the relevant Green’s strain vector
TOTAL LAGRANGIAN VIRTUAL WORK EQUATION
$+;(g+;(g)‘+;($) $+;(&)‘+;(g)*+;(g)* au au au au au au aw aw ~+~~~.~~~.-~~~~ ay ax ax ay ax ay ax ay au aw au -+-+_._+-._+-.a2 ax ax do aw au -+-+-._+-._+-.-~ a2 ay ay
au au a0 a2 ax a2 au av av a.2 ay a2
aw ax aw ay
The basis of the formulation is the virtual work equation for a continuum written in a total Lagrangian coordinate system under the assumption of small strains and (2) conservative loading as [ 141
I
aw a2 aw
a2
Introducing the Von Karman assumptions [ 131 which imply that derivatives of u and D with respect to x, y and z are small and noting that w is independent of z allows Green’s strain to be rewritten in terms of the midplane deformations of eqn (l), as,
de=u dv =
where V is the undeformed volume, u the Piola-Kirchhoff stress vector, dr the virtual Green’s strain vector due to the virtual displacements du, p the mass density, q the body forces per unit mass and p surface tractions acting over an undeformed area a. For the plate formulation this virtual work equation is to be rewritten in terms of area integrals over the midsurface. The internal virtual work d Wi can be written as an integral over the plate area A and thickness t as, dWi =
where the linear midplane strains are
o_
-
11 1 g
ay 1
”
(dzu, + dC,u, + dyxyrxy (7)
A
!!!+aU ay ax
the bending strain
dC=u du =
Substituting the strain expression from eqn (3) into eqn (7) and integrating over the thickness allows d Wi to be rewritten only as an area integral giving,
ax
6
I
+ d yXZrXZ + d yyZryr) dz da.
aa
the inplane strain
(6)
d Wi =
IA
dgT& da
(8)
in which the stress resultant vector 6 is,
ae,
ax
&O = -
ae, ay
ia+& . ay ax
and the shear strain
WI which contains the following components, see Fig. 2, inplane 4,
ai e
6
o=
!I -ax a* --
x
e
ay Y
=
bending
IN,, N,, NxyI’ =
1dz 1T m
(ax, uyy.7xy
Geometrically
205
nonlinear plate behaviour using a Mindlin formulation
q becomes simply, (124
i = [WA WA pqzr,0, OIT
where typically pqXr is the force per unit area in the x direction. Finally the third integral of eqn (6) which represents the external virtual work due to tractions may be split into two terms relating to the plate surfaces (a,), z = ITt/2 and the plate edges (a.), see Fig. 3 as, dW>X=~Sdu’Pdn+~~du’P.do
(13a)
where the tractions on the surfaces z = c f/2 are Fig. 2. Stress resultants.
p = [PX,P,,
WI
Pzl’
shear and the global edge tractions P, are given in terms of the local edge tractions P., P, and P, as,
The generalised strain vector P corresponding to the stresses B of eqn (8) is,
(13c) where c = Cos B and s = Sin 0, see Fig 3. Y4 P"
where the linear component 8’ is, 0
il iI
go=
lP QO
0
(lob)
EJ
I
-----_)
and the nonlinear component fL is,
pL=
4
L
o
Substituting for u, u and w from eqn (1) into eqn (13a)
WC) and integrating the edge area integral explicitly through
0
The second integral in eqn (6) which represents the external virtual work due to body forces may be rewritten as, dW:x =
+duq, +dwqz)dz da
(11)
the thickness gives the work term d W:, as the sum of an area integral over a, and a line integral over an edge length s. dWb=~iT@dct~diT6,dr
(12a) where the midplane displacement vector i is,
and the generalised body force vector q is
iii, I’
(t4b)
in which typically # = 5 (r/2)P,, where - and t refer to forces on the top and bottom surfaces respectively. The generalised edge forces are
= [I
UW
(14a)
where B = [PX,Pw P*, k,
substituting for u, u and w from eqn (1) and integrating explicitly through the thickness allows eqn (11) to be expressed as an area integral giving,
ii = [ff, 8, I+,8x,ey1=
x
Fig. 3. Edgetractions.
II2
(P,‘, P,: P,‘, - Px’z, - Py’z) dz 42
T 1 .
(14c)
The virtual work eqn (6) may now be written entirely in terms of “mid plane” quantities as,
t/2
=
if
[I
_,,*
dqx, %*4.9 - qxz,- qyz)dz I =
(12c)
the body force vector pq is constant over the thickness
CONSlTNTtW
RELATIONS
For a linear elastic isotropic, homogeneous material with the stress in the z-direction o, = 0 the stresses u of
206
A.
PICA
et al. VARIATION OF STRAIN
eqn (5) are related to the Green’s strains E of eqn (3) by, a=Dr
(16a)
where the matrix D of constitutive coefficients is, iY Vl E
0 0
0 0
0 0
OOI-VO 2
D=m
00
Before proceeding with the discretisation of the virtual work eqn (15) it is necessary to consider further the variation of strain di due to the virtual displacements dii. Generally dP is given as the sum of the variation of the linear and nonlinear generalised strains as dB =di”tdBL.
0
(16b)
Since 6’ is a linear function of ii a
!$EO
0
000
ax dP” =
0
0
0
0
0
0
I
a I?( i$ ax !
0
d
0
0
/o -i
0
0
jo
0
0
P/ ax
0
Ja
0
a
0
-ay oio _a a ._...__...._.~___________“_____“x. 0
or
0
-I
1ay
0
-1
If the displacement gradients with respect to the lateral displacements t are,
in which A
/o
._____________I.________________________~
(17a)
Er Dp=j-qD;
0
oyo
oy
where E is Young’s modulus and tr is Poisson’s ratio. Substituting for E from eqn (3) and integrating explicitly through the thickness t enables the stress resultants b of eqn (9) to be written in terms of the generalised strains P of eqn (10) as,
(18)
-
Et’
.
D~=12(1_v2)D;
-
Gt
&=;I2
(17c)
e=
a3 ar;l T [ ax ’ ay 1
(20)
where
rlv
then the variation of the nonlinear component of the inplane strain is obtained from eqn (4d) in terms of the virtual gradients dg as,
07
G=E
2(1 t v) (17d)
depr. = Ae de
(2la)
where
and cx is a shear correction factor to allow for nonuniform shear stress distribution.
(21b)
and dezd
auatiT
[ ax ’ ay
(2lc)
1
FINITE ELEMENT FORMULATlON OF FQUILIBRIUM EQUATIONS The displacement ii within an element is given as a function of n discrete nodal displacements as,
(22a) LAGRANGIAN
IOLI
OS plus hlearchlcc, bubble functlon N,=ll-~‘111-~‘1
HETEROSIS
I OH I
05 OL wth w,. : 0
hlearchlcol
Fig. 4. Mindlin plate elements.
where Nj are shape functions, see Fig. 4, Is, is a 5 x matrix and the nodal displacements 6, are, si = [/ii, Bi,CJi,&i, 8,)I“.
5
unit
(224
Geometrically
nonlinear plate behaviour using a Mindlin formulation
For the Lagrangian and Heterosis elements a hierarchical formulation is adopted. Equation (22a) may be conveniently rewritten as
where
Boi = ii=Nii
207
(22c)
1
B& 0 0 B& [ 0 B&
where
md
in which component matrices are inplane (22d) aNi aNi
1 J ay-x
The virtual displacements di are written in terms of the nodal virtual displacements dg from eqn (22~) as di=Ndg
(2W
bending
where
IIsn St
dg=d
1
B:i =
f
alvi
(23b)
O -ay
-ax
shear The displacement gradients 8 of eqn (20) may now be written in terms of the nodal displacements 6 and cartesian derivatives of the shape functions as
aNi
,.,&=L-Ni ’
(244
8=G6
where
[$ The strain matrix
G=[C,,G*...Gnl and
(244
1
0
(26e)
ahri
-Ni
1
(26f)
’
BL consists of the nodal submatrices BL = [BL,, BLZ . . . B,_,I
Wd
where
oo%oo
(244
(26h)
The virtual gradients of equation (21~) are now written in terms of virtual nodal displacements as
Substituting eqns (23a) and (25) into (19) and (21a) respectively yields the strain variation in terms of the virtual nodal displacements dg as
G, = [
(-)()aNiO()’ ay
de=Gdg
The generaked Green’s strain vector P of eqn (lOa) is given in terms of nodal displacements 6, displacement gradients Ab and Cartesian derivatives of N as,
d = &+;B~(b)]g
(2W
where & is the strain matrix giving the linear straini lPorboand r,” and BL, which is linearly dependent upon 6, gives the non-linear strains 4‘. Consequently the nonlinear strains are quadratically dependent upon the nodal displacements 6. The constant matrix BO is written in terms of nodal submatrices as,
BO=PO,,Bm...Bonl
dP=Bdg
(25)
(26b)
(27a)
in which B = Bot BL(S).
(27b)
The virtual work eqn (15) is discretised, for an element, by substituting eqns (23a) and (27a) for dQ and di respectively giving (2ga) where the equivalent nodal load vector R due to body forces and tractions is NT@da + 1. NT fie ds. R=j-ANT4da+j-s
(28b)
Since the nodal virtual displacements d& are arbitrary
208
A. PICAet al.
the element nonlinear equilibrium equations become,
& = [o*, l72 (29)
e(6)= ( B*&da-R=O.
. .
. &IT.
(W
This enables the equilibrium eqns (29) to be expanded as
The load vector R may also contain nodal point loads. Equation (29) can have the dual role of representing either the element, or in an assembled form the total equihbrium equation. It is a nonlinear equation in 6 since B and C?are linear and quadratic functions of 8 respectively.
da-R=0
(35)
SOLUTION To N0NLlNJ.M EQuILlaRNM EQUATIONS The solution algorithms for the assembled nonlinear equilibrium eqns (29) are based on the Newton-Raphson method which consists of a series of linear solutions. If an initial estimate 6, for the total displacements gives residual (unbalanced) forces #(&) # 0 than an improved value W+r is obtained by equating to zero the linearised Taylor’s series expansion of S(L+,) in the neighbourhood of 6, as
#(&+I) z $(SI) + KTAL
= 0
(30)
Differentiating & with respect to 8, yields a typical term in KT as,
where (aRJa8,) = 0 for conservative loading. It is shown in the appendix that the tangent stiffness matrix can be written as,
where Kr is known as the assembled tangent stiffness matrix evaluated at Si and given by (31) Equation (30) is the linear incremental equilibrium equation which gives the linearised approximation to the relation between the residual forces and incremental displacements, A& at a point Si on the equilibrium path. The improved solution is then found as, St+, = 6, t A&.
Kr=/~BTBBdo+~_-GT[;Y
:]Gda.
(37)
Substituting for B, b and C from eqns (27b), (17b) and (24b) respectively gives, KT=&tKL+Krr
(38d
where the constant linear elastic stiffness matrix Ko is, &=/-&TB&dn
(38b)
(32)
Equations (30) and (32) represent the Newton-Raphson solution to the nonlinear eqns (29). The terms $($) and ASi give measures of the convergence of the solution To improve numerical stability and to give intermediate results, the load R is usually applied in increments. Since the Newton-Raphson method (N - R) requires repeated calculation and inversion of Kn a modified N-R method (KT2) is also employed whereby KT is calculated once only on the second iteration of each load increment.
the initial displacement matrix KL, which is quadratically dependent upon displacements 8, is KL =
IA
B,,‘BBL da t
IA
BL’&, da t
IA
BL=&B,_ da (38~)
and the initial stress matrix Id, is
K,=jjT[zy 2]Gde.
(384
TANGEhT SllFFNEss MATRIX The tangent stiffness matrix given in eqn (31) may be written as,
KT =
(33)
For computational efficiency it is necessary to consider the various nodal substitkess matrices contained in eqn (38a), this avoids much zero multiplication. Using eqn (26c) for the nodal submatrix Boi allows K,,, to be expressed in terms of inplane, bending and shear components as Koij =
where m is the number of degrees of freedom per element. For convenience the eight stress resultants of eqn (9) are rewritten as,
(39) where ij = 1 to n. From eqn (26c) and (26h) the nodal initial displacement
Geometricallynonlinearplate behaviourusing a Mindlinformulation
209
SI'RESS CALCULATION
matrix KLy is,
For all the Mindlin elements stresses are extrapolated from Gauss point values in a manner consistent with the integration rule used for the tangent stiffness matrix calculation. The extrapolation which results is given in Table 2. in which ki = AeG,
Examples (i) Preliminary. In the examples deflections, stresses and loads are given in non-dimensional form as follows,
Wb)
where the matrix Gib is,
aNi Gib =
_$ i
[ ay
w-w
0 0’ 1
The predominance of zeros in the submatrix Gi of eqn (24c) indicates that the nodal initial stress matrix Iij is best calculated explicitly giving the only non zero term at the 3 x 3 position in the 5 x 5 matrix as,
w = $(0,0)/t Q = (qa’)/(Et’) P = (pa*)/(Et’)
Central deflection Load (a) uniformly distributed (b) point Extreme fibre stress
flit-&Y) = u(x,y)a2/(Et2)
where a = characteristic plate dimension; q = uniformly distributed load; p =point load; i = x or y; U(X,y)= maximum extreme fibre stress at x,y. All examples use the KT2 modified Newton-Raphson approach and convergence is checked using a total residual norm criteria as,
(41)
EQUIVALENTNODALLOADSDUETOSTlWSSEs
To evaluate the residual nodal force vector +4(S) of eqn (29) the equivalent nodal forces due to the stress resultants 6 may be written for a typical node i as, Pi =
(RpOi)T&p da (k,,,’ kP + (B&)& •t (B&)T& I -CAL
(42)
INTEGRATION
The integrals in eqns (29) and (39-41) are evaluated using numerical integration. The orders of the Gauss quadratures used, see Table 1, are selective and depend upon whether the iteAmunder considera$on is a function of the membrane (D,, &,,). bending (I&,,&b) or shear (D,, &,) behaviour. The equivalent nodal load vector R due to the external forces is calculated using a 3 x 3 ru1e.t
where y = convergence tolerance (%). Figures in brackets in Tables 3-18 are percentage errors. (ii) Clamped square plate. This example has become a “standard” used by many authors[8,9,12] to check GNL thin plate formulations. The analytical thin plate solution is given by Levy [IS] who solved Von Karman’s equations using a double Fourier series; this solution is quoted as having a possible error of less than 2%. Since the problem is symmetric only a quarter mesh need be used, see Fig. 5. Tables 3-5 show results for central deflection, central and edge stress respectively. These provide a comparison between the four Mindlin elements and two thin plate elements obtained from a previous study[l2]. The thin plate elements are the Irons-Razzaque non-conforming triangular element[l6] (RI) with linear inplane behaviour and the slope conforming rectangular element introduced by Bogner et a/.[171 with parabolic inplane behaviour (BFS).
Table 1. Integrationrules LN
QS
QL
QH
MEMBRANE
2x2
2x2
3x3
3x3
BENDING
2x2
2x2
3x3
3x3
SHEAR
1x1
2x2
2x2
2x2
Table 2. Stress extrapolation LN
QL
OH
MEMBRANE
BILINEAR
BILINEAR
BIQUADRATIC
BIQIJADRATIC
BENDING
BILINEAR BILINEAR
BIQUADRATIC
BIQUADRATIC
SHEAR
CONSTANT
t2 X2 for Linear element. CAS Vol. II.No. SE
QS
BILINEAR BILINEAR
BILINEAR
210
Geometrically nonlinear plate behaviour using a Mindlin formulation Table 3. Clamped square plate-central deflections W
l-
FINITE ELEMENT
qs
LN
QH
qL
BFS
RI
0.2349
(0.87)
0.2351
(0.78)
0.2363
(0.08)
0.2361
(0.38)
0.2387
(0.73)
0.2361
(0.38)
0.4705
(0.10)
0.4673
(0.79)
0.4699
(0.23)
0.4687
(0.49)
0.4717
(0.16)
0.4685
(0.52)
0.6990
(0.57)
0.6887
(0.90)
0.6915
(0.50)
0.6902
(0.69)
0.6916
(0.48)
0.6900
(0.73)
(1.00)
0.9015
(1.15)
0.9008
(1.23)
0.9012
(1.18)
95.
0.912
0.9191
(0.78)
0.9003
(1.29)
0.9029
134.9
1.121
1.1318
(0.96)
1.1041
(1.51)
1.1063
(1.31)
1.1050
(1.43)
1.1025
(1.65)
1.1047
(1.46)
1.3350
(0.91)
1.2990
(1.81)
1.3009
(1.67)
1.2997
(1.76)
1.2961
(2.03)
1.2992
(1.80)
1.5347
(0.90)
1.4913
(1.95)
1.4928
(1.85)
1.4916
(1.93)
1.4879
(2.18)
1.4908
(2.00)
1.7439
(1.74)
i.6774
(2.13)
1.6786
(2.06)
1.6775
(2.13)
1.6744
(2.31)
1.6761
(2.21)
1.9109
(0.47)
1.8682
(1.77)
1.8555
(2.44)
1.8545
(2.50)
1.8529
(2.58)
1.8524
(2.61)
Table 4. Clamped square plate--centre stresses ux(O,O)
ri
LOAO Q
I?t
FINITE ELEMENT
EXACT (15)
BFS
RI
QH 2.6144
(0.55)
2.6890
(3.42)
2.6328
5.4520
(4.85)
5.4140
(4.11)
5.4047
(1.26) (5.47)
8.2908
(3.63)
8.0205
(0.26)
8.3372
(4.21)
17.7!
2.6
38.3
5.2
63.4
8.0
95.
11.1
11.575
(4.28)
11.115
(0.14)
11.103
(0.03)
11.066
(0.30)
10.521
(5.22)
11.127
(0.24)
34.9
I1 I1
13.3
14.567
(9.53)
13.817
(3.88) / 13.827
(3.96)
13.789
(3.67)
12.971
(2.47)
13.866
(4.26)
84.
15.9
17.493
(10.02)
16.461
(3.53)
16.497
(3.75)
16.456
(3.49)
15.390
(3.21)
16.552
(4.10)
i9.2
20.446
(6.49)
19.160
(0.21)
19.225
(0.13)
19.178
(0.12)
17.885
(6.85)
19.296
(0.50)
18.
21.9
23.605.
(7.78)
21.902
(0.00)
21.994
(0.43)
21.938
(0.17)
20.438
(6.67)
22.081
(0.83)
-02.
-25.1
26.318
(4.85)
24.805
(1.18)
24.780
(1.27)
24.713
(1.54)_
23.020
(8.29)
24.881
(0.87)
12
45.
I!
n
8.5570
(6.96)
8.3528
(4.41)
8.3258
(4.07)
Table 5. Clamped square plate-edge stresses u*((a/2,0)
r
LOAD Q --
1 ‘i 17.7:
38.3
LN
QS
5.48
3.3663
(38.6)
il.52
7.1122
(38.3)
5.3256
QL
f
(2.82)
11.155
(3.17)
5.3163
W
(2.99)
5.2991
RI5.2523
(3.30)
11.216
(2.64)
11.189
(2.87)
(4.16)
10.807
(6.19)
5.2915
(3.44)
11.137
(3.32)
63.4
18.03
11.198
(37.9)
17.515
(2.86)
17.726
(1.68)
17.697
(1.64)
16.680
(7.49)
17.556
(2.63)
95.
25.32
15.639
(38.2)
24.522
(3.15)
24.967
(1.39)
24.946
(1.48)
23.001
(9.16)
24.649
(2.65)
134.9
33.5
20.461
(38.9)
32.278
(3.65)
33.045
(1.36)
33.038
(1.38)
29.875
(10.8)
32.485
(3.03)
84.
42.4
25.597
(39.6)
40.710
(3.98)
41.885
(1.22)
41.900
(1.18)
37.257
(12.1)
40.961
(3.39)
245.
52.8
31.172
(41.0)
50.043
(5.22)
51.719
(2.05)
51.763
(1.96)
45.357
(14.1)
50.251
(4.83)
1318.
63.9
37.459
(41.4)
60.064
(6.00)
62.325
(2.46)
62.401
(2.34)
54.011
(15.5)
60.097
(5.95)
75.8
43.121
(43.11
70.897
(6.471
73.407
(3.16)
73.520
(3.01)
63.00
(16.9)
70.190
(7.40)
1402.
1
!
Y t r-------
,’
_,
/j
,’
,‘,
’
/’
’
’
,I
,’
x
I
I
I
I
I
I
L____
-_
t___2____-_
.
L
From the tables it can be concluded that the use of Mindlin formulation in the GNL context is an acceptable
@a_ /’
l
___-
_-__J
o= 300u-l 3,n E= 03 x10' Ibh' Y= 0316
t=
5. Square clamped and simply supported plate under uniformly distributed load-geometry and material properties.
Fii.
1
FINITE ELDIENT
EXACT (15)
T
alternative to the thin plate Kirchhoff formulation. Generally the 8 and 9 noded elements (QS, QL, QH, BFS) give comparable results for deflections and stresses whilst the elements having linear inplane behaviour (LN, RI) yield inferior stresses. A surprising result is the accuracy of the linear (LN) element for displacements. These results were obtained using a convergence tolerance y = 1%. (iii) Simply suppotied squareplate. Here the Mindlin formulation is compared to the thin plate Kirchhoff results obtained by Rushton[l8] who used a finite difference dynamic relaxation approach with a quoted error of less than 0.5%. The boundary conditions used are as follows, on x=ka/2;
u=v=w=tl,=O 0, free
(W
Geometrically nonlinear plate behaviourusing a Mindlin formulation
on y=+.a/2;
211
P
u=u=w=&=O 0, free.
Wb)
Plate geometry, material properties and mesh idealisation are shown in Fig. 5. The results for centre deflections and extreme fibre stresses are shown in Tables 6 and 7 respectively. A convergence tolerance of y = 1% was used. Again the 8 and 9 noded elements perform well in comparison with the finite difference results whilst the linear element gives less accurate, but acceptable, values. (iv) Clamped skew plates. Here the ability of the elements to mode distorted plates is explored. Firstly by reanalysing the clamped square plate, using the mesh shown in Fig. 6 (B = 0). This eleven Heterosis element mesh, for the quarter plate, gives displacements comparable with values obtained using the previous 16 element mesh. Inevitably, though, this coarse distorted mesh produces a degradation in stresses. Further distortion of the mesh allows two clamped skew plates @I= 30”, 45”)
Fig. 6. Clamped skew plate under uniformly distributed loadmesh details (see Fig. 5 for a, 1, E and v).
to be analysed. These latter problems have been discussed previously in Ref. [ 101where triangular Kirchhoff elements were employed, but no numerical results given. See Tables 8-10 for displacement, centre and edge stress results respectively. In these examples a convergence tolerance of y = 0.1% was used. (v) Clamped circular plate. Here the performance of
Table 6. Simply supported square plate--centre deflections W REF.
FINITE ELEMENT
1
(18)
LN I
1
I
QS
QL
I
I
QH I
0.335
0.2531 (5.40) 0.3478 (3.82)
0.3480 (3.84) 0.3478 (3.82)
0.818
0.8303 (1.51) 0.8184 (0.05)
0.8185 (0.06: 0.8184 (0.05)
1.47
1.4931 (1.57) 1.4655 (0.31)
1.4657 (0.29)
2.40
2.4380 (1.58) 2.3927 (0.30)
2.3932 (0.28) 2.3928 (0.30)
3.83
3.8985 (1.79) 3.8124 (0.46)
3.8134 (0.43) 3.8128 (0.45)
6.07
6.1571 (1.43) 6.0521 (0.29)
6.0539 (0.26) 6.0530 (0.28)
I
1.4655 (0.31)
Table 7. Simply supported square plate-centre stresses gx(0,O) FINITE ELEMENT
REF.
LOAD
Q
(18)' LN
9.16 36.6 146.5
QS
QL
QH
2.46
2.6354 (7.13) 2.6214 (6.56)
2.6029 (5.81) 2.6016 (5.76)
6.90
7.1837 (4.11) 7.0026 (1.49)
6.9826 (1.20) 6.9787 (1.14)
14.5
15.263
(5.26) 14.644
(1.00) 14.635
(0.93) 14.615
(0.79)
586.
30.0
31.396
(4.65) 30.183
(0.61) 30.188
(0.63) 30.130
(0.43)
2344.
65.2
68.428
(4.95) 65.673
(0.73) 65.756
(0.85) 65.661
(0.71)
9377.
148.3
(1.14) 149.87
(1.06)
155.24
(4.68) 149.66 (0.92)
Table 8. Clamped skew plates-centre
displacements W
FINITE ELEMENT
LOAD 17.79
149.99
6 = o" * 0.2291 (3.35)
B -5300
6 = 450
0.1412
0.0690
38.3
0.4576 (2.84)
0.2940
0.1473
63.4
0.6781 (2.43)
0.4585
0.2394
95.
0.8900 (2.41)
0.6309
0.3476
134.9
1.0947 (2.35)
0.8061
0.4703
184.
1.2906 (2.45)
0.9779
0.6013
245.
1.4838 (2.44)
1.1487
0.7388
318.
1.6708 (2.52)
1.3140
0.8762
402.
1.8488 (2.79)
1.4708
1.0083
*
error with respect to Ref. (15). see table 3.
I
212
A.
PICA et ai.
Table 9. Clamped skew elate-extreme fibre stresses at centre in direction normal to the edge FINITE ELEMENT
LOAD 6 =
o”
6 = 450
17.79
2.8929 (11.3)
1.5930 3.4775
38.3
6.0767 (16.9)
63.4
9.3357 (16.7)
5.7852
95.
12.590 (13.4)
a.5907
134.9
15.840 (19.1)
11.863
184.
19.074 (19.9)
15.425
245.
22.418 (16.7)
19.209
318.
25.845 (18.0)
23.016
402.
29.315 (16.8)
26.708
Table 10. Clamped skew plate-midside extreme fibre stresses in direction normal to the edge LOAD Q
FINITE ELEMENT 6 = 00
!
(12.2) (11.5)
B = 450
6 = 300 3.0991
1.6819
6.6972
3.6455 6.0640
(10.3)
10.974
(10.0)
16.038
(10.2)
21.925
12.875
(10.4)
28.524
17.386
(11.7)
35.959
22.761
(12.7)
44.020
(13.9)
52.452
9.1004
.....-. . Fish 28.868 35.504
’
.
-
b:2OOtn t =21n E :O 1x10@ lb/r? " =0.3
Fig. 7. Clamped circular and elliptical plates-geometry and material properties. the three “Q-family” elements (QS, QL, QH) is com-
pared with analytical Kiichhoff solutions for uniformly distributed[l9] and central point[20] loading. For both problems the curved boundary mesh shown in Fig. 7 is used. The results, using a convergence tolerance of O.l%, for displacements and stresses are given in Tables 11-15. Generally it is dithcult to distinguish a clear trend in the behaviour of the elements; whilst displacements compare well with the analytical solutions, the stresses are irregular, an element giving good accuracy for one load case and bad for the other. Nevertheless the Heterosis (QH) element consistently overestimates the edge stresses for both loading conditions. (vi) Clamped elliptical plate. Finally results are presented for the clamped elliptical plate under uni-
formly distributed load, shown in Fig. 7, these being compared with the analytical solution given in Ref. [19]. These results are shown in Tables 16-18 and were obtained using a 0.1% convergence tolerance. Centre displacements are in good agreement with the analytical results for all elements. For the QS and QL elements, the stresses are acceptable, the centre stresses being slightly more accurate. Again an irregularity is observed-the QH elements providing better stresses at the edge than at the centre. CONCLUSIONS
This paper has presented a geometrically nonlinear analysis of plates using a finite element Mindlin formulation. The paper contrasts the performance of the
Geometrically nonlinear plate behaviour using a Mindlin formulation
213
Table 11. Clamped circular plate (UDL)-centre displacements W LOAD
FINITE ELEMENT
REF.(19)
I Q
QS
QH
QL
1.
0.169
0.1614
(4.50)
0.1682
(0.47)
0.1644
(2.72)
2.
0.323
0.3116
(3.53)
0.3231
(0.03)
0.3167
(1.95)
3.
0.457
0.4453
(2.56)
0.4591
(0.45)
0.4514
(1.23)
6.
0.761
0.7566
(0.57)
0.7702
(1.21)
0.7626
(0.21)
10.
1.035
1.0417
(0.65)
1.0514
(1.58)
1.0459
(1.05)
15.
1.279
1.2950
(1.25)
1.3007
(1.70)
1.2974
(1.44)
Table 12. Clamped circular plate (UDL)-centre stress o,(O,O) r
,r REF. (19)
FINITE ELENENT
QL
QS
QH
1.
0.502
0.5018
(0.03)
0.5304
(5.66)
0.5119
2.
0.986
1.0044
(1.87)
1.0611
(7.61)
1.0270
(1.97) (4.15)
3.
1.421
1.4724
(3.62)
1.5524
(9.25)
1.5079
(6.17) (8.29)
6.
2.477
2.6070
(5.25)
2.7384
(10.5)
2.6824
10.
3.542
3.6774
(3.82)
3.8657
(9.14)
3.8077
(7.50)
15.
4.622
4.6595
(0.81)
4.9163
(6.38)
4.8567
(5.08)
Table 13. Clamped circular plate (UDLj-edge stresses u,(O,a) LOAD
Q
FINITE ELEMENT
REF.(19)
QL
QS
QH
1.
0.772
0.7831
(1.44)
0.7541
(2.32)
0.9352
(21.1)
2.
1.550
1.5550. (0.32)
1.5046
(2.93)
1.8553
(19.,7)
3.
2.302
2.2924
(0.42)
2.2265
(3.28)
2.7294
(18.6)
6.
4.324
4.2639
(1.39)
4.1799
(3.33)
5.0393
(16.5)
10.
6.517
6.4692
(0.73)
6.4094
(1.65)
7.5871
15.
8.694
8.8125
(1.36)
8.8271
(1.53)
10.267
(16.4) (18.1)
Table 14. Clamped circular plate (point load)-centre displacements W REF.(20)
T-
I ELEMENT FINIITI QL
QS
0.2130
0.2044
(4.04)
0.2149
(0.89)
0.2085
(2.11)
0.4052
0.3908
(3.56)
0.4088
(0.88)
0.3977
(1.86)
0.5705
0.5528
(3.10)
0.5754
(0.87)
0.5612
(1.62)
0.7123
0.6930
(2.71)
0.7184
(0.86)
0.7021
(1.44)
0.8354
0.8152
(2.41)
0.8427
(0.87)
0.8246
(1.29)
0.9442
0.9237
(2.17) 1
0.9527
(0.90)
I
Table 15. Clamped circular plate (point load)-edge stresses q(0.a) kINI_
LOAD REF. (20)
QH
QL
QS
Q 1.
0.4858
0.5318
(9.48)
0.4932
(1.52)
0.6767
(39.3)
2.
0.9592
1.0478
(9.24)
0.9765
(1.80)
1.3268
(38.3)
3.
1.3974
1.5240
(9.06)
1.4248
(1.96)
1.9194
(37.3)
4.
1.7988
1.9583
(8.87)
1.8355
(2.04)
2.4537
(36.4)
5.
2.1679
2.3553
(8.64)
2.2129
(2.07)
2.9377
(35.5)
6.
2.5110
2.7221
(8.41)
2.5635
(2.09)
3.3816
(34.7)
Linear, Serendipity, Lagrangian and Heterosis isoparametric elements by comparing results with those obtained from analytical and some other numerical thin plate solutions. The validity of the approach has been demonstrated for thin plates with rectangular, skew or curved boundaries. It was hoped that a discernable pattern of element behaviour would emerge from this investigation
but this proved not to be the case. For rectangular element meshes it would seem advisible to use the QH, Heterosis, element in view of its “safe” characteristics referred to in the introduction. However the presence of curved boundaries in the meshes led to inconsistent behaviour for all elements particularly with respect to stresses. Except in the case of the QH element this does not imply inaccurate results.
214
A. PICAet al. Table 16. Clamped elliptical plate--centre displacements W FINITE ELEMENT
REF.(19)
QL
QS 0.3449
(0.25)
0.3462
(0.64)
0.3449
(0.25)
0.600
0.6055
(0.92)
0.6078
(1.29)
0.6063
(1.04)
0.789
0.7985
(1.21)
0.8003
(1.43)
0.7990
(1.27)
1.163
1.1841
(1.82)
1.1845
(1.85)
1.1840
(1.81)
1.490
1.5085
(1.24)
1.5078
(1.19)
1.5080
(1.21)
1.790
1.7968
(0.38)
1.7952
(0.29)
1.7959
(0.33)
Table 17. Clamped elliptical plate-centre LOAD Q
QH
0.344
stresses +(O,O)
FINITE ELEMENT
REF. (19) _
QS
QL
QH
1.
0.862
0.9082
(5.36)
0.9021
(4.66)
0.9128
(5.89)
2.
1.550
1.6652
(7.43)
1.6575
(6.93)
1.6795
(8.36)
3.
2.077
2.2482
(8.24)
2.2384
(7.77)
2.2707
(9.32)
6.
3.217
3.4583
(7.50)
3.4450
(7.09)
3.4984
(8.75)
10.
4.333
4.5462
(4.92)
4.5328
(4.61)
4.6009
(6.18)
L 15.
5.490
5.5945
(1.90)
5.5873
(1.77)
5.6633
(3.16)
Table 18. Clamped elliptical plate-edge stresses ~~(0, a) FINITE ELEIlENT
REF. (19)
QS
The formulation presented herein has been extended to include shallow shells[21] allowing for the investigation of the postbuckling behaviour of plates with initial imperfections. This will be the subject of a further report. REFERENCES 1.
0. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, New York (1978). 2. E. Hmton, A. Razzaque, 0. C. Zienkiewicz and J. D. Davies, A simple finite element solution for plates of homogeneous, sandwich and cellular construction. Proc. Inst. Ciu. Engrs. Part. 2,59,43-65 (Mar. 1975). 3. T. J. R. Hughes, M. Cohen and M. Haroun, Reduced and selective integration techniques in the finite element analysis of plates. Nuclear Engng Design 46.203-222 (1978). 4. T. J. R. Hughes, R. L. Taylor and W. Kanoknukulchai, A simple and efficient finite element for plate bending. Int. J. N&I. Meth. Engng ll(lO), 1529-43(1977). 5. E. D. L. Puuh, E. Hinton and 0. C. Zienkiewicz. A studv of quadrilateral plate bending elements with ieduced . integration. Int. J. Num. Meth. Engng 12, 1059-79(1978). 6. E. Hinton and N. Bicanic, A comparison of Lagrangian and Serendipity Mindlin plate elements for free vibration. comput.Smlctulvs10(3), 483-493 (1979). 7. T. I. R. Hughes and M. Cohen, The “Heterosis” finite element for plate bending. Comput. Structures 9(5), 445-450 (1978). 8. D. W. Murray and E. L. Wilson, Finite element large deflection analysis of plates. A.S.C.E. %(EMl). 143-165 (Feb. 1%9). _ 9. P. G. Bergan and R. W. Clough, Large deflection analysis of plates and shallow shells using the finite element method. Int. 1. Num. Meth. Engng $543-556 (1973).
QL
!
QH
10. R. S. Srinivasan and W. Bobby, Nonlinear analysis of skew plates using the finite element method. Cornput. Structure 6, 199-202 (1976). 11. A. K. Noor, M. D. Mathers and M. S. Anderson, Exploiting symmetries for efficient postbuckling analysis of composite plates. J. AIAA l&24-32 (Jan. 1976). 12. R. D. Wood, The Application of Finite Element Methods lo Geometrically Nonlinear structural Analysis. University of Wales, Swansea, U.K. C/Ph/20/73 (1973): 13. Y. C. Funa. Foundations of Solid Mechanics. Prentice Hall (l%S). -’ 14. R. D. Wood and 0. C. Zienkiewicz, Geometrically nonlinear finite element analysis of beams, frames, &hes and axisymmetric shells. Comout. Structures 7.725-35 (19771. 15. S. ievy, Square plate with clamped edges unde; normal pressure producing large deflections. NACA, Tech. Note 847 (1942). 16. A. Razzaque, Program for triangular bending element with derivative smoothing. Znt. J. Num. Meth. Engng 5, 588-9 (1973). 17 F. K. Bogner, R. L. Fox and L. A. Schmit, The generation of inter-element compatible sti&ess and mass matrices by the use of interpolation formulae. Proc. Conf. Matrix Melh. Srruct. Anal. Wright-Patterson A.F.B. (1965). 18. K. R. Rushton, Large deflexion of plates with initial curvature. Int. J. Mech. Sci. 12, 1037-1051(1970). 19. N. A. Weil and N. M. Newmark, Large deflections of elliptical plates. J. Appl. Mech. 23, 21-26 (1956). 20. R. Schmit, Large deflections of a clamped circular plate. A.S.C.E. 94(EM6), 1603-6 (Dec. 1968). 21. A. Pica, Geometrically Nonlinear Analysis of Mindlin Plates by Finite Hement Method. University of Wales, Swansea, U.K. C/M/132/78(Mar. 1978).
Geometrically
Substituting eqns (27b) and (A6) into (Al) gives
APPENLNX
Substituting eqn (36) into eqn (33) enables Kr to be written completely as, Kr =
I
[B’L(&.) A
215
nonlinear plate behaviour using a Mindlin formulation
do
t IbI(
r,=&ITbd~+f~M(A[~+~‘]dn.
(*7)
(Al) Considering the last integral in eqn (A7) and noting that Bo is not a function of 8 gives, upon expansion,
where
c M(&)B da =
and M(e)
= (A3)
‘a Ulas,
‘a was,
The differentials with respect to 8, in eqn (A8), are the coefficients of the G matrix of eqn (25) typically, for degree of
a “”
““as,
freedom i, Since the stress resultants d are a linear function of the strains i. see eqn (17a), the variation of stress may be written as, dd=ljdG.
(A9) (A4)
Substituting for d0 from eqn (27a) gives dd in the terms of the variation of 8 as, d& = SB d8.
(AS)
It is evident from eqns (A2) and (AS) that L(B) = fiB.
(A6)
and
at?!!.
(ab >
ay=G
21.
(AlO)
Consequently eqn (A7) may be rearranged to give the expression for Kr shown in eqn (37).
FINITE ELEMENT ANALYSIS OF GEOMETRICALLY NONLINEAR PLATE BEHAVIOUR USING A MINDLIN FORMULATION A. PICAt Dipartimento di Scienze dell’informazione, Universita’ di Pisa, Italy and R. D.
WODDt
and E. HINTONt
Department of Civil Engineering, University College of Swansea, Swansea, Wales (Received I5 December 1978; received for publication 31 January 1979)
Ah&rue-The paper presents a finite element analysis of the geometrically nonlinear behaviour of plates using a Mindlin formulation with the assumption of small rotations. A comparison of the performance of Linear, Serendipity, Lagrangian and Heterosis elements is given for square, skew, circular and elliptical plates subjected to distributed and point loading. All results are presented numerically. Generally the higher order elements give good results although no particular element emerges as being “best”.
INlRODUCTlON
Linear elastic analysis of plates using the finite element displacement method may be formulated using either Kirchhoff [I] or Mindlin[Z] plate theory. In the former theory shear deformation is ignored, that is, normals to the midsurface remain normal after deformation and in a finite element context, shape functions have to ensure C’ continuity or otherwise a non-conforming formulation must be adopted. Alternatively in Mindlin plate theory, which allows for shear deformation, normals remain straight, but not necessarily normal to the midsurface after deformation. More significantly shape functions only require c” continuity. This latter requirement makes Mindlin plate elements particularly easy to formulate with the result that they have been successfully applied to many linear problems[3,4]. Such formulations readily accomodate thick or thin plates, variable thickness, curved boundaries and laminated materials. The initial confidence expressed in these elements was later subject to doubt on two accounts. Firstly the Serendipity elements were found to exhibit an overstiff “locking” behaviour in very thin plate situations. Secondly Lagrangian elements with (the necessary) selective integration schemes contained spurious zero energy modes not associated with rigid body movement. These latent problems, in an otherwise elegant formulation, led to extensive investigation[5,6] of Serendipity and Langrangian elements, with reduced and selective integration, in the static and vibration analysis context. Generally it emerged[5] that the 9 noded Lagrangian element with selective integration of the shear term was best, although even this occasionally failed, albeit under exceptional circumstances due to the propagation on excitation of a spurious zero energy mode. This remaining problem appears to have been resolved with the emergence of the combined SerendipitylLagrangian “Heterosis” element [7] which eliminates the offending zero energy mode.
Xecturer.
Geometrically nonlinear (GNL) analysis of plates has naturally followed the same development with the majority of applications being based on the Von Karman non-linear plate equations [8- 101.These admit a coupling between bending and membrane behaviour with the retention of the Kirchhoff normality constraint. Recently, however, the Mindlin plate formulation has appeared in the GNL context [ 111. The aim of this paper is to provide a GNL Mindlin formulation under the assumption of small rotations and to present results for Linear, Serendipity, Lagrangian and “Heterosis” elements. Analyses are carried out for clamped and simply supported square plates, skew plates and clamped circular and elliptical plates. Where appropriate, answers are compared to “exact” thin plate solutions. Because of the general absence of numerical results in the literature, no comparison has been possible with other linite element solutions with the exception of the results obtained using the thin plate programs developed by the second authorH21. Consequently this paper takes the opportunity to present all the results numerically to aid future studies. Numerical integration and nodal stress calculation are discussed.
PLATE DEFIXMATIONS
In Mindlin plate theory the displacements u = [u,u,wl’ at (x,y,z) are expressed as functions of the midplane (z = 0) translations il, r? and rt~and independent normal rotations 6, and f3,,see Fig. l(a), as, u = W,y)-Z&(X,Y)
(la)
U= 8 (X,Y)- 2% (X.Y)
(lb)
w = $(x,y)
(lc)
0, and 6,. are the rotations of the normals, to the undeformed midplane, in the xz and yz planes respectively. These normals are not necessarily normal to the midplane after deformation and consequently shear deformation is permitted, see Fig. l(b). 203
204
A. PICA et al. and finally the nonlinear component of inplane strain 1 a+ * 2%( ) ‘_ lak* 4 -* : -- ( ) 1 2 ay a* a+ -.. ax ay
(44
Equations (4) represent the generalised strains. [ai
Mad-plane
I b]
deformattons
Deformotton of In x-z plane
cross-sectIon
STRESSES
Fig. I. Deformations.
The Piola-Kirchhoff stress vector u associated with the Green’s strains Q of eqn (3) is,
STRAINS
For a Mindlin plate the relevant Green’s strain vector
TOTAL LAGRANGIAN VIRTUAL WORK EQUATION
$+;(g+;(g)‘+;($) $+;(&)‘+;(g)*+;(g)* au au au au au au aw aw ~+~~~.~~~.-~~~~ ay ax ax ay ax ay ax ay au aw au -+-+_._+-._+-.a2 ax ax do aw au -+-+-._+-._+-.-~ a2 ay ay
au au a0 a2 ax a2 au av av a.2 ay a2
aw ax aw ay
The basis of the formulation is the virtual work equation for a continuum written in a total Lagrangian coordinate system under the assumption of small strains and (2) conservative loading as [ 141
I
aw a2 aw
a2
Introducing the Von Karman assumptions [ 131 which imply that derivatives of u and D with respect to x, y and z are small and noting that w is independent of z allows Green’s strain to be rewritten in terms of the midplane deformations of eqn (l), as,
de=u dv =
where V is the undeformed volume, u the Piola-Kirchhoff stress vector, dr the virtual Green’s strain vector due to the virtual displacements du, p the mass density, q the body forces per unit mass and p surface tractions acting over an undeformed area a. For the plate formulation this virtual work equation is to be rewritten in terms of area integrals over the midsurface. The internal virtual work d Wi can be written as an integral over the plate area A and thickness t as, dWi =
where the linear midplane strains are
o_
-
11 1 g
ay 1
”
(dzu, + dC,u, + dyxyrxy (7)
A
!!!+aU ay ax
the bending strain
dC=u du =
Substituting the strain expression from eqn (3) into eqn (7) and integrating over the thickness allows d Wi to be rewritten only as an area integral giving,
ax
6
I
+ d yXZrXZ + d yyZryr) dz da.
aa
the inplane strain
(6)
d Wi =
IA
dgT& da
(8)
in which the stress resultant vector 6 is,
ae,
ax
&O = -
ae, ay
ia+& . ay ax
and the shear strain
WI which contains the following components, see Fig. 2, inplane 4,
ai e
6
o=
!I -ax a* --
x
e
ay Y
=
bending
IN,, N,, NxyI’ =
1dz 1T m
(ax, uyy.7xy
Geometrically
205
nonlinear plate behaviour using a Mindlin formulation
q becomes simply, (124
i = [WA WA pqzr,0, OIT
where typically pqXr is the force per unit area in the x direction. Finally the third integral of eqn (6) which represents the external virtual work due to tractions may be split into two terms relating to the plate surfaces (a,), z = ITt/2 and the plate edges (a.), see Fig. 3 as, dW>X=~Sdu’Pdn+~~du’P.do
(13a)
where the tractions on the surfaces z = c f/2 are Fig. 2. Stress resultants.
p = [PX,P,,
WI
Pzl’
shear and the global edge tractions P, are given in terms of the local edge tractions P., P, and P, as,
The generalised strain vector P corresponding to the stresses B of eqn (8) is,
(13c) where c = Cos B and s = Sin 0, see Fig 3. Y4 P"
where the linear component 8’ is, 0
il iI
go=
lP QO
0
(lob)
EJ
I
-----_)
and the nonlinear component fL is,
pL=
4
L
o
Substituting for u, u and w from eqn (1) into eqn (13a)
WC) and integrating the edge area integral explicitly through
0
The second integral in eqn (6) which represents the external virtual work due to body forces may be rewritten as, dW:x =
+duq, +dwqz)dz da
(11)
the thickness gives the work term d W:, as the sum of an area integral over a, and a line integral over an edge length s. dWb=~iT@dct~diT6,dr
(12a) where the midplane displacement vector i is,
and the generalised body force vector q is
iii, I’
(t4b)
in which typically # = 5 (r/2)P,, where - and t refer to forces on the top and bottom surfaces respectively. The generalised edge forces are
= [I
UW
(14a)
where B = [PX,Pw P*, k,
substituting for u, u and w from eqn (1) and integrating explicitly through the thickness allows eqn (11) to be expressed as an area integral giving,
ii = [ff, 8, I+,8x,ey1=
x
Fig. 3. Edgetractions.
II2
(P,‘, P,: P,‘, - Px’z, - Py’z) dz 42
T 1 .
(14c)
The virtual work eqn (6) may now be written entirely in terms of “mid plane” quantities as,
t/2
=
if
[I
_,,*
dqx, %*4.9 - qxz,- qyz)dz I =
(12c)
the body force vector pq is constant over the thickness
CONSlTNTtW
RELATIONS
For a linear elastic isotropic, homogeneous material with the stress in the z-direction o, = 0 the stresses u of
206
A.
PICA
et al. VARIATION OF STRAIN
eqn (5) are related to the Green’s strains E of eqn (3) by, a=Dr
(16a)
where the matrix D of constitutive coefficients is, iY Vl E
0 0
0 0
0 0
OOI-VO 2
D=m
00
Before proceeding with the discretisation of the virtual work eqn (15) it is necessary to consider further the variation of strain di due to the virtual displacements dii. Generally dP is given as the sum of the variation of the linear and nonlinear generalised strains as dB =di”tdBL.
0
(16b)
Since 6’ is a linear function of ii a
!$EO
0
000
ax dP” =
0
0
0
0
0
0
I
a I?( i$ ax !
0
d
0
0
/o -i
0
0
jo
0
0
P/ ax
0
Ja
0
a
0
-ay oio _a a ._...__...._.~___________“_____“x. 0
or
0
-I
1ay
0
-1
If the displacement gradients with respect to the lateral displacements t are,
in which A
/o
._____________I.________________________~
(17a)
Er Dp=j-qD;
0
oyo
oy
where E is Young’s modulus and tr is Poisson’s ratio. Substituting for E from eqn (3) and integrating explicitly through the thickness t enables the stress resultants b of eqn (9) to be written in terms of the generalised strains P of eqn (10) as,
(18)
-
Et’
.
D~=12(1_v2)D;
-
Gt
&=;I2
(17c)
e=
a3 ar;l T [ ax ’ ay 1
(20)
where
rlv
then the variation of the nonlinear component of the inplane strain is obtained from eqn (4d) in terms of the virtual gradients dg as,
07
G=E
2(1 t v) (17d)
depr. = Ae de
(2la)
where
and cx is a shear correction factor to allow for nonuniform shear stress distribution.
(21b)
and dezd
auatiT
[ ax ’ ay
(2lc)
1
FINITE ELEMENT FORMULATlON OF FQUILIBRIUM EQUATIONS The displacement ii within an element is given as a function of n discrete nodal displacements as,
(22a) LAGRANGIAN
IOLI
OS plus hlearchlcc, bubble functlon N,=ll-~‘111-~‘1
HETEROSIS
I OH I
05 OL wth w,. : 0
hlearchlcol
Fig. 4. Mindlin plate elements.
where Nj are shape functions, see Fig. 4, Is, is a 5 x matrix and the nodal displacements 6, are, si = [/ii, Bi,CJi,&i, 8,)I“.
5
unit
(224
Geometrically
nonlinear plate behaviour using a Mindlin formulation
For the Lagrangian and Heterosis elements a hierarchical formulation is adopted. Equation (22a) may be conveniently rewritten as
where
Boi = ii=Nii
207
(22c)
1
B& 0 0 B& [ 0 B&
where
md
in which component matrices are inplane (22d) aNi aNi
1 J ay-x
The virtual displacements di are written in terms of the nodal virtual displacements dg from eqn (22~) as di=Ndg
(2W
bending
where
IIsn St
dg=d
1
B:i =
f
alvi
(23b)
O -ay
-ax
shear The displacement gradients 8 of eqn (20) may now be written in terms of the nodal displacements 6 and cartesian derivatives of the shape functions as
aNi
,.,&=L-Ni ’
(244
8=G6
where
[$ The strain matrix
G=[C,,G*...Gnl and
(244
1
0
(26e)
ahri
-Ni
1
(26f)
’
BL consists of the nodal submatrices BL = [BL,, BLZ . . . B,_,I
Wd
where
oo%oo
(244
(26h)
The virtual gradients of equation (21~) are now written in terms of virtual nodal displacements as
Substituting eqns (23a) and (25) into (19) and (21a) respectively yields the strain variation in terms of the virtual nodal displacements dg as
G, = [
(-)()aNiO()’ ay
de=Gdg
The generaked Green’s strain vector P of eqn (lOa) is given in terms of nodal displacements 6, displacement gradients Ab and Cartesian derivatives of N as,
d = &+;B~(b)]g
(2W
where & is the strain matrix giving the linear straini lPorboand r,” and BL, which is linearly dependent upon 6, gives the non-linear strains 4‘. Consequently the nonlinear strains are quadratically dependent upon the nodal displacements 6. The constant matrix BO is written in terms of nodal submatrices as,
BO=PO,,Bm...Bonl
dP=Bdg
(25)
(26b)
(27a)
in which B = Bot BL(S).
(27b)
The virtual work eqn (15) is discretised, for an element, by substituting eqns (23a) and (27a) for dQ and di respectively giving (2ga) where the equivalent nodal load vector R due to body forces and tractions is NT@da + 1. NT fie ds. R=j-ANT4da+j-s
(28b)
Since the nodal virtual displacements d& are arbitrary
208
A. PICAet al.
the element nonlinear equilibrium equations become,
& = [o*, l72 (29)
e(6)= ( B*&da-R=O.
. .
. &IT.
(W
This enables the equilibrium eqns (29) to be expanded as
The load vector R may also contain nodal point loads. Equation (29) can have the dual role of representing either the element, or in an assembled form the total equihbrium equation. It is a nonlinear equation in 6 since B and C?are linear and quadratic functions of 8 respectively.
da-R=0
(35)
SOLUTION To N0NLlNJ.M EQuILlaRNM EQUATIONS The solution algorithms for the assembled nonlinear equilibrium eqns (29) are based on the Newton-Raphson method which consists of a series of linear solutions. If an initial estimate 6, for the total displacements gives residual (unbalanced) forces #(&) # 0 than an improved value W+r is obtained by equating to zero the linearised Taylor’s series expansion of S(L+,) in the neighbourhood of 6, as
#(&+I) z $(SI) + KTAL
= 0
(30)
Differentiating & with respect to 8, yields a typical term in KT as,
where (aRJa8,) = 0 for conservative loading. It is shown in the appendix that the tangent stiffness matrix can be written as,
where Kr is known as the assembled tangent stiffness matrix evaluated at Si and given by (31) Equation (30) is the linear incremental equilibrium equation which gives the linearised approximation to the relation between the residual forces and incremental displacements, A& at a point Si on the equilibrium path. The improved solution is then found as, St+, = 6, t A&.
Kr=/~BTBBdo+~_-GT[;Y
:]Gda.
(37)
Substituting for B, b and C from eqns (27b), (17b) and (24b) respectively gives, KT=&tKL+Krr
(38d
where the constant linear elastic stiffness matrix Ko is, &=/-&TB&dn
(38b)
(32)
Equations (30) and (32) represent the Newton-Raphson solution to the nonlinear eqns (29). The terms $($) and ASi give measures of the convergence of the solution To improve numerical stability and to give intermediate results, the load R is usually applied in increments. Since the Newton-Raphson method (N - R) requires repeated calculation and inversion of Kn a modified N-R method (KT2) is also employed whereby KT is calculated once only on the second iteration of each load increment.
the initial displacement matrix KL, which is quadratically dependent upon displacements 8, is KL =
IA
B,,‘BBL da t
IA
BL’&, da t
IA
BL=&B,_ da (38~)
and the initial stress matrix Id, is
K,=jjT[zy 2]Gde.
(384
TANGEhT SllFFNEss MATRIX The tangent stiffness matrix given in eqn (31) may be written as,
KT =
(33)
For computational efficiency it is necessary to consider the various nodal substitkess matrices contained in eqn (38a), this avoids much zero multiplication. Using eqn (26c) for the nodal submatrix Boi allows K,,, to be expressed in terms of inplane, bending and shear components as Koij =
where m is the number of degrees of freedom per element. For convenience the eight stress resultants of eqn (9) are rewritten as,
(39) where ij = 1 to n. From eqn (26c) and (26h) the nodal initial displacement
Geometricallynonlinearplate behaviourusing a Mindlinformulation
209
SI'RESS CALCULATION
matrix KLy is,
For all the Mindlin elements stresses are extrapolated from Gauss point values in a manner consistent with the integration rule used for the tangent stiffness matrix calculation. The extrapolation which results is given in Table 2. in which ki = AeG,
Examples (i) Preliminary. In the examples deflections, stresses and loads are given in non-dimensional form as follows,
Wb)
where the matrix Gib is,
aNi Gib =
_$ i
[ ay
w-w
0 0’ 1
The predominance of zeros in the submatrix Gi of eqn (24c) indicates that the nodal initial stress matrix Iij is best calculated explicitly giving the only non zero term at the 3 x 3 position in the 5 x 5 matrix as,
w = $(0,0)/t Q = (qa’)/(Et’) P = (pa*)/(Et’)
Central deflection Load (a) uniformly distributed (b) point Extreme fibre stress
flit-&Y) = u(x,y)a2/(Et2)
where a = characteristic plate dimension; q = uniformly distributed load; p =point load; i = x or y; U(X,y)= maximum extreme fibre stress at x,y. All examples use the KT2 modified Newton-Raphson approach and convergence is checked using a total residual norm criteria as,
(41)
EQUIVALENTNODALLOADSDUETOSTlWSSEs
To evaluate the residual nodal force vector +4(S) of eqn (29) the equivalent nodal forces due to the stress resultants 6 may be written for a typical node i as, Pi =
(RpOi)T&p da (k,,,’ kP + (B&)& •t (B&)T& I -CAL
(42)
INTEGRATION
The integrals in eqns (29) and (39-41) are evaluated using numerical integration. The orders of the Gauss quadratures used, see Table 1, are selective and depend upon whether the iteAmunder considera$on is a function of the membrane (D,, &,,). bending (I&,,&b) or shear (D,, &,) behaviour. The equivalent nodal load vector R due to the external forces is calculated using a 3 x 3 ru1e.t
where y = convergence tolerance (%). Figures in brackets in Tables 3-18 are percentage errors. (ii) Clamped square plate. This example has become a “standard” used by many authors[8,9,12] to check GNL thin plate formulations. The analytical thin plate solution is given by Levy [IS] who solved Von Karman’s equations using a double Fourier series; this solution is quoted as having a possible error of less than 2%. Since the problem is symmetric only a quarter mesh need be used, see Fig. 5. Tables 3-5 show results for central deflection, central and edge stress respectively. These provide a comparison between the four Mindlin elements and two thin plate elements obtained from a previous study[l2]. The thin plate elements are the Irons-Razzaque non-conforming triangular element[l6] (RI) with linear inplane behaviour and the slope conforming rectangular element introduced by Bogner et a/.[171 with parabolic inplane behaviour (BFS).
Table 1. Integrationrules LN
QS
QL
QH
MEMBRANE
2x2
2x2
3x3
3x3
BENDING
2x2
2x2
3x3
3x3
SHEAR
1x1
2x2
2x2
2x2
Table 2. Stress extrapolation LN
QL
OH
MEMBRANE
BILINEAR
BILINEAR
BIQUADRATIC
BIQIJADRATIC
BENDING
BILINEAR BILINEAR
BIQUADRATIC
BIQUADRATIC
SHEAR
CONSTANT
t2 X2 for Linear element. CAS Vol. II.No. SE
QS
BILINEAR BILINEAR
BILINEAR
210
Geometrically nonlinear plate behaviour using a Mindlin formulation Table 3. Clamped square plate-central deflections W
l-
FINITE ELEMENT
qs
LN
QH
qL
BFS
RI
0.2349
(0.87)
0.2351
(0.78)
0.2363
(0.08)
0.2361
(0.38)
0.2387
(0.73)
0.2361
(0.38)
0.4705
(0.10)
0.4673
(0.79)
0.4699
(0.23)
0.4687
(0.49)
0.4717
(0.16)
0.4685
(0.52)
0.6990
(0.57)
0.6887
(0.90)
0.6915
(0.50)
0.6902
(0.69)
0.6916
(0.48)
0.6900
(0.73)
(1.00)
0.9015
(1.15)
0.9008
(1.23)
0.9012
(1.18)
95.
0.912
0.9191
(0.78)
0.9003
(1.29)
0.9029
134.9
1.121
1.1318
(0.96)
1.1041
(1.51)
1.1063
(1.31)
1.1050
(1.43)
1.1025
(1.65)
1.1047
(1.46)
1.3350
(0.91)
1.2990
(1.81)
1.3009
(1.67)
1.2997
(1.76)
1.2961
(2.03)
1.2992
(1.80)
1.5347
(0.90)
1.4913
(1.95)
1.4928
(1.85)
1.4916
(1.93)
1.4879
(2.18)
1.4908
(2.00)
1.7439
(1.74)
i.6774
(2.13)
1.6786
(2.06)
1.6775
(2.13)
1.6744
(2.31)
1.6761
(2.21)
1.9109
(0.47)
1.8682
(1.77)
1.8555
(2.44)
1.8545
(2.50)
1.8529
(2.58)
1.8524
(2.61)
Table 4. Clamped square plate--centre stresses ux(O,O)
ri
LOAO Q
I?t
FINITE ELEMENT
EXACT (15)
BFS
RI
QH 2.6144
(0.55)
2.6890
(3.42)
2.6328
5.4520
(4.85)
5.4140
(4.11)
5.4047
(1.26) (5.47)
8.2908
(3.63)
8.0205
(0.26)
8.3372
(4.21)
17.7!
2.6
38.3
5.2
63.4
8.0
95.
11.1
11.575
(4.28)
11.115
(0.14)
11.103
(0.03)
11.066
(0.30)
10.521
(5.22)
11.127
(0.24)
34.9
I1 I1
13.3
14.567
(9.53)
13.817
(3.88) / 13.827
(3.96)
13.789
(3.67)
12.971
(2.47)
13.866
(4.26)
84.
15.9
17.493
(10.02)
16.461
(3.53)
16.497
(3.75)
16.456
(3.49)
15.390
(3.21)
16.552
(4.10)
i9.2
20.446
(6.49)
19.160
(0.21)
19.225
(0.13)
19.178
(0.12)
17.885
(6.85)
19.296
(0.50)
18.
21.9
23.605.
(7.78)
21.902
(0.00)
21.994
(0.43)
21.938
(0.17)
20.438
(6.67)
22.081
(0.83)
-02.
-25.1
26.318
(4.85)
24.805
(1.18)
24.780
(1.27)
24.713
(1.54)_
23.020
(8.29)
24.881
(0.87)
12
45.
I!
n
8.5570
(6.96)
8.3528
(4.41)
8.3258
(4.07)
Table 5. Clamped square plate-edge stresses u*((a/2,0)
r
LOAD Q --
1 ‘i 17.7:
38.3
LN
QS
5.48
3.3663
(38.6)
il.52
7.1122
(38.3)
5.3256
QL
f
(2.82)
11.155
(3.17)
5.3163
W
(2.99)
5.2991
RI5.2523
(3.30)
11.216
(2.64)
11.189
(2.87)
(4.16)
10.807
(6.19)
5.2915
(3.44)
11.137
(3.32)
63.4
18.03
11.198
(37.9)
17.515
(2.86)
17.726
(1.68)
17.697
(1.64)
16.680
(7.49)
17.556
(2.63)
95.
25.32
15.639
(38.2)
24.522
(3.15)
24.967
(1.39)
24.946
(1.48)
23.001
(9.16)
24.649
(2.65)
134.9
33.5
20.461
(38.9)
32.278
(3.65)
33.045
(1.36)
33.038
(1.38)
29.875
(10.8)
32.485
(3.03)
84.
42.4
25.597
(39.6)
40.710
(3.98)
41.885
(1.22)
41.900
(1.18)
37.257
(12.1)
40.961
(3.39)
245.
52.8
31.172
(41.0)
50.043
(5.22)
51.719
(2.05)
51.763
(1.96)
45.357
(14.1)
50.251
(4.83)
1318.
63.9
37.459
(41.4)
60.064
(6.00)
62.325
(2.46)
62.401
(2.34)
54.011
(15.5)
60.097
(5.95)
75.8
43.121
(43.11
70.897
(6.471
73.407
(3.16)
73.520
(3.01)
63.00
(16.9)
70.190
(7.40)
1402.
1
!
Y t r-------
,’
_,
/j
,’
,‘,
’
/’
’
’
,I
,’
x
I
I
I
I
I
I
L____
-_
t___2____-_
.
L
From the tables it can be concluded that the use of Mindlin formulation in the GNL context is an acceptable
@a_ /’
l
___-
_-__J
o= 300u-l 3,n E= 03 x10' Ibh' Y= 0316
t=
5. Square clamped and simply supported plate under uniformly distributed load-geometry and material properties.
Fii.
1
FINITE ELDIENT
EXACT (15)
T
alternative to the thin plate Kirchhoff formulation. Generally the 8 and 9 noded elements (QS, QL, QH, BFS) give comparable results for deflections and stresses whilst the elements having linear inplane behaviour (LN, RI) yield inferior stresses. A surprising result is the accuracy of the linear (LN) element for displacements. These results were obtained using a convergence tolerance y = 1%. (iii) Simply suppotied squareplate. Here the Mindlin formulation is compared to the thin plate Kirchhoff results obtained by Rushton[l8] who used a finite difference dynamic relaxation approach with a quoted error of less than 0.5%. The boundary conditions used are as follows, on x=ka/2;
u=v=w=tl,=O 0, free
(W
Geometrically nonlinear plate behaviourusing a Mindlin formulation
on y=+.a/2;
211
P
u=u=w=&=O 0, free.
Wb)
Plate geometry, material properties and mesh idealisation are shown in Fig. 5. The results for centre deflections and extreme fibre stresses are shown in Tables 6 and 7 respectively. A convergence tolerance of y = 1% was used. Again the 8 and 9 noded elements perform well in comparison with the finite difference results whilst the linear element gives less accurate, but acceptable, values. (iv) Clamped skew plates. Here the ability of the elements to mode distorted plates is explored. Firstly by reanalysing the clamped square plate, using the mesh shown in Fig. 6 (B = 0). This eleven Heterosis element mesh, for the quarter plate, gives displacements comparable with values obtained using the previous 16 element mesh. Inevitably, though, this coarse distorted mesh produces a degradation in stresses. Further distortion of the mesh allows two clamped skew plates @I= 30”, 45”)
Fig. 6. Clamped skew plate under uniformly distributed loadmesh details (see Fig. 5 for a, 1, E and v).
to be analysed. These latter problems have been discussed previously in Ref. [ 101where triangular Kirchhoff elements were employed, but no numerical results given. See Tables 8-10 for displacement, centre and edge stress results respectively. In these examples a convergence tolerance of y = 0.1% was used. (v) Clamped circular plate. Here the performance of
Table 6. Simply supported square plate--centre deflections W REF.
FINITE ELEMENT
1
(18)
LN I
1
I
QS
QL
I
I
QH I
0.335
0.2531 (5.40) 0.3478 (3.82)
0.3480 (3.84) 0.3478 (3.82)
0.818
0.8303 (1.51) 0.8184 (0.05)
0.8185 (0.06: 0.8184 (0.05)
1.47
1.4931 (1.57) 1.4655 (0.31)
1.4657 (0.29)
2.40
2.4380 (1.58) 2.3927 (0.30)
2.3932 (0.28) 2.3928 (0.30)
3.83
3.8985 (1.79) 3.8124 (0.46)
3.8134 (0.43) 3.8128 (0.45)
6.07
6.1571 (1.43) 6.0521 (0.29)
6.0539 (0.26) 6.0530 (0.28)
I
1.4655 (0.31)
Table 7. Simply supported square plate-centre stresses gx(0,O) FINITE ELEMENT
REF.
LOAD
Q
(18)' LN
9.16 36.6 146.5
QS
QL
QH
2.46
2.6354 (7.13) 2.6214 (6.56)
2.6029 (5.81) 2.6016 (5.76)
6.90
7.1837 (4.11) 7.0026 (1.49)
6.9826 (1.20) 6.9787 (1.14)
14.5
15.263
(5.26) 14.644
(1.00) 14.635
(0.93) 14.615
(0.79)
586.
30.0
31.396
(4.65) 30.183
(0.61) 30.188
(0.63) 30.130
(0.43)
2344.
65.2
68.428
(4.95) 65.673
(0.73) 65.756
(0.85) 65.661
(0.71)
9377.
148.3
(1.14) 149.87
(1.06)
155.24
(4.68) 149.66 (0.92)
Table 8. Clamped skew plates-centre
displacements W
FINITE ELEMENT
LOAD 17.79
149.99
6 = o" * 0.2291 (3.35)
B -5300
6 = 450
0.1412
0.0690
38.3
0.4576 (2.84)
0.2940
0.1473
63.4
0.6781 (2.43)
0.4585
0.2394
95.
0.8900 (2.41)
0.6309
0.3476
134.9
1.0947 (2.35)
0.8061
0.4703
184.
1.2906 (2.45)
0.9779
0.6013
245.
1.4838 (2.44)
1.1487
0.7388
318.
1.6708 (2.52)
1.3140
0.8762
402.
1.8488 (2.79)
1.4708
1.0083
*
error with respect to Ref. (15). see table 3.
I
212
A.
PICA et ai.
Table 9. Clamped skew elate-extreme fibre stresses at centre in direction normal to the edge FINITE ELEMENT
LOAD 6 =
o”
6 = 450
17.79
2.8929 (11.3)
1.5930 3.4775
38.3
6.0767 (16.9)
63.4
9.3357 (16.7)
5.7852
95.
12.590 (13.4)
a.5907
134.9
15.840 (19.1)
11.863
184.
19.074 (19.9)
15.425
245.
22.418 (16.7)
19.209
318.
25.845 (18.0)
23.016
402.
29.315 (16.8)
26.708
Table 10. Clamped skew plate-midside extreme fibre stresses in direction normal to the edge LOAD Q
FINITE ELEMENT 6 = 00
!
(12.2) (11.5)
B = 450
6 = 300 3.0991
1.6819
6.6972
3.6455 6.0640
(10.3)
10.974
(10.0)
16.038
(10.2)
21.925
12.875
(10.4)
28.524
17.386
(11.7)
35.959
22.761
(12.7)
44.020
(13.9)
52.452
9.1004
.....-. . Fish 28.868 35.504
’
.
-
b:2OOtn t =21n E :O 1x10@ lb/r? " =0.3
Fig. 7. Clamped circular and elliptical plates-geometry and material properties. the three “Q-family” elements (QS, QL, QH) is com-
pared with analytical Kiichhoff solutions for uniformly distributed[l9] and central point[20] loading. For both problems the curved boundary mesh shown in Fig. 7 is used. The results, using a convergence tolerance of O.l%, for displacements and stresses are given in Tables 11-15. Generally it is dithcult to distinguish a clear trend in the behaviour of the elements; whilst displacements compare well with the analytical solutions, the stresses are irregular, an element giving good accuracy for one load case and bad for the other. Nevertheless the Heterosis (QH) element consistently overestimates the edge stresses for both loading conditions. (vi) Clamped elliptical plate. Finally results are presented for the clamped elliptical plate under uni-
formly distributed load, shown in Fig. 7, these being compared with the analytical solution given in Ref. [19]. These results are shown in Tables 16-18 and were obtained using a 0.1% convergence tolerance. Centre displacements are in good agreement with the analytical results for all elements. For the QS and QL elements, the stresses are acceptable, the centre stresses being slightly more accurate. Again an irregularity is observed-the QH elements providing better stresses at the edge than at the centre. CONCLUSIONS
This paper has presented a geometrically nonlinear analysis of plates using a finite element Mindlin formulation. The paper contrasts the performance of the
Geometrically nonlinear plate behaviour using a Mindlin formulation
213
Table 11. Clamped circular plate (UDL)-centre displacements W LOAD
FINITE ELEMENT
REF.(19)
I Q
QS
QH
QL
1.
0.169
0.1614
(4.50)
0.1682
(0.47)
0.1644
(2.72)
2.
0.323
0.3116
(3.53)
0.3231
(0.03)
0.3167
(1.95)
3.
0.457
0.4453
(2.56)
0.4591
(0.45)
0.4514
(1.23)
6.
0.761
0.7566
(0.57)
0.7702
(1.21)
0.7626
(0.21)
10.
1.035
1.0417
(0.65)
1.0514
(1.58)
1.0459
(1.05)
15.
1.279
1.2950
(1.25)
1.3007
(1.70)
1.2974
(1.44)
Table 12. Clamped circular plate (UDL)-centre stress o,(O,O) r
,r REF. (19)
FINITE ELENENT
QL
QS
QH
1.
0.502
0.5018
(0.03)
0.5304
(5.66)
0.5119
2.
0.986
1.0044
(1.87)
1.0611
(7.61)
1.0270
(1.97) (4.15)
3.
1.421
1.4724
(3.62)
1.5524
(9.25)
1.5079
(6.17) (8.29)
6.
2.477
2.6070
(5.25)
2.7384
(10.5)
2.6824
10.
3.542
3.6774
(3.82)
3.8657
(9.14)
3.8077
(7.50)
15.
4.622
4.6595
(0.81)
4.9163
(6.38)
4.8567
(5.08)
Table 13. Clamped circular plate (UDLj-edge stresses u,(O,a) LOAD
Q
FINITE ELEMENT
REF.(19)
QL
QS
QH
1.
0.772
0.7831
(1.44)
0.7541
(2.32)
0.9352
(21.1)
2.
1.550
1.5550. (0.32)
1.5046
(2.93)
1.8553
(19.,7)
3.
2.302
2.2924
(0.42)
2.2265
(3.28)
2.7294
(18.6)
6.
4.324
4.2639
(1.39)
4.1799
(3.33)
5.0393
(16.5)
10.
6.517
6.4692
(0.73)
6.4094
(1.65)
7.5871
15.
8.694
8.8125
(1.36)
8.8271
(1.53)
10.267
(16.4) (18.1)
Table 14. Clamped circular plate (point load)-centre displacements W REF.(20)
T-
I ELEMENT FINIITI QL
QS
0.2130
0.2044
(4.04)
0.2149
(0.89)
0.2085
(2.11)
0.4052
0.3908
(3.56)
0.4088
(0.88)
0.3977
(1.86)
0.5705
0.5528
(3.10)
0.5754
(0.87)
0.5612
(1.62)
0.7123
0.6930
(2.71)
0.7184
(0.86)
0.7021
(1.44)
0.8354
0.8152
(2.41)
0.8427
(0.87)
0.8246
(1.29)
0.9442
0.9237
(2.17) 1
0.9527
(0.90)
I
Table 15. Clamped circular plate (point load)-edge stresses q(0.a) kINI_
LOAD REF. (20)
QH
QL
QS
Q 1.
0.4858
0.5318
(9.48)
0.4932
(1.52)
0.6767
(39.3)
2.
0.9592
1.0478
(9.24)
0.9765
(1.80)
1.3268
(38.3)
3.
1.3974
1.5240
(9.06)
1.4248
(1.96)
1.9194
(37.3)
4.
1.7988
1.9583
(8.87)
1.8355
(2.04)
2.4537
(36.4)
5.
2.1679
2.3553
(8.64)
2.2129
(2.07)
2.9377
(35.5)
6.
2.5110
2.7221
(8.41)
2.5635
(2.09)
3.3816
(34.7)
Linear, Serendipity, Lagrangian and Heterosis isoparametric elements by comparing results with those obtained from analytical and some other numerical thin plate solutions. The validity of the approach has been demonstrated for thin plates with rectangular, skew or curved boundaries. It was hoped that a discernable pattern of element behaviour would emerge from this investigation
but this proved not to be the case. For rectangular element meshes it would seem advisible to use the QH, Heterosis, element in view of its “safe” characteristics referred to in the introduction. However the presence of curved boundaries in the meshes led to inconsistent behaviour for all elements particularly with respect to stresses. Except in the case of the QH element this does not imply inaccurate results.
214
A. PICAet al. Table 16. Clamped elliptical plate--centre displacements W FINITE ELEMENT
REF.(19)
QL
QS 0.3449
(0.25)
0.3462
(0.64)
0.3449
(0.25)
0.600
0.6055
(0.92)
0.6078
(1.29)
0.6063
(1.04)
0.789
0.7985
(1.21)
0.8003
(1.43)
0.7990
(1.27)
1.163
1.1841
(1.82)
1.1845
(1.85)
1.1840
(1.81)
1.490
1.5085
(1.24)
1.5078
(1.19)
1.5080
(1.21)
1.790
1.7968
(0.38)
1.7952
(0.29)
1.7959
(0.33)
Table 17. Clamped elliptical plate-centre LOAD Q
QH
0.344
stresses +(O,O)
FINITE ELEMENT
REF. (19) _
QS
QL
QH
1.
0.862
0.9082
(5.36)
0.9021
(4.66)
0.9128
(5.89)
2.
1.550
1.6652
(7.43)
1.6575
(6.93)
1.6795
(8.36)
3.
2.077
2.2482
(8.24)
2.2384
(7.77)
2.2707
(9.32)
6.
3.217
3.4583
(7.50)
3.4450
(7.09)
3.4984
(8.75)
10.
4.333
4.5462
(4.92)
4.5328
(4.61)
4.6009
(6.18)
L 15.
5.490
5.5945
(1.90)
5.5873
(1.77)
5.6633
(3.16)
Table 18. Clamped elliptical plate-edge stresses ~~(0, a) FINITE ELEIlENT
REF. (19)
QS
The formulation presented herein has been extended to include shallow shells[21] allowing for the investigation of the postbuckling behaviour of plates with initial imperfections. This will be the subject of a further report. REFERENCES 1.
0. C. Zienkiewicz, The Finite Element Method. McGraw-Hill, New York (1978). 2. E. Hmton, A. Razzaque, 0. C. Zienkiewicz and J. D. Davies, A simple finite element solution for plates of homogeneous, sandwich and cellular construction. Proc. Inst. Ciu. Engrs. Part. 2,59,43-65 (Mar. 1975). 3. T. J. R. Hughes, M. Cohen and M. Haroun, Reduced and selective integration techniques in the finite element analysis of plates. Nuclear Engng Design 46.203-222 (1978). 4. T. J. R. Hughes, R. L. Taylor and W. Kanoknukulchai, A simple and efficient finite element for plate bending. Int. J. N&I. Meth. Engng ll(lO), 1529-43(1977). 5. E. D. L. Puuh, E. Hinton and 0. C. Zienkiewicz. A studv of quadrilateral plate bending elements with ieduced . integration. Int. J. Num. Meth. Engng 12, 1059-79(1978). 6. E. Hinton and N. Bicanic, A comparison of Lagrangian and Serendipity Mindlin plate elements for free vibration. comput.Smlctulvs10(3), 483-493 (1979). 7. T. I. R. Hughes and M. Cohen, The “Heterosis” finite element for plate bending. Comput. Structures 9(5), 445-450 (1978). 8. D. W. Murray and E. L. Wilson, Finite element large deflection analysis of plates. A.S.C.E. %(EMl). 143-165 (Feb. 1%9). _ 9. P. G. Bergan and R. W. Clough, Large deflection analysis of plates and shallow shells using the finite element method. Int. 1. Num. Meth. Engng $543-556 (1973).
QL
!
QH
10. R. S. Srinivasan and W. Bobby, Nonlinear analysis of skew plates using the finite element method. Cornput. Structure 6, 199-202 (1976). 11. A. K. Noor, M. D. Mathers and M. S. Anderson, Exploiting symmetries for efficient postbuckling analysis of composite plates. J. AIAA l&24-32 (Jan. 1976). 12. R. D. Wood, The Application of Finite Element Methods lo Geometrically Nonlinear structural Analysis. University of Wales, Swansea, U.K. C/Ph/20/73 (1973): 13. Y. C. Funa. Foundations of Solid Mechanics. Prentice Hall (l%S). -’ 14. R. D. Wood and 0. C. Zienkiewicz, Geometrically nonlinear finite element analysis of beams, frames, &hes and axisymmetric shells. Comout. Structures 7.725-35 (19771. 15. S. ievy, Square plate with clamped edges unde; normal pressure producing large deflections. NACA, Tech. Note 847 (1942). 16. A. Razzaque, Program for triangular bending element with derivative smoothing. Znt. J. Num. Meth. Engng 5, 588-9 (1973). 17 F. K. Bogner, R. L. Fox and L. A. Schmit, The generation of inter-element compatible sti&ess and mass matrices by the use of interpolation formulae. Proc. Conf. Matrix Melh. Srruct. Anal. Wright-Patterson A.F.B. (1965). 18. K. R. Rushton, Large deflexion of plates with initial curvature. Int. J. Mech. Sci. 12, 1037-1051(1970). 19. N. A. Weil and N. M. Newmark, Large deflections of elliptical plates. J. Appl. Mech. 23, 21-26 (1956). 20. R. Schmit, Large deflections of a clamped circular plate. A.S.C.E. 94(EM6), 1603-6 (Dec. 1968). 21. A. Pica, Geometrically Nonlinear Analysis of Mindlin Plates by Finite Hement Method. University of Wales, Swansea, U.K. C/M/132/78(Mar. 1978).
Geometrically
Substituting eqns (27b) and (A6) into (Al) gives
APPENLNX
Substituting eqn (36) into eqn (33) enables Kr to be written completely as, Kr =
I
[B’L(&.) A
215
nonlinear plate behaviour using a Mindlin formulation
do
t IbI(
r,=&ITbd~+f~M(A[~+~‘]dn.
(*7)
(Al) Considering the last integral in eqn (A7) and noting that Bo is not a function of 8 gives, upon expansion,
where
c M(&)B da =
and M(e)
= (A3)
‘a Ulas,
‘a was,
The differentials with respect to 8, in eqn (A8), are the coefficients of the G matrix of eqn (25) typically, for degree of
a “”
““as,
freedom i, Since the stress resultants d are a linear function of the strains i. see eqn (17a), the variation of stress may be written as, dd=ljdG.
(A9) (A4)
Substituting for d0 from eqn (27a) gives dd in the terms of the variation of 8 as, d& = SB d8.
(AS)
It is evident from eqns (A2) and (AS) that L(B) = fiB.
(A6)
and
at?!!.
(ab >
ay=G
21.
(AlO)
Consequently eqn (A7) may be rearranged to give the expression for Kr shown in eqn (37).