- Email: [email protected]
Four-node ‘flat’ shell element: Drilling degrees of freedom, membrane-bending coupling, warped geometry, and behavior
Computers & Strucrures Vol. 50, No. 4. pp. 549-W.
1994
Copyright Q 1994 Elswier Scii Ltd Printedin Great Britain.All rights rcwved co45-7949/w 56.00 + 0.00
Pergamon
FOUR-NODE ‘FLAT’ SHELL ELEMENT: DRILLING DEGREES OF FREEDOM, MEMBRANE-BENDING COUPLING, WARPED GEOMETRY, AND BEHAVIOR R. D.
CODK
Department of Engineering Mechanics and Astronautics, University of Wisconsin, Madison, WI 53706, USA (Received 18 August 1992) Alxstracf-A very simple way to obtain a quadrilateral shell element is to write membrane and bending stiffness matrices for a flat geometry and add them. This basic membrane-plus-bending element can be modified to incorporate membrane&ending coupling and to accommodate a geometry in which the four nodes are not coplanar. The present paper develops a 24 degree of freedom quadrilateral shell element and concurrently provides a status report on the formulation method by pointing out what works and
what does not. Numerical results are good: the element is not the best available four-node shell element in all test cases but is far from the worst.
INTRODUCMON
A finite element for shell analysis can be obtained by combining two flat elements, one for plane elasticity and the other for bending of flat plates. Early triangular elements were of this type. They did not work well in some problems because of poor membrane response and the absence of coupling between membrane and bending actions. Also, ridge-like intersections were awkward to model because there were but five degrees of freedom (DOF) per node. More recently, a membranebending coupling device has been devised [l, 21 and the advent of ‘drilling’ rotational DOF [3,4] has lessened the ridge intersection difficulty. Unfortunately, drilling DOF make membrane locking more likely, i.e. the model may become excessively stiff because rotational DOF associated with bending activate membrane strains as well as bending strains. Locking, however, can usually be reduced or even removed. As a result of all these efforts some very good triangular shell elements are now available [2, 5,6]. Similar procedures can be used to produce a four-node quadrilateral shell element having six DOF per node. However, the quadrilateral presents some difficulties that do not appear in the triangle, so that the quadrilateral is not a direct extension of the triangle. What follows is a report on the quadrilateral, with attention to procedures that succeed and those that fail. The principal shortcoming of the quadrilateral appears to be a tendency toward membrane locking, which however is not great, so that a good element is nevertheless produced.
ELEMENT
FORMULATION
Local coordinate system
The first step in element formulation is to place the element in a local coordinate system xyz that is oriented with reference to the element and is defined by using the global XYZ coordinates of element nodes (Fig. 1). Element properties are formulated in the local coordinate system. The last step before assembly of elements is transformation of nodal DOF from local to global directions. Use of a local coordinate system is convenient and also guarantees that the assembled finite element structure will display geometric invariance, regardless of the nature of the element stiffness formulation. Details of how the local coordinate system are established are not important. The procedure used here is similar to that described in [7]. The result is a local system xyz with xy the reference plane (Fig. 2). If the element is flat, its local nodal coordinates zi are all zero. If the element is warped, the zi are alternately + H and -H units above and below the reference plane xy. The x and y axes are parallel to sides of a rectangular element, i.e. x and y coincide with lines connecting midpoints of opposite sides. If the element is nonrectangular, x and y are approximately parallel to the lines connecting midpoints, with the longer of the two lines more heavily weighted. Except for membranebending coupling terms and the correction for a warped element, stiffness matrix formulation proceeds as if the element were flat, using 549
550
R. D.
COOK
brane locking when used in its curved form for shells, but more tendency toward shear locking when used in its flat form for membrane problems. However, shear locking can be avoided for most element shapes by defining the membrane shear strain as a constant over the quadrilateral [lo] (3) (a) (b) Fig. 1. (a) Quadrilateral shell element, showing global and local coordinate systems. (b) Nodal DQF in the local system XyZ. the plane quadrilateral that the element projects on the xy reference plane. Basic membrane formulation The membrane element is of the isoparametric type and is numerically integrated using a 2 x 2 Gauss rule. For membrane action there are three DOF at each comer i, consisting of x and y translations ui and r.+and drilling rotation Bzi,positive counterclockwise. Shape functions for the displacement field are obtained by starting with the 16 DOF of an element having both midside nodes and corner nodes but translational DOF only, then transforming in the manner described in [8]. In terms of the usual isoparametric coordinates t and r~, where - 1 G 5 < 1 and - 1 < r~Q 1, the resulting membrane displacement field that operates on comer DOF is
where A, is the area of a subtriangle. Here the quadrilateral is regarded as composed of the four overlapping subtriangles 123, 234, 341, and 412, with (~~r)~the shear strain in the k th subtriangle as defined by standard constant-strain angle operations [9]. In order to permit various options for shear strain to be tested, we define the final expression for yXYin the quadrilateral as Yxy =
B(r,)c,,
+
(1 -
/3)kJX~)(3,9
where (y,)(,, is produced by eqn (1) using standard manipulations, /I is a number 0 < /I -C 1, and (Ye,,) is the value of yXYdefined by eqn (3). Patch tests are passed for all values of 1. A modest but worthwhile improvement in membrane response results from the addition of two internal DOF a, and a*. Thus, u and v of eqn (la) and (lb) are augmented to become 24 =
(la)
+
(1 - {‘)a,
v = (lb) + (1 - q2)a2. Ni”i+(Y,4N*-Y2,N3)e,,
14= i
i=,
-~43hR3
+
(~43N,
-
(la)
~14WR4
(4)
W (5b)
Nodeless DOF a, and a2 are condensed before assembly of elements. Using eqns (5), plane rectangular elements give exact displacements in pure bending; using eqns (l), bending displacements are about 7% low. These results are unaffected by the value of /I used in eqn (4). In shell problems the effect of two internal DOF is less than 7%. Control of a mechanism
+
6,~
Ns -
-
XY NT N3
x23 Ns F’z,
+
+
h N, -
623 h’s
x41 %
Pz4
(lb)
where the first four Ni are the standard bilinear shape functions [e.g. Nr = (1 - <)(l - r/)/4], xi/ = Xi - Xj and y, = yi - yj are comer coordinate differences, and N,=h(l
-?)(l
-rt),
&=iLg(l +r)(l
--q2)
N,=+Jl
-r’)(l
+rt),
N,=&(l
-rt’).
-r)(l
Like all other elements having drilling DOF, a quadrilateral element has the instability in which all drilling DOF Bzihave the same value. A previously suggested control device [8] can be applied directly to the quadrilateral, but this leads to a very slight failure
(2)
Standard procedures [9] are now applied to produce a membrane stiffness matrix that operates on 12 nodal DOF. Drilling DOF are hierarchic, meaning that they can be omitted from the element formulation. If they are omitted the element has less tendency toward mem-
(a)
(b)
Fig. 2. (a) Reference element in the xy plane. (b) Element edge and DOF used to formulate membrane-bending coupling.
Four-node ‘flat’ shell element
of the patch test unless /? = 0 or /I = 1 in eqn (4). The q~d~laterai passes the patch test for all values of /? if instead the device is applied to the four overIapping subtriangles 123, 234, 341, and 412, and the four results are added. For the subtriangle spanning nodes 1, 2, and 3 for example, we add to the quad~lateral stiffness matrix the rank-one matrix k, = 10-6E.4tRR, where E = elastic modulus, A = subtriangle area, t = element thickness, and R
=
x23
Yz3
4A 4A
1 2 Y31 1 -xJZ ---YJZ --3 4A 4A 3 4A 4A
1
Plate-bending stiffness resists the transverse displacement DOF wi and the in-plane rotation DOF 0, and f?,. For plate action we choose the well-known DKT element [l 1, 121, arranged as a set of four overlapping subtriangles to avoid directional bias, and with each subtriangle having half the flexural rigidity of the actual material. This choice is made merely because the DKT element works well and is conveniently available. A good quadrilateral plate element could be used instead, probably with very similar results. c~upl~g
The coupling procedure is applied in various but similar forms in [l, 2,6]. The form used here is the very simple form described in 161,so it is presented ~ompac~y here. Along a typical element edge of length LI between nodes i and j, let s be an edgetangent coordinate and u, be the s-direction displacement, where u, varies linearly with s between nodal values u, and Use.Let the lateral displacement w along edge 0 be L,-s w=L,wi
q$
-%Wj+ ,
(0,
-
(c,),
e,>,
(7)
I
where Bi and 13~are nodal rotations in the xy plane, positive outward from the edge, and expressible in terms of fJIXi and f?,,and the angle $i between the x axis and direction ij. The average strain along edge v is taken as
+
2. (ei -
e,).
(9)
i
If standard strain transformations and eqn (9) are applied to each side of a triangular element we obtain relationships of the form c,=Ce
and
ez=SB,
(10)
where ss is a column matrix of three edge-parallel strains, the Cartesian strains are
C contains sines and cosines of edge orientations (pi , 8 contains six in-plane nodal rotation DOF, and S contains the information of eqn (9) for each edge along with sines and cosines of &. The expression L = C’S8 defines a constant strain state over the triangle. Terms of C’S are entered into the membrane strain-displacement matrix and make it possible for inextensional bending states to exist, thus removing a source of membrane locking. Unfortunately the foregoing procedure for a triangle is not directly transferrable to a quadrilateral. A least squares procedure that uses four edges to determine the three strains in L fails for a rectangle because it provides no info~ation about shear strain yXy.This defect can be cured by adding diagonals of the quadrilateral to the least squares procedure, but then a cylindrical shell element does not respond properly to moment vectors parallel to the axis of the cylinder. Both problems are solved by using the scheme of four overlapping triangles for all three strains in c = C’S& as described by yXyalone in eqn (3). Here the weighting by triangle areas A, is not necessary but provides a definite improvement in performance. The foregoing subtriangle procedure requires elevations hi at midpoints of both diagonals of the quadrilateral. They are computed as follows. Let h, , h2, h,, and hg be the hi at midpoints of sides 12, 23, 34, and 41. ~esumably these hi are available as input data. Let h, be hi at the midpoint of diagonal 24
where & and 5. A similar mid-diagonal are calculated where the latter integral contains the 6, expression of Mareuerre shallow shell theorvI. 1131. a The edee
= y I
1(6)
Bending stiffness
membrane-bending
elevation z may be taken as a parabola between nodes i and j with a mid-edge elevation hi relative to line g (Fig. 2). Thus eqn (8) yields
3
which operates on DOF [u, v, 0,, u2 u, 6,, u3 vj e,,]‘. The factor of 10e6 is rather arbitrary. A factor as large as lop2 changes element behavior very little [8]. The value 10m6 has almost no effect other than suppression of the mechanism.
551
qs are the 4 and tt coordinates of point calculation is used for h6, the other elevation, Isoparametric coordinates directly, e.g. CS=
(x2 + x4 w (XI + x2 + x3 + x4)/4 .
(12)
R. D. COOK
552
Fig. 3. View of adjacent shell elements 1 and 2 parallel to local xy planes.
The procedure of eqns (11) and (12) is approximate unless the element is rectangular or a parallelogram. A rigorous calculation of h, and h, would be much more tedious, and numerical evidence does not suggest that it would be worth the trouble.
Fig. 5. Shell roof test case. R = 25.0, L = 50.0, thickness = 0.25, E = 4.32 x lOa, v = 0, shell weight = 90.0 per unit area. A 2 x 2 mesh is shown. Vertical deflection at A is 0.3024 according to theory.
Warping correction
Thus far the element has been formulated using its projection, for which all nodes lie in the xy reference plane. If global coordinates do not place all element nodes in the same plane, element behavior is much too stiff in many problems unless modified by a warping correction. The correction used here is taken from eqns (9b), (SC), (lOc), and (11) of [7l. The correction is embodied in a matrix T such that the transformation TrIcT produces a z-direction stiffness from membrane stiffness coefficients while leaving existing membrane stiffness coefficients and bending stiffness coefficients unchanged. Matrix T is generated entirely from nodal coordinate data in the local xyz system. Membrane locking tendency
The origin of membrane locking is suggested by Fig. 3, where a nodal rotation 0, is associated with bending deformation in element 2 but has a drilling component 6, normal to the base plane of element 1; 0, is associated with membrane action in element 1. Accordingly eb is in part resisted by membrane stiffness, which is typically far greater than bending stiffness and hence causes bending deformations to be much too small. This ‘locking’ tendency vanishes as elements become coplanar. The same difficulty appears in a triangular shell element [6]. The element of [6] contains a form of constant-strain triangle and requires two stabilization matrices: one is described by eqn (6) and the other is
a rank-two matrix that operates on only the drilling DOF. The latter matrix is the sole cause of membrane locking in shells, but must be present when elements are coplanar. This matrix can be multiplied by a factor that is unity for coplanar elements but approaches zero as the curvature of the element increases. Thus it is easy to avoid membrane locking in the triangular shell element. The corresponding remedy for the quadrilateral would use /_I= 1 in eqn (4) for coplanar elements and arrange for /I to approach zero as elements become more highly curved. However, edge tractions on elements produce nodal torques corresponding to the 0,, when /I = 1, but zero torques when #I = 0. Hence adjacent elements having different values of /I would be incompatible (with the triangular element, similar adjustments of the second stabilization matrix present no such difficulty). At this time the writer does not know of a satisfactory cure, but fortunately the quadrilateral element does not have a great propensity for membrane locking. It may be noted that use of one-point integration for only the membrane terms associated with the tlli, with four-point integration for other terms, produces elements that either fail the patch test or behave badly otherwise.
1.0
.A__
OF
1.0
I I-
L
“.”
lW
OF
0.9
I.1 1.9
Diaphragm
3.1 3.9
2.1 2.9
4.1 4.9
5.1
.I
I.0 t 1.0
Fig. 4. Plane beams of Poisson ratio v = 0.3 loaded in pure bending. Drawing is not to scale.
Fig. 6. Pinched cylinder test case. R = 300.0, L = 600.0, thickness = 3.0, E = 3.0 x 106, v =0.3, F = 1.0. A 2 x 2 mesh is shown. Displacement at A is 0.18248 x 10e4 according to theory.
Four-node ‘tlat’ shell element
J
553
Fig. 8. Cantilevered strip with end-to-end twist of 90”. L = 12, b = 1.1, r = 0.0032, E = 29 x 106,v = 0.22. A 2 x 4 mesh is shown. Load F = 10e6 applied as forces 0.25 x 10m6, 0.50 x 10e6, and 0.25 x 10e6 at the three end nodes. Displacement at C is 5256 x 10e6 according to theory.
Free
+F
X Fig. 7. Hemispherical shell test base. One quadrant with a 4 x 4 mesh is shown. Radius = 10.0, thickness = 0.04, E = 6.825 x lo’, v = 0.3, F = 1.0 for the quadrant modeled. Displacement at A is 0.0924 according to theory. NUMERICAL RRSULTS
Membrane behavior The two beam problems of Fig. 4 are standard test cases [8,14]. The arrangement of trapezoidal elements is a test case in which four-node quadrilaterals often display shear locking. Expressed as the ratio of computed tip deflection w to exact tip deflection w,,, the results are rectangular
mesh, fl = 0:
w/w,, = 1.000
rectangular
mesh, /I = 1:
w/w,, = 1BOO
distorted mesh, fl = 0:
w/w,, = 0.051
distorted mesh, /I = 1:
w/w, = 0.897.
Clearly, for rectangular elements j does not matter, while badly shaped elements suffer unless fl = 1. With j? = 1, the present membrane formulation appears to have almost the same behavior as the membrane element suggested in [14]. Shell test cases The standard ‘obstacle course’ test cases [15] are depicted in Figs 5-8. Data for these problems are stated in the captions. In Figs 5 and 6, nodes at a diaphragm support are assumed free to rotate about a normal to the diaphragm, in accord with the simple support condition recommended in [16]. In Fig. 7 a small hole of radius 0.175 is placed at the top, where all DOF are restrained. All meshes ate uniform
Table 1. Shell test case results reported as the ratio of computed deflection to accepted theoretical deflection Element formulation
Mesh
Flat /3=0
Flat /?=I
5.085 1.548 1.084 1.014
CUNtd
&Ned
j=o
/I=1
4.941 1.541 1.085 1.014
2.951 1.109 1.010 0.995
2.471 1.080 1.009 0.996
0.012 0.098 0.630 0.938
0.012 0.097 0.627 0.931
2.125 1.150 0.941 1.016
0.015 0.178 0.775 0.990
Hemisphere, Fig. 7 1x1 2x2 4x4 8x8
0.117 0.852 1.002 1.001
0.056 0.228 0.444 0.948
1.533 0.654 0.586 0.971
0.057 0.168 0.606 0.962
Twisted strip, Fig. 8 1x1 2x2 2x4 2x8
0.916 0.859 0.978 0.999
0.122 0.608 0.969 0.993
Shell rooJ Fig. 5 1x1 2x2 4x4 8x8 Pinched cylinder, Fig. 6 1X1
2x2 4x4 8x8
t Element sides are straight in a rectangular mesh.
r;
1;
554
R. D.
COOK
Table 2. Shell test case results using the ‘curved softened’ TRIANGULAR Mesh
Shell roof
IX1
moo
2x2 4x4 8x8
1.406 0.880 0.939
element of [6]
Pinched cylinder
Hemisphere
Twisted strip
0.121 0.959 0.974 0.999
0.213 0.873 0.978 1.007
0.775 0.962t 0.992$
t 2 x 4 mesh. $2 x 8 mesh.
except those used for the hemisphere, where smaller elements appear at the top. In Fig. 8, note that f = 0.0032 rather than the less demanding case of t = 0.32 that also appears in [15]. Results appear in Table 1. The ‘flat’ formulation is obtained by omitting the membrane-bending coupling device, i.e. by setting hi = 0 in eqn (9) for all midpoints. We see from Table 1 that /I = 1 provides a stiffer formulation that /_l= 0, but the tendency toward membrane locking provided by /? = 1 is small. For the shell roof and pinched cylinder cases, the best element formulation appears to be curved with /I = 1. For the pinched cylinder the curved element with fl = 0 appears to converge faster; however in this problem the zone of greatest deformation is very localized around the load points, so the apparent good convergence properties of this element may instead be a symptom of a too-flexible formulation. In the hemisphere test case there is considerable bending and rigid body motion of elements but little membrane strain. For this case the best element is flat and without drilling DOF. Similar behavior is observed with triangular elements [2]: flat constantstrain triangles without drilling DOF are very accurate in the hemisphere problem. The twisted strip problem clearly shows the need for the warping correction if elements are thin. Without the warping correction, deflection ratios in the last line of Table 1 become 0.0001 and 0.0003. The regular meshes in Figs 5-7 do not result in warped elements. With the pinched cylinder of Fig. 6, numerical tests in which nodes were relocated so that elements became warped showed that the warping correction sometimes stiffened and sometimes softened the response but always made the computed result more nearly exact. In comparison with other available four-node quadrilaterals [17,18] the present curved element with /I = 1 appears to be more accurate in the shell roof and pinched cylinder problems but less accurate in the hemisphere problem. The twisted strip problem is not used in [17, 181. For comparison with results given by a triangular element [6] we refer to Table 2. Triangular elements are produced by dividing each quadrilateral along its shorter diagonal. A two-element patch of triangular elements has the same number of DOF as one quadrilateral element, so that the corresponding meshes in Tables 1 and 2 use the same number of DOF.
CONCLUDINGREMARKS This paper has considered shell elements obtained by the very simple process of combining standard membrane and bending formulations with a device for membrane-bending coupling and a device for inclusion of warping effects. The process leads to reasonably good four-node quadrilaterals and threenode triangles. Quadrilaterals have better membrane behavior than triangles and are free of the directional bias that usually accompanies triangles, e.g. when triangles form a cross-hatch pattern rather than a herringbone or union-jack pattern. Also, quadrilaterals do not have the ‘free comer’ problem present in triangles, in which a curved element may display large and erroneous displacements if it has a loaded node that is not shared by another element. If used without internal DOF, static condensation is avoided, while a quadrilateral built of four non-overlapping triangles would have six internal DOF to be condensed before assembly of elements and would require additional data to properly locate the internal node on the shell surface. On the other hand, two triangles are needed to replace one quadrilateral, so there is more interelement membrane-bending coupling with triangles than with quadrilaterals for a given number of nodes, the smaller span of the triangle reduces the tendency toward membrane locking. If the tendency for membrane locking in the quadrilateral can be reduced it may prove to be a very effective element.
REFERENCES 1. H. Stolarski, T. Relytschko, N. Carpenter and J. M. Kennedy, A simple triangular curved shell element. Engng Comput. 1, 210-218 (1984). 2. N. Carpenter, H. Stolarski and T. Belytschko, Improvements in 3-node triangular shell elements. Int. J. Numer. Meth. Engng 23, 1643-1667 (1986). 3. D. J. Allman, A compatible triangular element including vertex rotations for plane elasticity analysis. Comput. Strut. 19, l-8 (1984). 4. N. Carpenter, H. Stolarski and T. Relytschko, A flat triangular shell element with improved membrane interpolation. Comm. Appl. Numer. Meth. 1, 161-168 (1985).
5. J. Fish and T. Belytschko, Stabilized rapidly convergent 18-degrees-of-freedom flat triangular shell element. Inr. J. Numer. Meth. Engng 33, 149-162 (1992). 6. R. D. Cook, Further development of a three-node triangular shell element. Inr. J. Numer. Meth. Engng. 36, 1413-1425 (1993).
Four-node ‘flat’ shell element 7. B. P. Naganarayana and G. Prathap, Force and moment corrections for the warped four-node quadrilateral plane shell element. Comput. Struct. 33, 1107-1115 (1989). 8. R. H..MacNeal and R. L. Harder, A refined four-noded membrane element with rotational degrees of freedom. Comput. Struct. 28, 75-84 (1988). 9. R. D. Cook, D. S. Malkus and M. E. Plesha, Concepts and Applications of Finite Element Analysis, 3rd Edn. John Wiley, New York (1989). 10. R. D. Cook, Avoidance of parasitic shear in plane element. AXE J. Struct. Div. 101, 1239-1253 (1975). Il. J. A. Stricklin, W. E. Haisler, P. R. Tisdale and R. Gunderson, A rapidly converging triangular plate element. AIAA Jnl7, 180-181 (1969). 12. C. Jeyachandrabose, J. Kirkhope and C. R. Babu, An alternative formulation for the DKT plate-bending element. ht. J. Numer. Meth. Engng 21, 1289-1293 (1985). 13. K. Marguerre, Zur theorie der gekriimmten plate grosser formlnderung. Proc. Fifih Int. Cong. for Appl.
14.
15.
16.
17.
18.
555
Me& (Edited by J. P. Den Hartog and H. Peters), pp. 93-101. John Wiley, New York (1939). M. Iura and S. N. Atluri, Formulation of a membrane finite element with drilling degrees of freedom. Comput. Mech. 9, 417-428 (1992): _ T. Belytschko, B. L. Wong and H. Stolarski, Assumed strain stabilization procedure for the 9-node Lagrange shell element. Int. J. Numer. Meth. Engng 28, 385414 (1989). 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). W. K. Liu, E. S. Law, D. Lam and T. Belytschko, Resultant-stress degenerated-shell element. Comput. Meth. Appl. Mech. Engng 55, 259-300 (1986). K. Y. Sxe and C. L. Chow, A mixed formulation of a four-node Mindlin shell/plate with interpolated covariant transverse shear strains. Comput. Struct. 40, 775-784 (1991).
1994
Copyright Q 1994 Elswier Scii Ltd Printedin Great Britain.All rights rcwved co45-7949/w 56.00 + 0.00
Pergamon
FOUR-NODE ‘FLAT’ SHELL ELEMENT: DRILLING DEGREES OF FREEDOM, MEMBRANE-BENDING COUPLING, WARPED GEOMETRY, AND BEHAVIOR R. D.
CODK
Department of Engineering Mechanics and Astronautics, University of Wisconsin, Madison, WI 53706, USA (Received 18 August 1992) Alxstracf-A very simple way to obtain a quadrilateral shell element is to write membrane and bending stiffness matrices for a flat geometry and add them. This basic membrane-plus-bending element can be modified to incorporate membrane&ending coupling and to accommodate a geometry in which the four nodes are not coplanar. The present paper develops a 24 degree of freedom quadrilateral shell element and concurrently provides a status report on the formulation method by pointing out what works and
what does not. Numerical results are good: the element is not the best available four-node shell element in all test cases but is far from the worst.
INTRODUCMON
A finite element for shell analysis can be obtained by combining two flat elements, one for plane elasticity and the other for bending of flat plates. Early triangular elements were of this type. They did not work well in some problems because of poor membrane response and the absence of coupling between membrane and bending actions. Also, ridge-like intersections were awkward to model because there were but five degrees of freedom (DOF) per node. More recently, a membranebending coupling device has been devised [l, 21 and the advent of ‘drilling’ rotational DOF [3,4] has lessened the ridge intersection difficulty. Unfortunately, drilling DOF make membrane locking more likely, i.e. the model may become excessively stiff because rotational DOF associated with bending activate membrane strains as well as bending strains. Locking, however, can usually be reduced or even removed. As a result of all these efforts some very good triangular shell elements are now available [2, 5,6]. Similar procedures can be used to produce a four-node quadrilateral shell element having six DOF per node. However, the quadrilateral presents some difficulties that do not appear in the triangle, so that the quadrilateral is not a direct extension of the triangle. What follows is a report on the quadrilateral, with attention to procedures that succeed and those that fail. The principal shortcoming of the quadrilateral appears to be a tendency toward membrane locking, which however is not great, so that a good element is nevertheless produced.
ELEMENT
FORMULATION
Local coordinate system
The first step in element formulation is to place the element in a local coordinate system xyz that is oriented with reference to the element and is defined by using the global XYZ coordinates of element nodes (Fig. 1). Element properties are formulated in the local coordinate system. The last step before assembly of elements is transformation of nodal DOF from local to global directions. Use of a local coordinate system is convenient and also guarantees that the assembled finite element structure will display geometric invariance, regardless of the nature of the element stiffness formulation. Details of how the local coordinate system are established are not important. The procedure used here is similar to that described in [7]. The result is a local system xyz with xy the reference plane (Fig. 2). If the element is flat, its local nodal coordinates zi are all zero. If the element is warped, the zi are alternately + H and -H units above and below the reference plane xy. The x and y axes are parallel to sides of a rectangular element, i.e. x and y coincide with lines connecting midpoints of opposite sides. If the element is nonrectangular, x and y are approximately parallel to the lines connecting midpoints, with the longer of the two lines more heavily weighted. Except for membranebending coupling terms and the correction for a warped element, stiffness matrix formulation proceeds as if the element were flat, using 549
550
R. D.
COOK
brane locking when used in its curved form for shells, but more tendency toward shear locking when used in its flat form for membrane problems. However, shear locking can be avoided for most element shapes by defining the membrane shear strain as a constant over the quadrilateral [lo] (3) (a) (b) Fig. 1. (a) Quadrilateral shell element, showing global and local coordinate systems. (b) Nodal DQF in the local system XyZ. the plane quadrilateral that the element projects on the xy reference plane. Basic membrane formulation The membrane element is of the isoparametric type and is numerically integrated using a 2 x 2 Gauss rule. For membrane action there are three DOF at each comer i, consisting of x and y translations ui and r.+and drilling rotation Bzi,positive counterclockwise. Shape functions for the displacement field are obtained by starting with the 16 DOF of an element having both midside nodes and corner nodes but translational DOF only, then transforming in the manner described in [8]. In terms of the usual isoparametric coordinates t and r~, where - 1 G 5 < 1 and - 1 < r~Q 1, the resulting membrane displacement field that operates on comer DOF is
where A, is the area of a subtriangle. Here the quadrilateral is regarded as composed of the four overlapping subtriangles 123, 234, 341, and 412, with (~~r)~the shear strain in the k th subtriangle as defined by standard constant-strain angle operations [9]. In order to permit various options for shear strain to be tested, we define the final expression for yXYin the quadrilateral as Yxy =
B(r,)c,,
+
(1 -
/3)kJX~)(3,9
where (y,)(,, is produced by eqn (1) using standard manipulations, /I is a number 0 < /I -C 1, and (Ye,,) is the value of yXYdefined by eqn (3). Patch tests are passed for all values of 1. A modest but worthwhile improvement in membrane response results from the addition of two internal DOF a, and a*. Thus, u and v of eqn (la) and (lb) are augmented to become 24 =
(la)
+
(1 - {‘)a,
v = (lb) + (1 - q2)a2. Ni”i+(Y,4N*-Y2,N3)e,,
14= i
i=,
-~43hR3
+
(~43N,
-
(la)
~14WR4
(4)
W (5b)
Nodeless DOF a, and a2 are condensed before assembly of elements. Using eqns (5), plane rectangular elements give exact displacements in pure bending; using eqns (l), bending displacements are about 7% low. These results are unaffected by the value of /I used in eqn (4). In shell problems the effect of two internal DOF is less than 7%. Control of a mechanism
+
6,~
Ns -
-
XY NT N3
x23 Ns F’z,
+
+
h N, -
623 h’s
x41 %
Pz4
(lb)
where the first four Ni are the standard bilinear shape functions [e.g. Nr = (1 - <)(l - r/)/4], xi/ = Xi - Xj and y, = yi - yj are comer coordinate differences, and N,=h(l
-?)(l
-rt),
&=iLg(l +r)(l
--q2)
N,=+Jl
-r’)(l
+rt),
N,=&(l
-rt’).
-r)(l
Like all other elements having drilling DOF, a quadrilateral element has the instability in which all drilling DOF Bzihave the same value. A previously suggested control device [8] can be applied directly to the quadrilateral, but this leads to a very slight failure
(2)
Standard procedures [9] are now applied to produce a membrane stiffness matrix that operates on 12 nodal DOF. Drilling DOF are hierarchic, meaning that they can be omitted from the element formulation. If they are omitted the element has less tendency toward mem-
(a)
(b)
Fig. 2. (a) Reference element in the xy plane. (b) Element edge and DOF used to formulate membrane-bending coupling.
Four-node ‘flat’ shell element
of the patch test unless /? = 0 or /I = 1 in eqn (4). The q~d~laterai passes the patch test for all values of /? if instead the device is applied to the four overIapping subtriangles 123, 234, 341, and 412, and the four results are added. For the subtriangle spanning nodes 1, 2, and 3 for example, we add to the quad~lateral stiffness matrix the rank-one matrix k, = 10-6E.4tRR, where E = elastic modulus, A = subtriangle area, t = element thickness, and R
=
x23
Yz3
4A 4A
1 2 Y31 1 -xJZ ---YJZ --3 4A 4A 3 4A 4A
1
Plate-bending stiffness resists the transverse displacement DOF wi and the in-plane rotation DOF 0, and f?,. For plate action we choose the well-known DKT element [l 1, 121, arranged as a set of four overlapping subtriangles to avoid directional bias, and with each subtriangle having half the flexural rigidity of the actual material. This choice is made merely because the DKT element works well and is conveniently available. A good quadrilateral plate element could be used instead, probably with very similar results. c~upl~g
The coupling procedure is applied in various but similar forms in [l, 2,6]. The form used here is the very simple form described in 161,so it is presented ~ompac~y here. Along a typical element edge of length LI between nodes i and j, let s be an edgetangent coordinate and u, be the s-direction displacement, where u, varies linearly with s between nodal values u, and Use.Let the lateral displacement w along edge 0 be L,-s w=L,wi
q$
-%Wj+ ,
(0,
-
(c,),
e,>,
(7)
I
where Bi and 13~are nodal rotations in the xy plane, positive outward from the edge, and expressible in terms of fJIXi and f?,,and the angle $i between the x axis and direction ij. The average strain along edge v is taken as
+
2. (ei -
e,).
(9)
i
If standard strain transformations and eqn (9) are applied to each side of a triangular element we obtain relationships of the form c,=Ce
and
ez=SB,
(10)
where ss is a column matrix of three edge-parallel strains, the Cartesian strains are
C contains sines and cosines of edge orientations (pi , 8 contains six in-plane nodal rotation DOF, and S contains the information of eqn (9) for each edge along with sines and cosines of &. The expression L = C’S8 defines a constant strain state over the triangle. Terms of C’S are entered into the membrane strain-displacement matrix and make it possible for inextensional bending states to exist, thus removing a source of membrane locking. Unfortunately the foregoing procedure for a triangle is not directly transferrable to a quadrilateral. A least squares procedure that uses four edges to determine the three strains in L fails for a rectangle because it provides no info~ation about shear strain yXy.This defect can be cured by adding diagonals of the quadrilateral to the least squares procedure, but then a cylindrical shell element does not respond properly to moment vectors parallel to the axis of the cylinder. Both problems are solved by using the scheme of four overlapping triangles for all three strains in c = C’S& as described by yXyalone in eqn (3). Here the weighting by triangle areas A, is not necessary but provides a definite improvement in performance. The foregoing subtriangle procedure requires elevations hi at midpoints of both diagonals of the quadrilateral. They are computed as follows. Let h, , h2, h,, and hg be the hi at midpoints of sides 12, 23, 34, and 41. ~esumably these hi are available as input data. Let h, be hi at the midpoint of diagonal 24
where & and 5. A similar mid-diagonal are calculated where the latter integral contains the 6, expression of Mareuerre shallow shell theorvI. 1131. a The edee
= y I
1(6)
Bending stiffness
membrane-bending
elevation z may be taken as a parabola between nodes i and j with a mid-edge elevation hi relative to line g (Fig. 2). Thus eqn (8) yields
3
which operates on DOF [u, v, 0,, u2 u, 6,, u3 vj e,,]‘. The factor of 10e6 is rather arbitrary. A factor as large as lop2 changes element behavior very little [8]. The value 10m6 has almost no effect other than suppression of the mechanism.
551
qs are the 4 and tt coordinates of point calculation is used for h6, the other elevation, Isoparametric coordinates directly, e.g. CS=
(x2 + x4 w (XI + x2 + x3 + x4)/4 .
(12)
R. D. COOK
552
Fig. 3. View of adjacent shell elements 1 and 2 parallel to local xy planes.
The procedure of eqns (11) and (12) is approximate unless the element is rectangular or a parallelogram. A rigorous calculation of h, and h, would be much more tedious, and numerical evidence does not suggest that it would be worth the trouble.
Fig. 5. Shell roof test case. R = 25.0, L = 50.0, thickness = 0.25, E = 4.32 x lOa, v = 0, shell weight = 90.0 per unit area. A 2 x 2 mesh is shown. Vertical deflection at A is 0.3024 according to theory.
Warping correction
Thus far the element has been formulated using its projection, for which all nodes lie in the xy reference plane. If global coordinates do not place all element nodes in the same plane, element behavior is much too stiff in many problems unless modified by a warping correction. The correction used here is taken from eqns (9b), (SC), (lOc), and (11) of [7l. The correction is embodied in a matrix T such that the transformation TrIcT produces a z-direction stiffness from membrane stiffness coefficients while leaving existing membrane stiffness coefficients and bending stiffness coefficients unchanged. Matrix T is generated entirely from nodal coordinate data in the local xyz system. Membrane locking tendency
The origin of membrane locking is suggested by Fig. 3, where a nodal rotation 0, is associated with bending deformation in element 2 but has a drilling component 6, normal to the base plane of element 1; 0, is associated with membrane action in element 1. Accordingly eb is in part resisted by membrane stiffness, which is typically far greater than bending stiffness and hence causes bending deformations to be much too small. This ‘locking’ tendency vanishes as elements become coplanar. The same difficulty appears in a triangular shell element [6]. The element of [6] contains a form of constant-strain triangle and requires two stabilization matrices: one is described by eqn (6) and the other is
a rank-two matrix that operates on only the drilling DOF. The latter matrix is the sole cause of membrane locking in shells, but must be present when elements are coplanar. This matrix can be multiplied by a factor that is unity for coplanar elements but approaches zero as the curvature of the element increases. Thus it is easy to avoid membrane locking in the triangular shell element. The corresponding remedy for the quadrilateral would use /_I= 1 in eqn (4) for coplanar elements and arrange for /I to approach zero as elements become more highly curved. However, edge tractions on elements produce nodal torques corresponding to the 0,, when /I = 1, but zero torques when #I = 0. Hence adjacent elements having different values of /I would be incompatible (with the triangular element, similar adjustments of the second stabilization matrix present no such difficulty). At this time the writer does not know of a satisfactory cure, but fortunately the quadrilateral element does not have a great propensity for membrane locking. It may be noted that use of one-point integration for only the membrane terms associated with the tlli, with four-point integration for other terms, produces elements that either fail the patch test or behave badly otherwise.
1.0
.A__
OF
1.0
I I-
L
“.”
lW
OF
0.9
I.1 1.9
Diaphragm
3.1 3.9
2.1 2.9
4.1 4.9
5.1
.I
I.0 t 1.0
Fig. 4. Plane beams of Poisson ratio v = 0.3 loaded in pure bending. Drawing is not to scale.
Fig. 6. Pinched cylinder test case. R = 300.0, L = 600.0, thickness = 3.0, E = 3.0 x 106, v =0.3, F = 1.0. A 2 x 2 mesh is shown. Displacement at A is 0.18248 x 10e4 according to theory.
Four-node ‘tlat’ shell element
J
553
Fig. 8. Cantilevered strip with end-to-end twist of 90”. L = 12, b = 1.1, r = 0.0032, E = 29 x 106,v = 0.22. A 2 x 4 mesh is shown. Load F = 10e6 applied as forces 0.25 x 10m6, 0.50 x 10e6, and 0.25 x 10e6 at the three end nodes. Displacement at C is 5256 x 10e6 according to theory.
Free
+F
X Fig. 7. Hemispherical shell test base. One quadrant with a 4 x 4 mesh is shown. Radius = 10.0, thickness = 0.04, E = 6.825 x lo’, v = 0.3, F = 1.0 for the quadrant modeled. Displacement at A is 0.0924 according to theory. NUMERICAL RRSULTS
Membrane behavior The two beam problems of Fig. 4 are standard test cases [8,14]. The arrangement of trapezoidal elements is a test case in which four-node quadrilaterals often display shear locking. Expressed as the ratio of computed tip deflection w to exact tip deflection w,,, the results are rectangular
mesh, fl = 0:
w/w,, = 1.000
rectangular
mesh, /I = 1:
w/w,, = 1BOO
distorted mesh, fl = 0:
w/w,, = 0.051
distorted mesh, /I = 1:
w/w, = 0.897.
Clearly, for rectangular elements j does not matter, while badly shaped elements suffer unless fl = 1. With j? = 1, the present membrane formulation appears to have almost the same behavior as the membrane element suggested in [14]. Shell test cases The standard ‘obstacle course’ test cases [15] are depicted in Figs 5-8. Data for these problems are stated in the captions. In Figs 5 and 6, nodes at a diaphragm support are assumed free to rotate about a normal to the diaphragm, in accord with the simple support condition recommended in [16]. In Fig. 7 a small hole of radius 0.175 is placed at the top, where all DOF are restrained. All meshes ate uniform
Table 1. Shell test case results reported as the ratio of computed deflection to accepted theoretical deflection Element formulation
Mesh
Flat /3=0
Flat /?=I
5.085 1.548 1.084 1.014
CUNtd
&Ned
j=o
/I=1
4.941 1.541 1.085 1.014
2.951 1.109 1.010 0.995
2.471 1.080 1.009 0.996
0.012 0.098 0.630 0.938
0.012 0.097 0.627 0.931
2.125 1.150 0.941 1.016
0.015 0.178 0.775 0.990
Hemisphere, Fig. 7 1x1 2x2 4x4 8x8
0.117 0.852 1.002 1.001
0.056 0.228 0.444 0.948
1.533 0.654 0.586 0.971
0.057 0.168 0.606 0.962
Twisted strip, Fig. 8 1x1 2x2 2x4 2x8
0.916 0.859 0.978 0.999
0.122 0.608 0.969 0.993
Shell rooJ Fig. 5 1x1 2x2 4x4 8x8 Pinched cylinder, Fig. 6 1X1
2x2 4x4 8x8
t Element sides are straight in a rectangular mesh.
r;
1;
554
R. D.
COOK
Table 2. Shell test case results using the ‘curved softened’ TRIANGULAR Mesh
Shell roof
IX1
moo
2x2 4x4 8x8
1.406 0.880 0.939
element of [6]
Pinched cylinder
Hemisphere
Twisted strip
0.121 0.959 0.974 0.999
0.213 0.873 0.978 1.007
0.775 0.962t 0.992$
t 2 x 4 mesh. $2 x 8 mesh.
except those used for the hemisphere, where smaller elements appear at the top. In Fig. 8, note that f = 0.0032 rather than the less demanding case of t = 0.32 that also appears in [15]. Results appear in Table 1. The ‘flat’ formulation is obtained by omitting the membrane-bending coupling device, i.e. by setting hi = 0 in eqn (9) for all midpoints. We see from Table 1 that /I = 1 provides a stiffer formulation that /_l= 0, but the tendency toward membrane locking provided by /? = 1 is small. For the shell roof and pinched cylinder cases, the best element formulation appears to be curved with /I = 1. For the pinched cylinder the curved element with fl = 0 appears to converge faster; however in this problem the zone of greatest deformation is very localized around the load points, so the apparent good convergence properties of this element may instead be a symptom of a too-flexible formulation. In the hemisphere test case there is considerable bending and rigid body motion of elements but little membrane strain. For this case the best element is flat and without drilling DOF. Similar behavior is observed with triangular elements [2]: flat constantstrain triangles without drilling DOF are very accurate in the hemisphere problem. The twisted strip problem clearly shows the need for the warping correction if elements are thin. Without the warping correction, deflection ratios in the last line of Table 1 become 0.0001 and 0.0003. The regular meshes in Figs 5-7 do not result in warped elements. With the pinched cylinder of Fig. 6, numerical tests in which nodes were relocated so that elements became warped showed that the warping correction sometimes stiffened and sometimes softened the response but always made the computed result more nearly exact. In comparison with other available four-node quadrilaterals [17,18] the present curved element with /I = 1 appears to be more accurate in the shell roof and pinched cylinder problems but less accurate in the hemisphere problem. The twisted strip problem is not used in [17, 181. For comparison with results given by a triangular element [6] we refer to Table 2. Triangular elements are produced by dividing each quadrilateral along its shorter diagonal. A two-element patch of triangular elements has the same number of DOF as one quadrilateral element, so that the corresponding meshes in Tables 1 and 2 use the same number of DOF.
CONCLUDINGREMARKS This paper has considered shell elements obtained by the very simple process of combining standard membrane and bending formulations with a device for membrane-bending coupling and a device for inclusion of warping effects. The process leads to reasonably good four-node quadrilaterals and threenode triangles. Quadrilaterals have better membrane behavior than triangles and are free of the directional bias that usually accompanies triangles, e.g. when triangles form a cross-hatch pattern rather than a herringbone or union-jack pattern. Also, quadrilaterals do not have the ‘free comer’ problem present in triangles, in which a curved element may display large and erroneous displacements if it has a loaded node that is not shared by another element. If used without internal DOF, static condensation is avoided, while a quadrilateral built of four non-overlapping triangles would have six internal DOF to be condensed before assembly of elements and would require additional data to properly locate the internal node on the shell surface. On the other hand, two triangles are needed to replace one quadrilateral, so there is more interelement membrane-bending coupling with triangles than with quadrilaterals for a given number of nodes, the smaller span of the triangle reduces the tendency toward membrane locking. If the tendency for membrane locking in the quadrilateral can be reduced it may prove to be a very effective element.
REFERENCES 1. H. Stolarski, T. Relytschko, N. Carpenter and J. M. Kennedy, A simple triangular curved shell element. Engng Comput. 1, 210-218 (1984). 2. N. Carpenter, H. Stolarski and T. Belytschko, Improvements in 3-node triangular shell elements. Int. J. Numer. Meth. Engng 23, 1643-1667 (1986). 3. D. J. Allman, A compatible triangular element including vertex rotations for plane elasticity analysis. Comput. Strut. 19, l-8 (1984). 4. N. Carpenter, H. Stolarski and T. Relytschko, A flat triangular shell element with improved membrane interpolation. Comm. Appl. Numer. Meth. 1, 161-168 (1985).
5. J. Fish and T. Belytschko, Stabilized rapidly convergent 18-degrees-of-freedom flat triangular shell element. Inr. J. Numer. Meth. Engng 33, 149-162 (1992). 6. R. D. Cook, Further development of a three-node triangular shell element. Inr. J. Numer. Meth. Engng. 36, 1413-1425 (1993).
Four-node ‘flat’ shell element 7. B. P. Naganarayana and G. Prathap, Force and moment corrections for the warped four-node quadrilateral plane shell element. Comput. Struct. 33, 1107-1115 (1989). 8. R. H..MacNeal and R. L. Harder, A refined four-noded membrane element with rotational degrees of freedom. Comput. Struct. 28, 75-84 (1988). 9. R. D. Cook, D. S. Malkus and M. E. Plesha, Concepts and Applications of Finite Element Analysis, 3rd Edn. John Wiley, New York (1989). 10. R. D. Cook, Avoidance of parasitic shear in plane element. AXE J. Struct. Div. 101, 1239-1253 (1975). Il. J. A. Stricklin, W. E. Haisler, P. R. Tisdale and R. Gunderson, A rapidly converging triangular plate element. AIAA Jnl7, 180-181 (1969). 12. C. Jeyachandrabose, J. Kirkhope and C. R. Babu, An alternative formulation for the DKT plate-bending element. ht. J. Numer. Meth. Engng 21, 1289-1293 (1985). 13. K. Marguerre, Zur theorie der gekriimmten plate grosser formlnderung. Proc. Fifih Int. Cong. for Appl.
14.
15.
16.
17.
18.
555
Me& (Edited by J. P. Den Hartog and H. Peters), pp. 93-101. John Wiley, New York (1939). M. Iura and S. N. Atluri, Formulation of a membrane finite element with drilling degrees of freedom. Comput. Mech. 9, 417-428 (1992): _ T. Belytschko, B. L. Wong and H. Stolarski, Assumed strain stabilization procedure for the 9-node Lagrange shell element. Int. J. Numer. Meth. Engng 28, 385414 (1989). 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). W. K. Liu, E. S. Law, D. Lam and T. Belytschko, Resultant-stress degenerated-shell element. Comput. Meth. Appl. Mech. Engng 55, 259-300 (1986). K. Y. Sxe and C. L. Chow, A mixed formulation of a four-node Mindlin shell/plate with interpolated covariant transverse shear strains. Comput. Struct. 40, 775-784 (1991).