Comment on evaluation of shear forces and reactions from transverse shear deformations by using isoparametric quadratic mindlin plate elements

Compurers & Smcoues Printed in the U.S.A.

Vol. 19, No. 5/6, pp. W-903,

19&1

004s7949184 Fqamon

TECHNICAL

53.00 + .oo F?es Ltd.

NOTE

COMMENT ON EVALUATION OF SHEAR FORCES AND REACTIONS FROM TRANSVERSE SHEAR DEFORMATIONS BY USING ISOPARAMETRIC QUADRATIC MINDLIN PLATE ELEMENTS J. hOusEKt Swiss Federal Institute of Technology, Lausanne, Switzerland (Received 20 September 1983; received for publication 5 December

INTRODUCI’ION

A recent paper by Mukhopadhay[l] is concerned with evaluation of shear forces and reactions in plates in bending by using isoparametric quadratic elements. He presents an algorithm based on direct evaluation of shear forces from the element transverse shear deformations and applies the Kirchhoff formula known from the classical theory of thin plates for calculation of the plate reactions. The purpose of this note is: -to show that shear forces obtained from shear deformations of an isoparametric quadratic element are not reliable if the ratio I/h where h is thickness and 1 is element length or width, whichever is larger, is not sufficiently low, but that a significant improvement may be obtained by conveniently reducing transverse shear stiffness of the element. -to speculate on possible practical reasons for using the Kirchoff formula, inconsistent with the Mindlin plate theory. EVALUATION OF SHEAR FORCES FROM COMPUTED TRANSVERSE SHEAR DEFORMATIONS

Since the isoparametric plate element formulation takes into account the transverse shear deformation, it appears most natural to determine the transverse shear forces Q,, Q, directly from the corresponding transverse shear deformations ‘y, = &v/ax - @,, Y,,~= dw/& - S,,, where w stands for transverse displacement and @,, 19 are average rotations of the normals to the plate midsu rfyace. A transverse shear rigidity of 5/6 Gh is applied in the isotropic case to allow for non uniform shear distribution along the plate thickness h. The formulation is straightforward and it is likely that in the first years following publication of the Ahmad shell-plate, isoparametric element [2] its practical application was attempted (and abandoned!) by several research workers. Indeed, very soon it has been noticed that if displacement w and rotations 8,, 19, are interpolated within the element from nodal degrees of freedom, by using as usual, identical interpolation for each field

where the “shape” functions Ni = Ni(& ?I)are given polynomials of natural coordinates <,q, the transverse shear distribution produced is far from realistic. The reason for this behaviour is due to the fact that the degree of polynomial used to interpolate w is not sufficient to reproduce correctly the displacement field that corresponds to rotations 8,) 8, produced by bending of the plate. This error, naturally, tends to vanish if the degree of interpolation polynomials is progressively increased. In practice

TProfessor.

1983)

however the 8 node quadratic element is a near 100% choice because of its computational efficiency and its simplicity. It is well known that the displacement finite element procedure may be interpreted as a least squares error procedure where the measure of error is a weighted function of the error in stresses. As noted by Hermann[3], approximations which minimize the squared error tend to oscillate about the exact values, thus the stresses obtained from a finite element analysis would exhibit such an oscillatory nature. Care should therefore be taken to apply some technique of smoothing the stresses. Moreover, it may be expected that there exist unique “optimal stress’ points within certain finite element models at which the stresses or stress resultants have higher accuracy than at other points. In isoparametric elements the problem of finding such points is closely connected with the reduced integration technique, the beneficial effect of which (improved accuracy, significant attenuation of numerical dithculties known as “shear locking”) is now largely known and well understood (see e.g. [4,5] for references). More specifically, Barlow[6] has shown that for the 8 node quadratic element the 2 x 2 Gaussian integrating points are the optimum stress sampling points. Furthermore, Hinton and Campbell[7j have shown that the “local smoothing” consisting of a bi-linear extrapolation of the 2 x 2 Gauss point stress values to the element nodes, which are subsequently averaged at nodes common to several elements, is the natural method for sampling finite element stresses produced by quadratic elements using reduced integration. This, indeed, is probably also the most widely used technique for evaluating bending and twisting moments in isoparametric plate elements. The Barlow’s demonstration of the optimum stress points, unfortunately, does apply to transverse shear forces of the Mindlin mate model. Indeed. as deformations in the theory of elasti&y are obtained in terms of first derivatives of the cartesian displacement components, Barlow is concerned with such points which provide best accuracy for the first derivatives of a field. More specifically, he looks for such fixed points within the element where the 8 node Serendipity interpolation function, represented by 6 terms of a complete quadratic polynomial plus the two cubic terms <‘r~and rr12produces identical values of the first derivatives (with respect to < and 11) as a complete 10 terms cubic polynomial shearing the 8 given nodal values. He fmds that these are precisely the 2 x 2 Gauss points. The Mindlin theory, however, evaluates transverse shear by using underived fields of independent rotations S,,S,. Consequently the Marlow optimum points apply only to bending and twisting moments (calculated from the first derivatives of 8, and S,) but not to shear forces of an element. There are, however, other arguments for using the 2 x 2 Gauss points as sampling points for evaluating of shear forces. If the customary expression ]I” BTDB dxdv of stiffness matrix of an element is under-integrated with 2 x 2

