cofnpufers & slmcirrres Vol. 9, pp. 183-189 @ Pepmon Pnrs Ltd., 1978. Printed inGreat Britain
LARGE DEFLECTION ANALYSIS OF SKEW PLATES BY LUMPED TRIANGULAR ELEMENT FORMULATION D. N. BuRAGoHmt and S. C. PATODIS Department of Civil Engineering, Indian Institute of Technology, Bombay 400 076, India (Received 25 Febmary
1977; received for publication
16 September 1977)
Abstract-A finite difference scheme with triangufar mesh is presented for the analysis of skew plate problems with large defIections. The suggested fo~ulation is ~de~~ent of the boundary condition and uses energy principles to derive a set of nonlinear algebraic equations which are solved by using Ne~o~~ph~n iterative method with incremental loading. The inves?igation is concerned with the behaviour of constant thickness clamped and simply supported isotropic skew plates with immovable edges and subjected to uniformly distributed transverse toad. The effects of skew on plates with large deflections are investigated and comparisons are made with existing results; good agreement is shown. NOMENCLATURE
rectangular Cartesian coordinates outward normal directions to three sides of the triangle side of the skew plate the lengths of the element sides i-j, j-k and k - i respectively coordinates of point i displa~ment of point on middle plane of plate paraUe1to x, y and z axes, respectively bending moments and twisting moment inplane forces and shear area of the triangle flexural rigidity of the plate = Et’/12(1- v)* Young’s modulus of elasticity Poisson’s ratio thickness of the plate transverse deflection at the centre of the plate intensity of uniformly distributed load skew angle
boundary condition. The formulation is developed for geometric nonlinear analysis and is called as lumped triangular element (LTE) formulation. Graphical results are presented for the central deflection, bending moments and membrane forces for several skew angles, for a wide range of loads and for two edge conditions.
Referred to a fixed system of coordinates x and y and following von Karma& large deflection theory of plates, the displacements u, D, w define strains as[4]:
Symbols not given here are defined wherever they occur first.
Solution of large deflection problems of skew plates variously known as parallelogram plate, swept plate, rhombic plate and oblique plate. has not received much attention due to relatively difficult mathematical modelling. The nonlinear behaviour of isotropic clamped skew plates subjected to uniform pressure has been investigated only recently by Kennedy and Simon[l] using the perturbation technique and by Alwar !nd Rao[2] using the dynamic relaxation technique. In this paper, a simple and yet sufficiently accurate method, based on the finite difference energy principles, is presented for the analysis of skew plate problems. The six degrees of freedom constant strain triangle is used for cumulating membrane energy. The finite difference network adopted here for calculating bending energy is similar to the one suggested by YoshidaB] for plate bending analysis. The nodal degrees of freedom are the two membrane displacements II and u and the transverse displacement w. Additional degrees of freedom in the form of mid side normal rotations are introduced at the boundaries to make the formulation independent of the tAssistant Professor. SResearch Scholar.
in which the first column represents linear terms and the second gives nonlinear contribution in terms of the partial derivatives of the displacements u, u and w with respect to the coordinates x and y. Here subscript, “p” stands for inplane, “6” for bending and superscript, “L” stands for linear and “N” for nonlinear or large deflection contribution. A “comma” denotes derivative with respect to the subscript. The corresponding state of stress can be defined as
{4 = fmf
(2)
where [Dl is the etasticity matrix containing appropria~ material properties. Total strain energy U, by definition is U=~l{r)‘{~}dV=fI{r)~[D]{c}dV.
(3)
Now if the strain {E) can be related to displacements (8) by a relation matrix [B) as M= IBK@
(4)
D. N. BURAGOHAINand S. C. PATODI
184
The slopes of jk and ki along the 6 direction can be expressed in terms of tangential and normal rotations of the sides as
we can write,
U=~f{S}‘[8]T[D][8]{S}dV
(W,t)jk=~(wj-wk)sin~,-(w.,)jkcos6, I
(5)
= ; {#IKlGG
where [K] is known as stiffness matrix and as given by
and (W.f)ki
To estimate the strain energy or to calculate the strain displacement relation matrix [B], the continuum is divided into a set of rigid triangIes, connected by flexible nodes. Consider a typical interior element ijk with the elements whose one side is common with one of the three sides of ijk. The three neighbouring elements are ij1,jkm and kin as shown in Fig. 1. 1.4,u and w displacements are considered at each node to calculate membrane and bending energy. At the boundary lines, however, mid side normal rotation is also considered in addition to the nodal freedom to make the formulation independent of the boundary condition. Let 5, u and 6 be outward normat direction of the three sides of the triangle ijk. Other notations used are shown in Fig. 1. The second derivatives of normal dispIacement w in
(8)
=
$twiN
wk)
sin a1 - (wJ)ki cosa,.
(9)
After knowing (w.~)~, second derivative of w along 4 direction can be written as
Substituting eqns (8) and (9) in eqn (7) and then substituting eqn (7) in eqn (10) and after rearrangement, finally we can write
+ ski
cos* @&(w.r)ki - cos @p,sin Q>, wi
-cos~2sin~~wj+(cos~,sin~~ + cos & sin %)wkl*
(11)
x-y directions are obtained by first finding these derivaSimildy w,,.,, and w,r< can be expressed in terms of tives along the outward normal directions and then by the three normal displacements at i, j, k and three normal suitable transformation from local to global axes. rotations at the mid sides of the triangle. Let Q, R and S be the centres of sides i-j, j-k and k-i Now for a particular interior side ij, the corresponding respectively and corresponding normal rotations at Q, R normal rotation (w& can be expressed in terms of and S are denoted by (w& (w,,),, and (w+z)kirespecnodal displacements wi, wi, wk and w,. The variation of tively and let 6 axis, when produced, cut SR at P. first derivative of normal displacement w along side ij is Referring to Fig. 1, we can express the slope in the .$ constant and normal to it is at least linear. Hence normal die&on at point P in terms of slopes of jk and ki along displacement variation along the orthogonal axes 5 and y’ the 5 direction as (side i-j) will be (w+)p = t
[a& cos~&“,&k
+-hi
w=a,+a&+u,y’+a&*,
cos%tw,&d (7)
WI
We can find constants at, a2, . . . etc. by solving a set of
I
rx
Fig. 1. Lumpedtriangularelement.
185
Largedeflectionanalysis of ,skewplates simultaneous equations which will arise if the nodal conditions are inserted. At mid point Q, 8 = 0 and y’ = 0 and hence wI = a2 which after solution of the equations can be written as, (W&k
= (Y* =
(Fl + Fz)
Wj +
(Fl + F2)
Wj + F3 Wk + F4 Wl
(13) where
displacement vector {a} by a coefficient matrix [G] as
We know derivatives of w along the 5; q and 4 directions at mid side points Q, R and S respectively. For an interior side ij of the triangle ijk, at mid point Q, we can write
F, = -e,
( W.,)ij= ( W,e)fjCOS
Ilr F = XIC2Yt - x:Yk ’ a+,xk(xk -xl)’
sin @,
@c + y
and
(20) (w,J,, = (w,& sin Cp,- will, wicos 9,
where @‘cis the inclination of r axis with x axis. Substituting for (w*)ij from eqn (13) and after rearranging the terms we can write
and
Similar to eqn (13), expressions for other normal rotations (w.& and (w,~)~ can be obtained in cyclic order. Substrtuting values of rotations thus obtained in the expressions obtained for curvatures, we can finally write in matrix form
( w,)$~= (F,
cos
itt,
+
F2 cos
+ (F, COS@, - F2 COS@e +
sin Qp, -
Wk
W$
Wm
Wn}*
(15)
For a boundary element, say, with side jk along the boundary, the displacement vector would be: I&[’ =[a
Wj
wk
‘+‘I
(W.,)jk
WnI*
(16)
Thus we can
cb, + cos @$aij)Wi
F2 sin Cp, -
+ & sin @@wk +
Wj
t sin ~~aij)Wj
and
+ (F,
(Sb}T= {Wi
sin ~~aij~w~
FJ cos Q,w, + Fl, cos @,w,
(w,~), = (F, sin Qi, f, Fz sin
where (BL] is the coefficient matrix and can be automatically generated in the computer with the help of eqns (I 1) and (13) and {&} is the element displacement vector for bending and is given by
Qt, -
cos @Jaij)wi
F4 sin @,w,.
(21)
write for point Q, M, = [Gl&J
(22)
where [Glo is 2 x 6 matrix and can be obtained from*eqn (21). Finally the nonlinear strain displacement matrix [BbN] can be obtained at point Q from, WbNlo =
b%[~la
(23)
The curvatures in x-y diction can be obtained from the curvatures in the &q-i direction by suitable transformation. If [IL] is the transformation matrix, we can write
where [A], is the [A] matrix defined at point Q. Similarly [BbN] matrix can be obtained at other mid side points.
and hence
~~p~a~estrain d~~lacement matrix Considering triangular eiement ijk with nodal displacements B and v, strain displacement relation matrix can be obtained by using finite difference or fmite element method which can be finally written as [41:
[%I = P%Jb%~ = M~~l(4J
(17)
where matrix [&,I is of the order of 3 x 6 and is called as strain displacement matrix for bending. Nonlinear strain displacement matrix Nonlinear inplane strain contributions can be written as (eqn 1)
(18) where derivatives of w can be related to the bending
1 = IeJf4J
04
where bi = yj - yk, ct = xk -xr and other terms can be obtained . ... by a cyclic permutation of subscripts in the order rjR.
186
D.
N.
BURAGOHAIN
Govemi~g collinear e4ujiib~u~ equufio~~
If {$} represents the vector of the sum of the internal and external forces, we can derive the equilibrium equation by following the principle of virtual work and can finally write f4] {$}=\[I?]‘{u}dV-{R}=O
and
S. C.
PATODI
multiplication for c~culating element stiffness matrices. [K,] in eqn (27) is called as nonlinear stiffness matrix and is due to large displacements, given by
(30)
(25)
where {R) represents the external forces due to imposed loads and [B] is the strain disptacement matrix and is dependent on displacement {S} for a large deflection problem. Taking appropriate variation of eqn (25) with respect to (6}, the t~~ential stiffness matrix [K, 1 can be defined as t41, d{$) = W,ld@)
(26)
I&l = [&.I + K.J1+ [&I.
(23
where
Here [KJ is linear stiffness matrix and is given by
where [BbN] is calculated from eqn (23). Since [BbN] is linear in x and y and is known at mid side points of the triangle, nonlinear stiffness matrix can be calculated by integrating directly through its mid side values. For example,
(31) where [B] is defined at mid side points Q, R and S. [K,] in eqn (27) is known as initial stress matrix and is given by
[K-1 =I [%IIGI'PS][GJ dV where [I&,] and [Q,] are the elasticity matrices respectively for inplane and bending and are given by E
w%l=(f_vfj
lv
0
v *
0
I
and
14f
=
1y
12(1- v2) 0” ; [
0
(1 _;),2
I
.
where [G] is defined as per eqn (22) and [S] is the stress matrix and is given by (33)
[ 0 0 (I- 2412
Et3
(29)
where N;, NYand NXyare membrane forces. Now if the Newton-Raphson iterative procedure [4] is used for the solution of a set of equations represented by eqn (25), we can obtain correction in the trial solution for (i + 1)th iteration by (34)
A{S)i+r= - [KJi-‘I#}i For calculating linear stiffness matrix, [B,] and [BJ in eqn (28) are defined from eqns (24) and (17) respectively. Since these matrices are inde~ndent of the position within the element, integration can be replaced by area
and improved trial solution can be obtained from (35)
{Sh+l= {&Ii+ A{S)i+l. 24 SKEW
ANGLE
TOTAL-,-1
0
600
DEFLECTION
1600
2400
po*/ot
potot CENTRAL
6-O’
LINEAR f
_i
af
(32)
b)
STRESS
cr,
AT
Fig. 2. Variationof deflectionand stress with load, 6 = 0”.
tENTRE
3200
187
Largedeflectionanalysis of skew plates
In order to test the accuracy and usefulness of the method described, it is applied to the large deflection analysis of skew plates (with aspect ratio unity) subjected to uniform lateral pressure. Due to diagonal symmetry when all edges are clamped or simply supported, only one quarter of the plate is analysed and 6 x 6 discretization is used which results into 102 unknowns. Incremental Newton-Raphson iterative method is used for the solution of equations and Euclidean norm displacement criterion [5] is used with 1% cut-off for checking convergence. Results are presented for following two edge conditions: Fixed skew pl#fes
Results for large deflection of skew plates for skew angles B = o”, 15”, 30” and 45” are shown graphicaliy in Fig. 2-5 respectively. Four load increments, each of 800, are used to reach the final load value of 3200. Deflection
and stress variation are presented at the centre of the plate. Values of central deflection obtained by the present formulation (LTE) are compared with the values obtained by Kennedy and Simon[l] and Alwar and Rao [2]. Membrane, bending and total stresses computed at the centre of the plate are compared with those obtained by Alwar and Kao [21. It is clear from the figures that the deflection profiles obtained are almost identical to those obtained by Alwar and Rao[2] using dynamic relaxation technique with 20 x 20 mesh with a maximum difference of 5.31% for 45” skew. The values for 0” and 45” skew compare well with those obtained by Kennedy and Simon[l] using perturbation technique whereas vahtes obtained for 15”and 30” skews are comp~atively on higher side. Exceilent agreement is shown for the extreme fibre stress a,, calculated at the centre of the plate, with the values obtained by AIwar and Kao[2] for 0”, 15”, 30“ and 45” skew angles.
LINEAR
-
0
600
1600
2400
3200
po*/ot 01 CENTRAL
LTE
patot b)
DEFLECTION
STRESS
U7
AT
CENTRE
Fii. 3. Variationof deflectionand stress with load, 8 = 15”. 24r SKEW
0
600
1600
2400
3200
0
ANGLE
800
paVot 0)
CENTRAL
6 03s
1600
2400
po?ot DEFLECTION
b)
STRESS
ff+
AT
Fig. 4. Variationof deflectionand stress with load, B= 30”.
CENTRE
3200
D. N. BURAGOHAIN and S. C.
188 1.00
PATODI
24
SKEW
ANGLE
6
.45*
20 LINEAR
0.75
16 ~
“r; “
/
‘c as5 3u
*$a
12
S “0 x b
0
0.25 4
0
-
LTE
-
ALWAR
1600
800
2400
3200
0
800
1600
po~/ot .I
b)
DEFLECTION
CENTRAL
2400
3200
po*/ot STRESS
Or
AT
CENTRE
Fig. 5. Variation of deflection and stress with Ioad, B = 4.5’.
For the limiting case of 0” skew, the deflection results were also compared with Levy’s classical soiution[6] and it was found that the 6 x 6 LTE gives upper bound values with maximum discrepancy of 4.24% at the centre of the plate for the load value of 3200. Simply supported skew plates with immovable edges Variations of deflections, moments and membrane forces with skew angles and loads are presented in Fig. 6-g respectively. Only LTE results are presented since large deflection solutions for simpiy supported skew plates subjected to uniformfy ~s~but~ load are not readily available in literature. It can be observed from Fii. 6(a) that the large deflection curves tend to become linear and the maximum deflection at the centre decreases with the increase in skew. It is also interesting to examine the deflection ratio (w/w,) variation with change in skew angle.
0
800
1600 po4/0t
2400
Deflection profIles obtained, for a load value of 3200, along y axis are shown in Fig. 6(b). Figures 7 and 8 show the effect of skew on stress resultants. At the centre of a simply supported plate, the bending moment M, decreases very slowly and MY increases slowlv with the increase of skew. The rate of change of moments is of course different with respect to load. For any load, the membrane force N, falls steeply with increase in skew whereas NYincreases very slowly upto 15”skew and then decreases for higher skew angles.
A lumped triangular element formulation, developed herein based on the energy approach for the analysis of large deflection problems of skew plates, proves to be a convenient and effective matrix method for obtaining deflections and stresses. A comparison with the dynamic relaxation technique for the examples given herein in-
3200 8
4
8
DISTANCE o 1 CENTRAL
bf
DEFLECTION Fii.
6. Defiection
DEFLECTION
profileswith skew angle.
W/W,
2 ALONG ALONG
8 Y Y
AXIS
Large deflection analysis of skew plates
)_ !&7 I6yJb e
X/
c
k--a+
0.61
a i
IS SKEW a)
1
30 ANGLE
BENDING
povot
y
I
0
189
43
I
0
MOMENT
bf
hi,
30
IS SKEW
B
BENDING
ANGLE
d
0
MOMENT
My
Fit. 7. Variation of moments at centre with skew angle and load.
L0AO.U.
2.5
2.5
1
0
I5 SKEW a)
MEMBRANE
30 ANGLE FORCE
45
0
B
SKEW
N,
30
IS
b 1 MEMBRANE
ANGLE FORCE
45 8 NY
Fig. 8. Variation of membrane forces at centre with skew angle and load.
dicate very good a~eement even for relatively coarse division of the plate. The suggested formulation claims to be independent of the boundary condition (unhke finite difference formulation), involves very few unknowns, simplify the linear matrix manipulation by avoiding numerical integration, produces symmetric matrices, assures convergence to the true solution with finer discretization for the mesh generated by three sets of parallel lines and can be effectively processed on a digital computer. Acknowfedgement-The authors acknowledge with thanks the grant of a Senior Research FeLlowshipto Mr. S. C. Patodi by the Council of Scientific and Jndustrial Research.
1. J. 8. Kennedy and N. G. Simon, Linear and no&near analysis
of skewed plates. 3. Appt. Mech. Trans. ASME 34, 271-277 11967). 2. R. S. Alwar and N. R. Rao, Large elastic deformations of
clamped skewed plates by dynamic relaxation. Comprrt. stnfct. 4,381-398 (1974). 3. Y. Yoshida, Discrete triangular approximation of moment and displacement surfaces for plate bending analysis. pp. 139-182, U.S. Japan Conference on Recent Advances in Matrix Methods of Structural Analysis (1%9). 4. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill, London (1971). 5. C. Brebbii and J. Connor, Geometrically nonlinear finite element analysis. Proc. ASCE J. Engng Mech. Div. 95,463-483 (1%9). 6. S. Levy, Square plate with clamped edges under normal pressure producing large deSections. NACA Technical Report 740 (1942).