899

900

Technical Note

Gauss points, the result is equal to that obtained by exact integration of the expression IsA fiTD&xdy in which the components of B are bi-linear approximations of corresponding component of B in the least square error sence[7]. Thus the values of shear obtained at 2 x 2 Gauss points are bi-linearly smoothed, hence more reliable. If the plate is not thick enough and the parasitic shear energy large, these values may still be in large error. On the other hand, in some special cases the results produced at 2 x 2 Gauss point may be exact even if the plate is thin. Indeed, it may be shown[8] that if a quadratic Timoshenko beam element is underintegrated with 2 Gauss points, then bending moment and transverse shear at the 2 Gauss points and dispiacement and rotation at the element ends match the exact values in cases where exact solution would require a cubic variation for the rotation and a quartic one for the transverse displacement. This surprising behaviour still holds true if one considers the case of cylindrical bending (v =0) of an under-integrated quadratic Mindlin plate element. Preceding considerations are heIpfu1 to understand why a procedure consisting in evaluation of shear forces from deformation at 2 x 2 Gauss points, their subsequent bilinear extrapolation to element nodes and averaging of values obtained from adjacent element may, with the same finite element mesh, produce results the quality of which ranges from excellent to catastrophic, depending on the load, boundary conditions and thickness involved. Figures 1-4 show some typical results obtained for simply supported7 and clamped square plates (Poisson ratio v = 0.3, span-to-thickness ratio a/h = 100) by using 3 x 3 quadratic elements over a symmetric quarter of the plate. It may be seen that the values prcduced (represented by small squares in the graphics), although not systematically wrong, bear no resemblance to the actual dis~bution in general. It is known that despite reduced integration, very thin elements can still fail because of numerical difficulties and the remedy used[9, IO] is to arbitrarily reduce the transverse shear stiffness of the elements. Thus, computed transverse shear deformation is bigger than strict theory requires yet ne~~bly small. A square plate with a span-to-thickness ratio of 100 and a 3 x 3 grid of quadratic elements over a symmetric quarter, does not need such an intervention as far as displacements and moments are concerned. The medicine may, however, be profitable to shear forces in the elements. Figures l-4 and Tables 1 and 2 show that reducing the transverse shear rigidity S = Gh/1.2 by a factor c of 10-20

L--a4

/

Fig. I. Distribution of Q, along y = 0 of a uniformly loaded simply supported square plate. 3 x 3 quadratic elements used with various values of the correction factor c.

if sufhcient to produce a significant improvement of transverse shear yet not too large to result in inacceptable changes in dehection and bending moment. Observing that the bending rigidity of a plate is a cubic function of the thickness h, .I%~ D = 12(1 - VZ)’ while the shear rigidity is linear in h enables to show that reducing S by a factor of c modifies the bending to shear stiffness ratio to value which corresponds to a I/;: times thicker plate, h’=hJI;.

(3)

This reduces significantly the parasitic shear energy while flexural stiffness is kept unchanged. With h =a/lOO, the element side of I = a/6 and the correction factor c = 10, the i to h’ ratio is equal to L = -!h’ h$c

= 5.270.

Assuming that a I/h’ ratio of 5 is necessary to keep the parasitic transverse shear at an acceptable level, yields the following simple expression for recommended value of the

TBoundary condition w = 8, = 0 is used for a good agreement with the classical theory of thin plates.

Table I. Simply supported square plate (u/h = 100, v = 0.3, w = I& = Bt = 0 at all edges). Variation of results with transverse shear correction factor c Uniformly

distributed

Cont. central

load

load

I

)

w center

M, cen-

Qx mid-

w center

: pa"/D

ter: pa*

edge: pa

:

edge: P/a

1

0.00406

0.0494

0.340

0.01154

5.111

10

0.00408

0.0495

0.337

0.01174

0.335

20

0.00410

0.0495

0.337

0.01192

I Theory

1)

Qx mid-

Paz/D

0.00406

0.0479

Theory of Kirchhoff,

2)

0.338

0.366

I 0.01160~~ m

Theory of Mindlin.

0.417

901

Technical Note correction factor c: c = (1/5@,

(4)

where I is the actual thickness and I the length or width of the element, whichever is larger. EVALUATION

OF REACTIONS AT SUPPORTED BOUNDARIES Mukhopadhay postulates[l] that shear force at a supported edge is different from reaction and that the latter should be determined by using the following expressions (Kirchhoff forces) known from the classical theory of plates:

This, of course, is incorrect since the Mindlin theory does not suffer from limitations of the classical theory and, in particular, allows to specify at a plate edge three independent boundary conditions. Thus e.g. at a free edge x = x0, = 0, as compared to classical theory with Q,=M,=M, only two con&tions specified as M, = R, = 0. Although the Kimhhoff forces have no place in the Mindlin theory of moderately thick plates, their application may in some cases be justified as a technique to circumvent difficulties arising from the slow convergence of results obtained in the neighborhood of an edge. This problem occurs, e.g. for simply supported edges with unrestrained rotations 8, at which w=M,=M,,=O.

@a)

The slow convergence is due to extremely steep variation of

-0.5

IQY.Pa

-0.L

-0.3 -0.2 -0.5 -0.1

a

0.0 0.0 0.1 0.2

!

1

q c-1 0 C=lO

1

3x3 quadratic

1.0

c

-10

A

C

120

-Theory

A C-20

i

a

!

0.3 ! 0.4

I

0.5

j

0.6

i Fig. 2. Distribution of QYalong y = 0 of a simply supported square plate with a point load at center. 3 x 3 quadratic elements used with various values of the correction factor c.

Fig. 3. Distribution of Q, along y = 0 of a uniformly loaded clamped plate. 3 x 3 quadratic elements used with various values of the correction factor c.

Table 2. Square plate with clamped edges (a/h = 100, v = 0.3, w = 8, = 8, = 0 at all edges). Variation of results with transverse shear correction factor c Uniformly S F

c =

w center

:

pa'/D

distributed

Mx center

:

load

Mx mid-

pa2

edge

:

pa'

Concentrated Q, midedge

:

pa

w center

:

Paz/D

central

Mx midedge

:

P

load Q

mid-

edXge : P/a

1

0.00121

0.0240

-0.0471

0.470

0.00537

-0.1253

-0.238

10

0.00129

0.0245

-0.0487

0.439

0.00575

-0.1248

0.670

20

0.00131

0.0246

-0.0487

0.435

0.00596

-0.1242

0.702

Theory

0.00127

0.0.229

-0:0513

0.441

0.00561;)) m

-0.1258

0.794

1)

Theory of Kirchhoff,

2)

Theory of Mindlin.

902

Technical Note

2.

0 c =lO

I

3.0

A c -20 -lJ=w

i

" LO--

cl .

so-/

00054 I

c

I

7.0- I' 8.0 Aa

06

05

9.0i

07

08

09

I

Fig. 4. Distribution of Q, along y = 0 of a clamped plate with a point load at center. 3 x 3 quadratic elements used with various values of the correction factor c.

Fig. 5. Distribution of A4, along y/a = 0.725 in a uniformly loaded simply supported square plate.

twisting moments and shear forces (Figs. 5 and 6)t as compared to the case of a simply supported edge with constraints tangent rotations 8 specified as

with (6a) is equal to the Kirchhoff force (Sa) of a classical thin plate (see Table 2 in [13]). Consequently, an application of Kirchhoff forces in the Mindlin plate context may be regarded as an artifice which enables converting shear reactions obtained with constraint tangential rotations (6b) to those associated with vanishing twisting moments (ba)--case suffering from poor convergence. If, however, the physical nature of supports is such that boundary rotations are actually constrained, using of Kirchhoff forces is irrelvant in either of the both plate theories. If the plate becomes thick, Mindlin theory with conditions (6a) produces increasingly larger deflections and bending moments (Table 3). Thus using (6b) instead of (6a) may give with (5a, b) correct reactions but deflections and moments will be in error. Undoubtedly, for an average program user, such a technique is rather misleading and dangerous. Let us finally suggest to those who are not discouraged and wish to experiment with the “Kirchhoff technique”, that evaluating of partial derivatives of yXY= iWYzin (Sa,b) as products of second derivatives of rotabons 8, and 8, by the plate shear stiffness (as used in Ref. [l]) is, probably, not the best approach. It is likely that differentiating the locally smoothed bi-linear fields of Mx, = M,,, given by 2 x 2 values at Gauss sampling points, would be computationally more efficient and would yield more accurate results.

w=M,=@,=O.

(6b)

This condition corresponds to that used in the Kirchhoff theory where w = 0 implies awlat = 8, = 0. It is worth mentioning that either of boundary conditions (6a,b) can be physically justified depending on the exact nature of the support. Except for the neighborhood of the edges, the results produced within the plate with either boundary condition are nearly equal if the plate is thin. The shear forces at supported edges (hence reactions) are, however, different. Whereas the midedge shear of 0.338pa obtained with (6b) is equal to that of the classical theory of thin plates, the shear force of 0.438 pa (Table 3) obtained

t Results presented have been obtained by using a non-uniform mesh of 10 X 10 cubic isoparametric elements at a symmetric quadrant together with an alternative method of evaluating shear forces as described in Ref. [ 111. Similar results obtained by segmentation method for a/h= 50 and v=O.3 may also be found in a recent paper[l2] by Kant and Hinton.

Table 3. Simply supported square plate, a/h = 10, Y = 0.15. Comparison of results of the theories of Mindlin and Kirchhoff Theory

Boundary

w=M Mindlin

n

50

w = Mn = M

Kirchhoff

w center

conditions

simply

t

:

=o

nt=

0

supported

pa'/D

Mx center

:

pa2

Q, midedge

:

0.004236

0.04236

0.338

0.004665

0.04576

0.438

0.004062

0.04236

0.338

pa

Vx midedge

:

0.438

pa

Technical Note I, :pa O.fI36

at alledges

0.6

0.7

ae

09

Fig. 6. Distribution of Q, along y/a = 0.5 in a uniformly loaded simply supported square plate. REFERENCES 1 M.

Mukhopadhay, Analysis of plates using isoparametric quadratic element-shear, reaction, patch loading and some convergence studies. Comput. Structures 17, 587-597 (1983).

903

2. S. Ahmad, B. M. Irons and 0. C. Zienkiewicz, Analysis of thick and thin shell structures by curved finite elements. Znt. J. Num. Meth. Engng 2, 419-451 (1970). 3. I. Hermann, Interpretation of tinite element procedure in stress error minimisation. J. AXE 98 (EM%. 1331-1336 (1972). 4. 0. C. Zienkiewicz, The Finite Element Method, 3rd Edn. McGraw-Hill. London (1977). 5. R. D. Cook, doncepts Mb Applications of Finite Element Analysis, 2nd edn. Wiley, New york (1981). 6. J. Barlow, Optimal stress locations in finite element modes. Int. J. Num. Meth. Engng, 10, 243-251 (1976). 7. E. Hinton and J. Campbell, Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Num. Meth. Engng 8, 461480 (1974). 8. J. Jirousek, On the capacity of the isoparametric quadratic beam element to solve accurately a uniformly loaded beam segment. Comput. Structures (submitted for publication). 9. E. D. L. Pugh, E. Hinton and 0. C. Zienkiewicz, A study of quadrilateral plate bending elements with reduced integration. Znt. J. Num. Meth. Engng 12, 1059-1079 (1978). 10. T. J. R. Hughes, R. L. Taylor and W. Kanoknukulchai, A simple and efficient finite element for plate bending. Int. J. Num. Meth. Engng 11, 1529-1543 (1977). 11. J. Jirousek and A. Bouberguig, A contribution to the evaluation of shear forces and reactions of Mindlin plates by using isoparametric elements. Comp. Structures (to appear). 12. T. Kant and E. Hinton, Mindlin plate analysis by segmentation method. J. Engng. Tech. ASeE, lob 537-556 (1983). 13. J. Jirousek, Basis for development of large finite elements locally satisfying all field equations. Comput. Merh. Appl. Mech. Engng 14, 65-92 (1978).