A comparison of rotation-free triangular shell elements for unstructured meshes

Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015 www.elsevier.com/locate/cma

A comparison of rotation-free triangular shell elements for unstructured meshes Mattias Ga¨rdsback, Gunnar Tibert

*

KTH Mechanics, Royal Institute of Technology, Osquars backe 18, SE-100 44 Stockholm, Sweden Received 25 January 2007; received in revised form 25 June 2007; accepted 28 June 2007 Available online 20 July 2007

Abstract Many engineering applications require accurate and rapidly computed thin-shell elements. Rotation-free (RF) shell elements include the bending behaviour of thin shells without introducing any additional degrees of freedom compared to a membrane element. Instead, constant curvatures are approximated from the out-of-plane displacements of a patch of usually four triangular elements. A consequence of this is that the accuracy for irregular meshes has been unsatisfactory. The aim of this study is to find an RF shell element which is accurate also for unstructured meshes. The main difference between existing elements is whether they assume two-dimensional constant curvatures over the patch or use superposition of one-dimensional constant curvatures for the three pairs of triangles. The first assumption fulfils constant curvatures for a Kirchhoff plate exactly, whereas the second and most common assumption only approximates constant curvatures. The first assumption is significantly more resistant to element shape distortions, whereas the second assumption is slightly faster to compute and more appropriate on boundaries where one or more elements are missing or several neighbouring elements share a side. The combination is significantly more accurate for irregular meshes than other comparable RF elements for linear benchmark tests.  2007 Elsevier B.V. All rights reserved. Keywords: Shell element; Rotation-free; Constant curvatures; Unstructured mesh

1. Introduction Thin-film structures in space and sheet metal forming are examples of important engineering applications that require simple and rapidly computed thin-shell elements. Rotation-free (RF) shell elements is a family of thin-shell elements which includes the plate bending behaviour without introducing any rotational degrees of freedom (d.o.f.). Thus, it should be possible to use less d.o.f. for many thinshell problems and avoid troubles with ill-conditioned stiffness matrices caused by the introduction of rotational d.o.f. However, it has been argued in [1] and shown in e.g. [2] that the accuracy of RF elements is unsatisfactory for unstructured meshes. Therefore, the aim of this article *

Corresponding author. Tel.: +46 8 7907961; fax: +46 8 7969850. E-mail address: [email protected] (G. Tibert). URL: www.mech.kth.se (G. Tibert).

0045-7825/$ - see front matter  2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2007.06.017

is to derive an RF shell element that is accurate also for unstructured meshes. The RF shell elements are formed as a superposition of a simple membrane element and an RF bending element. Kirchhoff theory is used to approximate constant curvatures from the out-of-plane displacements of a patch of elements, consisting of the main element and its three neighbours in the case of triangular elements. For the plate part, many RF shell elements approximate the constant curvatures over the whole patch [3], over the whole patch using numerical integration [4,5] or over the edges [6–8]. We will show that the first approach is advantageous inside an unstructured mesh and that the last is advantageous on the boundary. The novelties of this study are: (i) the cause of the sensitivity of RF shell elements to shape distortions are for the first time explained and the Phaal–Calladine element, insensitive to such distortions, is re-introduced; (ii) the

5002

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

effects of different RF approximations of the boundary conditions are evaluated. 1.1. Prior work on rotation-free shell elements Barnes [9] presents the theory, but no numerical results, of an element where the average rotations of the elements sharing an edge are used to approximate constant curvatures normal to that edge. The superposition of the contributions to the generalised curvatures, from the curvatures of the three edges, are used to form the stiffness matrix. Chan and Davies [10] use the transversal displacements of an element to approximate the rotation angle of the element, and assume a mirror image element on the other side of the edge. These elements are investigated numerically by Hampshire et al. [11], who also extend the ideas by Chan and Davies by using the adjacent element instead of assuming a mirror image element. Phaal and Calladine [12] use a complete quadratic polynomial in two dimensions to assume constant curvatures over the whole patch of elements. The three unknown generalised curvatures are determined using the three approximations of the hinge angles. Their plate bending element is refined to cope with initially non-coplanar element patches, and combined with a constant strain triangle (CST) [13] to form a shell element [3]. It has been argued [6,14] that this formulation is inconvenient, mainly because of the fictitious nodes that must be introduced outside the boundary. On˜ate and Cervera [6] use integration by parts of the curvatures to derive their plate bending element. This leads to a formulation where the integral of the generalised curvatures over the element area is replaced with a sum of the integrals of the normal to each edge multiplied with the deflection gradient over the same edge. Combined with the CST element, a refined non-coplanar version of the element in [6] has been formulated as a shell element by On˜ate and Za´rate [15], which has been extended to be valid for large strain plasticity, and kinked and branched surfaces by Flores and On˜ate [16,17]. They have also suggested to use overlapping isoparametric linear strain triangles (LST) for better membrane behaviour [5,18]. Sabourin and Brunet [19] use the stiffness matrix proposed by Barnes [9], the relation between the rotation angles and displacements by Chan and Davies [10], and boundary conditions corresponding to the ones by On˜ate and Cervera [6] in their plate bending element, which is also extended to shells [7]. Guo et al. [8] use the DKT6 plate element developed by Batoz et al. [20], which is identical to Morley’s constant-moment plate bending element [21]. The rotational d.o.f. are approximated using the rotation angles of the two elements sharing each edge. Other, completely different RF formulations do exist. Cirak et al. [22] use subdivision surfaces, which are extended to large deflections [23]. However, subdivision surfaces is not a standard FE and consequently their formulation is not directly compatible with standard FE tools. Rio et al. [4] develop a shell element which takes into

account both large displacements and coupling effects between membrane and bending behaviour. This element uses isoparametric interpolation in an attempt to model constant curvatures normal to each edge, which leads to problems in finding suitable integration points. Hauptmann and Schweizerhof [24] have developed an RF element degenerated from solid elements, where the displacements of the top and bottom surfaces are required. RF elements have also been developed for quadrilateral elements [14]. However, this study is restricted to triangular elements since triangles are easier to fit to an arbitrary geometry. 1.2. Merits and demerits of rotation-free shell elements RF shell elements should only be used for thin shells. As with all finite-element formulations they have merits and demerits. The merits of RF shell elements are: • The total number of d.o.f. necessary to obtain a certain accuracy may be less than half for an RF element compared to the corresponding elements that use rotational d.o.f., i.e. Morley’s constant-moment plate-bending triangle and DKT, for many benchmark examples [2,17]. • The thinner the shell becomes, the greater the ratio between the membrane and flexural stiffness becomes. For an extremely thin shell, the stiffness matrix for an RF element will be identical to that of the membrane element, and better conditioned than that of shell elements that include rotational d.o.f. [25]. • Rotations are related to out-of-plane displacements by first derivatives. Since these d.o.f. are not fully independent, the stiffness matrix may be poorly conditioned, especially for fine meshes, if rotations are used [12]. • Rotations and higher-order variables are not commutative in three dimensions. Therefore, shell elements employing these d.o.f. must track strains and rigid-body modes in the case of large-deflection analysis [12]. • RF elements are well-suited for parallel computing since many similar and simple calculations must be performed [11]. • Rotational and higher-order d.o.f. are sometimes difficult to apply when symmetric boundary conditions are used, and when elements are joined at corners and intersections [12]. In addition, contact conditions are simplified when rotational d.o.f. are removed [2]. and the demerits are: • The accuracy of a thin-shell is often almost the same as the accuracy of the membrane part, and simple membrane elements, such as the CST or overlapping LST [5,17,18], are used together with an RF bending part. • The plate bending parts of most RF shell element formulations are sensitive to element shape distortion.

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

• The bandwidth of the total stiffness matrix is increased compared to the membrane element. Compared to many classical shell elements the bandwith is not decreased at the same rate as the number of d.o.f. [4]. • Direct combination of RF shell elements and other finite elements is not always straightforward, and requires special multipoint constraints [17]. Despite these demerits, the present study will show that the element shape distortion sensitivity can be removed from the list of demerits above. 2. Theory of rotation-free shell elements This section presents various RF shell element formulations in a unified way that enables fair and accurate comparison. It is shown in [26] that most elements in this study yield similar results, but here the differences and similarities will be highlighted. 2.1. Formulation of the stiffness matrix The RF elements belong to the family of flat shell elements formed by the superposition of a membrane element and a plate bending element. The strain energy, U, is expressed as the superposition of the energies due to the membrane and bending strains: Z Z 1 T 1 T U ¼ Um þ Ub ¼ e rm dV þ e rb dV : ð1Þ 2 m 2 b For linearly elastic materials, the relation between stresses, r, and strains, e, is given by Hooke’s law: r ¼ Ee;

ð2Þ

where the constitutive matrix for thin shells, using the assumption of plane stress, is 2 3 1 m 0 E 6 7 E¼ ð3Þ 4 m 1 0 5; 1  m2 1m 0 0 2

where w is the displacement in the z-direction and j the generalised curvatures. In FE calculations, em and j are interpolated from the element d.o.f. as em ¼ Bm d;

ð6Þ

j ¼ Bb d;

ð7Þ

where d are the element d.o.f. The plate bending moments are expressed as 0 1 mx B C ð8Þ mb ¼ @ my A ¼ Dj; mxy where D = t3E/12. The element stiffness matrix Ke then becomes Z Z ð9Þ Ke ¼ Km þ Kb ¼ tBTm EBm dA þ BTb DBb dA: The formulation of the membrane stiffness will be briefly discussed in Section 2.8.2, but is not within the scope of this study. Constant curvatures are assumed for all the investigated RF elements since interpolation polynomials of a higher order than two requires either more d.o.f. per element or more elements in the element patch. The contribution from the bending part to the element tangent stiffness matrix becomes Kb ¼ ABTb DBb :

ð10Þ

Thus, the only difference regarding the bending part of the stiffness matrix of the various RF shell element formulations is in the curvature approximation Bb. 2.2. Triangular shape functions Triangular shape functions, Fig. 1, are used by e.g. Morley [21] and On˜ate et al. [5,6]. The transversal displacement, w, can be approximated from the nodal transversal displacements, wi, as

c2

where E is the Young’s modulus and m is the Poisson’s ratio. The linear membrane strains, em, are 0 1 0 ou 1 ex ox B e C B ov C em ¼ @ y A ¼ @ oy A; ð4Þ ou ov cxy þ ox oy where u and v are the in-plane displacements. Thin plates are described by Kirchhoff plate theory, with bending strains 0 o2 w 1 0 1  ox2 jx B o2 w C B C C eb ¼ zj ¼ z@ jy A ¼ zB ð5Þ @  oy 2 A; ow jxy 2 oxoy

5003

w3

n2 b2

γ2

l2 s2

s1

L2

γ1 L1

w1 γ3 b3

l3

n3

y

L3

c3

n1 l1

s3

b1

w2 c1

x Fig. 1. Definition of triangular shape functions.

5004

w¼

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015 3 X

ð11Þ

wi Li ;

i¼1

where the linear shape functions Li for a triangle are given by Li ¼

ai þ bi x þ c i y ; 2A

ð12Þ

where a1 ¼ x2 y 3  x3 y 2 , b1 ¼ y 2  y 3 and c1 = x3  x2. Constants a2 , b2 , c2 , a3 , b3 and c3 are obtained by cyclic permutation of the subscripts. By repetitive use of the Pythagorean Theorem, the area of the triangle is obtained as 2A ¼ c3 b2  c2 b3 ¼ c1 b3  c3 b1 ¼ c2 b1  c1 b2

The rotation angle ai is then approximated as the directional derivative in the opposite direction to the normal ni ai ¼ $w  ðni Þ: For side 1, the rotation angle becomes a1 ¼

w1 w2 w3  cos b3  cos b2 : h1 h2 h3

ð19Þ

The remaining rotations are obtained analogously, with distances and angles defined in Fig. 2. Thus, the rotation angles are described by

ð13Þ a ¼ Vw;

and the length of side i is qffiffiffiffiffiffiffiffiffiffiffiffiffiffi li ¼ b2i þ c2i :

ð14Þ

ð20Þ

where aT = (a1. . .a6), wT = (w1. . .w6) and 2

The normal, in the plane, to side i is     cos ci 1 bi ni ¼ ¼ ; li c i sin ci

ð15Þ

where ci is the angle between the positive x-axis and the normal ni. It follows that the derivative with respect to the direction of ni satisfies o o o cos ci þ sin ci : ¼ oni ox oy

1 h1

6 cos b 3 6 h1 6 6 cos b 2 6 h1 6 V¼6 6 0 6 6 cos / 3 6 4 s1  cosq c2

0

0

1 h2

 cosh3b1

0

0

 cosh2b1

1 h3

0

0

 cosr2w3

 cosr3w3

1 h4

0

0

 coss3/1

0

1 h5

0

0

0

2

3

7 07 7 7 07 7 7: 07 7 7 07 5

ð21Þ

1 h6

ψ3

β2

β3

r2

r3

h4

2 γ2

2

H

3

h1

h2 h6 h3

β1

γ1

6

ð22Þ

For i = 1, 2 and 3, ai is computed with the shape functions of the main triangle and for i = 4, 5 and 6, ai is defined with the shape functions of triangle j, where j = i  3. The three nodes of triangle j are numbered as in Fig. 3. Repetitive use of Eqs. (12) and (15) in Eq. (22) results in

4

ψ2

φ1

0

 3  X ow oLj oLj ai   þ sin ci ¼ cos ci : oni ox oy j¼1

4

φ3

 cosh3b2

Alternatively, the triangular shape functions can be used to obtain identical expressions of the rotation angles [21]. Let Eq. (16) operate on Eq. (11) to obtain:

The rotation angles of two adjacent elements, ai and ai+3, can approximately describe the bending behaviour over the common edge i. The rotation angles, a, are just rigid-body rotations and can be determined in different ways, resulting in equal expressions. Phaal and Calladine [12], Sabourin and Brunet [19] and Guo et al. [8] all end up with the same matrix, which can be derived by noting that the gradient of the flat triangle is [14]   w1 w2 w3 $w ¼  n1 þ n2 þ n3 : ð17Þ h1 h2 h3

3

 cosh2b3

 cosq c1

1

ð16Þ

2.3. Rotation angles and hinge angles

5

ð18Þ

5

h5 s 3

q1

s1

1

6 1

Fig. 2. Definition of angles and distances.

q2

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

5 2

3

z*

4

3 3

1

5005

ns

3

2

A

A

1

A 1

1 1 2

2 1

2

2

θ1

w1* α1

3

A

α4

h1

w2* = w3* =0

w4* n

h4

3

Fig. 4. Coordinate system used to determine the one-dimensional constant curvature (dashed line) for a pair of triangles.

6 Fig. 3. Node numbering in the patch and in each triangle.

2

l1 2A

6 c c þb b 6 12 12 6 2Al2 6 c c þb b 6 13 13 6 2Al3 6 V¼6 6 0 6 6 c2 c2 þb2 b2 6 1321 3 6 2A l2 4 33 33

0

0

l2 2A c2 c3 þb2 b3 2Al3

c3 c1 þb3 b1 2Al1 c3 c2 þb1 b2 2Al2 l3 2A

0

0

0

0

c12 c13 þb12 b13 2A1 l1

c11 c13 þb11 b13 2A1 l1

l1 2A1

0

0

c22 c23 þb22 b23 2A2 l2

0

l2 2A2

0

0

0

c2 c1 þb2 b1 2Al1

c2 c3 þb2 b3

c31 c33 þb31 b33

2A3 l3

2A3 l3

0

where n is the coordinate in the normal direction and w* is the part of the deflection, w, that contributes to the curvatures. The curvature for pair 1 is obtained from the second derivative of Eq. (26):

3

7 0 7 7 7 0 7 7 7 7; 0 7 7 7 0 7 7 5

ð23Þ

l3 2A3

where Aj is the area of triangle j, bji and cji are defined by the shape functions of triangle j. The change in the hinge angle around edge i is obtained as the sum of the two rotation angles, ai and ai+3: 0 1 a1 þ a4 B C b ¼ @ a2 þ a5 A ¼ Ww; ð24Þ a3 þ a6 where row i of W is the sum of rows i and i + 3 of V. 2.4. Polynomials to describe constant curvatures A complete quadratic polynomial in two dimensions describes the out-of-plane displacement, w, of a flat patch with six nodes [12]: w ¼ c1 þ c2 x þ c3 y þ c4 x2 þ c5 xy þ c6 y 2 :

ð25Þ

The first three terms describe the rigid-body motion, whereas the last three terms represent two-dimensional constant curvatures for a Kirchhoff plate by Eq. (5) [12]. However, all six terms are required in order to have an expression that is invariant with respect to the coordinate system. The one-dimensional constant curvature for the pair of triangles sharing edge 1, Fig. 4, can be described by a second order polynomial of the relative vertical displacement w* [8]: nðn  h4 Þ  nðh1 þ nÞ  w þ w; w ðnÞ ¼ h1 ðh1 þ h4 Þ 1 h4 ðh1 þ h4 Þ 4

ð26Þ

jn;1 ¼

d2 w 2 ¼ ða1 þ a4 Þ; h1 þ h4 dn2

ð27Þ

where ai ¼ wi =hi was used. Analogous expressions can be obtained for sides 2 and 3. Both polynomials, Eqs. (25) and (26), assume initially planar patches. However, using the assumption that equal changes of the hinge angles cause equal strains extends their applicability to curved shells. In the case of a non-flat shell for the first polynomial, the patch is first flattened to the xy-plane by a coordinate transformation, and then back-transformed to its original coordinates in the end. For the second polynomial, the angles ai in Eq. (27) are computed from their initial positions. There exists no exact relation between the two-dimensional generalised curvatures j and the one-dimensional curvatures for a pair of triangles jn. Instead, Sabourin and Brunet [19] and Guo et al. [8], assume superposition of the contributions. For edge i o o o ¼ cos ci þ sin ci ; ox oni osi o o o ¼ sin ci þ cos ci ; oy oni osi ow ¼ 0; osi

ð28aÞ ð28bÞ ð28cÞ

where si is the coordinate along edge i. The superposition becomes j ¼ Rjn ; where 2 6 6 R¼6 6 4

ð29Þ

b21

b22

b23

l21

l22

l23

c21

c22

c23

l21

l22

l23

2 b1l2c1 1

2 b2l2c2 2

3

2 b3l2c3 3

7 7 7: 7 5

ð30Þ

5006

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

2.5. Morley’s constant-moment plate bending element Many of the RF shell elements resemble Morley’s constant-moment plate bending element [21]. The difference is that the out-of-plane rotations are approximated using the out-of-plane displacements of the adjacent nodes. Therefore, a study of RF elements should start with this element. The derivation that follows is adjusted to resemble RF elements. The generalised curvatures are related to the curvatures along the normal directions according to Eq. (30). The onedimensional constant curvature for pair 1, Fig. 4, can be obtained from Eq. (27). However, a4 is not available for Morley’s element since it is not an RF element. Instead, the rotation of the shell normal at edge 1, h1, is obtained from the first derivative of Eq. (26) evaluated at n = 0 [8]: h1 

dw ð0Þ h4 h1 ¼ a1  a4 : dn h1 þ h4 h1 þ h4

2 ða1  h1 Þ: h1

Analogous expressions can be obtained for sides 2 and 3. The relationship between the one-dimensional curvatures and the original d.o.f., q, is written as ð33Þ

where qT ¼ ðw1 ; w2 ; w3 ; h1 ; h2 ; h3 Þ, 2 l1 3 0 0 A 6 7 HM ¼ 4 0 lA2 0 5 0 and

2

0

0

A  AþA 2

0

0

A  AþA 3

3 7 7 7: 5

ð38Þ

The d.o.f. in Eq. (36), q can then be related to the transversal d.o.f., w, by q ¼ Mw; where " M¼

I

ð39Þ

0

Q V

# ð40Þ

;

where I is the 3 · 3 identity matrix and 0 is a 3 · 3 zero matrix. Finally, the curvatures are given by Bb ¼ RHM TM:

ð32Þ

jn ¼ HM Tq;

0

ð41Þ

ð31Þ

Eqs. (27) and (31) are used to obtain the following expression for the one-dimensional constant curvature: jn;1 ¼

where hT ¼ ðh1 ; h2 ; h3 Þ and 2 1 A A 0 0  AþA 1 1 6 AþA 2 6 A Q¼6 0 0 0 AþA2 4 3 A 0 0 0 AþA3

ð34Þ

l3 A

2.6.2. Element by Sabourin and Brunet Sabourin and Brunet [19] also approximate the generalised curvatures as a superposition of the one-dimensional curvatures for the three pairs of triangles, Eq. (29). As Eq. (27) describes constant curvature for pair 1, the relationship between the one-dimensional curvatures and the rotation angles becomes jn ¼ HS b; where

2

l1 AþA1

6 6 HS ¼ 6 0 4 0

ð42Þ

0 l2 AþA2

0

0

3

7 0 7 7; 5

ð43Þ

l3 AþA3

with b from Eq. (24). The generalised curvatures for the plate are thus described by

v11 6 T ¼ 4 v21

v12 v22

v13 v23

1 0

0 1

v31

v32

v33

0

0

3 0 7 0 5; 1

Bb ¼ RHS W:

ð44Þ

ð35Þ

This element is identical to the element by Guo et al. [8], but more straightforward to derive.

where vij are elements in V. The generalised curvatures for Morley’s element become

2.6.3. Element by Phaal and Calladine Phaal and Calladine [12] use the complete quadratic polynomial in two dimensions, Eq. (25). A relationship between the unknown transversal displacements, w, and the generalised curvatures, j, is obtained if the x- and ycoordinates of the six nodes of the flattened patch are inserted into the quadratic polynomial:

j ¼ RHM Tq:

ð36Þ

2.6. Rotation-free plate bending elements 2.6.1. Element by Guo et al. Guo et al. [8] use DKT6 by Batoz [20] as the origin of the plate part of their element. The DKT6 element is identical to Morley’s constant-moment plate bending element [21]. From Eq. (31), a relationship between the rotations and the deflection angles can be derived h ¼ Qa;

ð37Þ

w ¼ r1 þ r2 þ r3 þ Cj;

ð45Þ

where rT1 ¼ c1 ð1 . . . 1Þ, rT2 ¼ c2 ðx1    x6 Þ, rT3 ¼ c3 ðy 1    y 6 Þ and 2 2 3 x1 x22 x23 x24 x25 x26 16 7 ð46Þ CT ¼  4 y 21 y 22 y 23 y 24 y 25 y 26 5: 2 x1 y 1 x2 y 2 x3 y 3 x4 y 4 x5 y 5 x6 y 6

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

The hinge angles are obtained from Eq. (24). Note that W(r1 + r2 + r3) = 0 since the first three coefficients in Eq. (25) correspond to rigid-body motion. Then, a relation between the hinge angles and the curvatures is obtained by inserting Eq. (45) into Eq. (24) b ¼ WCj:

ð47Þ

Taking the inverse and using Eq. (24) again gives the following expression for the curvatures: 1

Bb ¼ ðWCÞ W:

ð48Þ

2.6.4. Element by On˜ate and Cervera The shortest formulation is the one by On˜ate and Cervera [6]. First, the generalised curvatures over an element are expressed as the mean of the curvatures over the element 2 3 2  oxo 2 Z Z 6 7 1 6  o22 7w dx dy; j¼ ð49Þ oy 5 4 A A 2 o 2 oxoy where A is the area of the element. Integration by parts yields Z 1 j¼ T$w ds; ð50Þ A C where C is the boundary of the element, 2 3 nx 0 6 7 T ¼ 4 0 ny 5 ny nx and $w ¼

ow ox ow oy

!



nx ¼ ny

ny nx



ow on ow os

ð51Þ

Bb ¼ RHM VO ;

ð54Þ

where VO contains the three last rows of V in Eq. (21) or (23). 2.7. Boundary conditions As the RF elements also use adjacent nodes in the computation of the curvatures, handling of boundary conditions (b.c.), where at least one node is missing, becomes slightly more complicated than for elements with rotational d.o.f. The b.c. for a simply supported edge uses the free edge b.c. in combination with pinning of the boundary nodes. A clamped edge is obtained by pinning the nodes at the edge and applying the symmetry b.c. Therefore, it is only necessary to describe the free and symmetry b.c. Phaal and Calladine use a fictitious mirror node to replace the missing node outside the boundary. For the symmetry b.c., the contributions to the curvatures from the mirror node are added to the corresponding node in the main element. For the free edge b.c., the contributions are subtracted. At a free edge, the moments mxy and either of mx or my become equal to zero for the whole element and not just along the boundary. On˜ate and Cervera [6], Sabourin and Brunet [19] and Guo et al. [8] apply their b.c., here denoted the common b.c., directly on the rotation angles a. For a free boundary no bending occurs, so the rotation angle of the missing element is opposite to that of the main element: aiþ3 ¼ ai :

ð55Þ

For a symmetric or clamped boundary, the rotation angle of the missing element is equal to that of the main element:

!

aiþ3 ¼ ai : ð52Þ

:

5007

Integration by parts implies in this case that the twodimensional generalised curvatures are approximated as a superposition of the contributions from the three edges. The contribution to the curvatures from ow/os is null for a flat triangle. Therefore, 0 ow 1 0 1 n2x;i on 3 X li B B ow1 C C ow C ð53Þ ¼ RHM B j¼ @ n2y;i A @ on2 A: oni A i¼1 ow 2nx;i ny;i on 3

To compute ow/oni at the edge, On˜ate and Cervera [6] use the average from the two elements sharing the edge. In theory, a more accurate alternative would be to use Eq. (26) to ensure correct one-dimensional constant curvatures for each pair of triangles, similarly as in [27]. In [5] an isoparametric approach with quadratic shape functions is used both for the bending and membrane part. Finally, using the average of the normal gradients at the boundary, and observing that the contributions from the central triangle then cancels out [6], gives

ð56Þ

The use of Eq. (55) gives zero moment around the free edge, but not for the whole element since the other hinge angles contribute to all moments mb. Contrary, the free edge b.c. by Phaal and Calladine always gives mxy and either mx or my equal to zero for the whole element, which produces a more flexible solution. It can be argued which is more correct. Flores and On˜ate [5,16,27,18] propose to set the moment normal to the edge to zero by adding a term that sets the corresponding curvatures to zero also for the formulation by On˜ate and Cervera. This boundary condition is hereafter denoted ‘‘Flores–On˜ate’’. Superposition of constant curvatures over the edges is advantageous for the symmetry b.c. since a more balanced expression for the curvatures, and consequently the stiffness matrix, is obtained. This can be visualised for a symmetric mesh as in Fig. 5a, and this is useful to explain the mesh sensitivity of the formulation by Phaal and Calladine. The curvature normal to the edge becomes, when the symmetry b.c. is applied, for Phaal and Calladine: o2 w w3  w2 ¼2 ; 2 oy h2

ð57Þ

5008

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

2 3

5

1

4

where the node numbers in Fig. 5b were used. For the common b.c. the expression becomes: o2 w w3 þ w5  w1  w2 ¼ : oy 2 h2

ð58Þ

As a consequence, the right corner is stiffer than the left corner of the mesh in Fig. 5a when the b.c. by Phaal and Calladine is applied, whereas the corners are equally stiff when the common b.c. is applied. 2.8. Rotation-free shell formulation 2.8.1. Bending part For both flat plates and initially curved shells, the rotation angles a of the element patch are calculated relative to their initial positions. Hence, the curvatures for a shell element can be calculated similar to those of a plate element. First, the patch of elements with six transversal d.o.f. are split into four triangles, each with three transversal d.o.f., Fig. 6. Local coordinate systems are introduced for each of the four triangles and the rotation angles of the adjacent triangles are calculated in the local frames. For a shell, Eq. (20) is therefore changed to a ¼ Vs ws ;

ð59Þ

where wTs ¼ ðw21 ; w1 ; w31 ; w32 ; w2 ; w12 ; w13 ; w3 ; w23 ; w14 ; w25 ; w36 Þ and

n2

y2

w32

w3

z 2 x2

z1

y1 x1

n w12 w1 w13

y x

z

w41

n1

w31

w2

v21 0 v31 0

w21

0

0

0 v12 0

0 v13 0

0 v22 0 0 v32 0

0 v23 0 0 v33 0

0

0 v42 v43

0 0 0 0 0 v61 v62 0

0 0

0 0

0

0

0

3

7 7 7 7 7: 0 0 v44 0 0 7 7 7 0 v53 0 v55 0 5 0 0 0 0 v66 0 0

0 0

0 0

ð60Þ wji

Finally, the transversal displacements in the frame of triangle j must be expressed in terms of the displacements in the global frame by the relation wji ¼ ui  nj ;

ð61Þ

where uTi ¼ ðui ; vi ; wi Þ and nj is the normal to triangle j, Fig. 6. The transformation from the twelve transversal d.o.f. in the four subsystems to the 18 translational d.o.f., corresponding to the six nodes, in the global coordinate system is given as ws ¼ Tn d;

ð62Þ

where d T ¼ ðuT1 ; uT2 ; uT3 ; uT4 ; uT5 ; uT6 Þ and 0 2 n n n3 0 0 0 0 0 0 B 0 0 0 n3 n n1 0 0 0 B B B 0 0 0 0 0 0 n1 n n2 TTn ¼ B B0 0 0 0 0 0 0 0 0 B B @0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1

0

0

0

0 0

0 0

n1

0

0 0

n2 0

C C C C C; 0C C C 0A n3 0 0

ð63Þ where all the matrix elements are of size 3 · 1. The shell element by Guo et al. [8] is not formulated exactly like this, even though this is more exact since the rotation angles a are computed from the initially curved configuration. No more calculations are required for the bending part of their shell elements or the shell version of the element by On˜ate and Cervera [15]. To compute the shell element by Phaal and Calladine [3], the outer nodes of the curved shell must be transformed back to the plane to define C, which is replaced with CX in Eq. (46), where X is given by 1 0 1 1 1 0 0 0 0 0 0 0 0 0 B0 0 0 1 1 1 0 0 0 0 0 0C C B C B C B 0 0 0 0 0 0 1 1 1 0 0 0 T C: ð64Þ B X ¼B C B0 0 0 0 0 0 0 0 0 1 0 0C C B @0 0 0 0 0 0 0 0 0 0 1 0A 0

0

0

0 0

0

0

0 0

0

0

1

3 2

w

n3 z3

x3

6 0 6 6 6 0 Vs ¼ 6 6 0 6 6 4 v51 0

2

Fig. 5. Symmetric mesh (a) and node numbering of a non-complete patch of elements (b).

w52

0 v11 0

y3 w63

Fig. 6. For a shell the transversal d.o.f. are defined in the local coordinate systems.

2.8.2. Membrane part To obtain a complete shell element the bending part must be combined with a membrane part. Since the bending part only requires three translational d.o.f. per node, and three nodes per element, the membrane part should be equally simple. Originally, it was suggested that a CST

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

2.5

2.0

(κ x +κ y )/2

element should be used [9]. Recently it has been suggested that the six node patch can be used to compute overlapping isoparametric LST elements [5]. Elements using LST membranes are prone to membrane locking, which can be diminished by assumed strain methods with reduced numerical integration as shown in [5,18]. As the aim of this study is to investigate the sensitivity of the RF bending part to mesh distortions, only the CST element is used as the membrane part for both the RF elements and Morley’s element to isolate the effects of the mesh distortion.

1.5

1.0 Phaal and Calladine Sabourin and Brunet Oñate and Cervera Flores and Oñate

0.5

3. Results In this section, the results from classical linear benchmark examples were used to investigate the differences between the elements.

5009

0 0.5

1.0

1.5

2.0

Λ Fig. 8. The average of jx and jy as a function of the elongation K of the outer elements.

3.1. Computational cost for stiffness matrix generation

3.2. Approximation of constant curvatures The ability of the RF elements to reproduce constant curvatures is a relevant first test of their accuracy. A circular plate is bent to form a paraboloid with known constant curvature, jx = jy = 1. The curvatures are then computed over a single patch of elements, where the heights, h, of the outer elements are multiplied by the elongation factor, K, Fig. 7. The element by Phaal and Calladine reproduces the curvatures exactly, whereas the other elements diverge from the exact value as the outer elements get more distorted, Fig. 8. Even though Flores and On˜ate [5] use quadratic shape functions, constant curvatures are not

obtained since the isoparametric approximation is inaccurate when positions of the side nodes are far from the average of the vertex nodes. The formulations by On˜ate and Cervera [6] and Flores and On˜ate [5] are more sensitive to distortions than that by Sabourin and Brunet [19]. The formulation by Guo et al. [8] is identical to that by Sabourin and Brunet, so only the latter will be considered hereafter. 3.3. Square plates As a first test of the bending behaviour simply supported and clamped square plates, Fig. 9, subjected to either a distributed load q or a point load P at the centre, were investigated. Only one quarter of the square plate was analysed because of symmetry. Parameters were taken from [28], i.e. L = 1 m, E = 17.472 GPa, t = 0.1 mm, m = 0.3, q = 0.1 mN/m2 and P = 0.4 mN. For a symmetric mesh all the elements give the same results in the interior, Figs. 9–13. Differences are due to different b.c. The b.c. by Phaal and Calladine is extremely mesh dependent because of the interpretation of the symmetry b.c. For the simply supported square plate subjected to a point load, the common b.c. shows the best convergence for both meshes, whereas for

centre point

Fig. 7. A circular plate, bent to form a paraboloid with known curvature, is used to test the ability to reproduce constant curvatures for a distorted patch of elements.

symmetry

symmetry

boundary

boundary

symmetry

boundary

centre point

boundary

A comparison of the computational cost for generating the element stiffness matrices presented in the previous section was performed. The CPU time to generate 1000 random shell element in Matlab on a common workstation was measured: 1.13 s for Phaal–Calladine, 0.98 s for Guo et al. 0.96 s for Sabourin-Brunet and 0.95 s for On˜ate-Cervera. Included is the time to compute the contributions from both the bending and membrane parts. The code was not optimized, and consequently, the numbers should just be used for comparison between these elements. The generation time for the element by Phaal and Calladine is 15–20% slower than for the other elements.

symmetry

Fig. 9. Mesh of one quarter of the simply supported and clamped plates: (a) mesh A and (b) mesh B.

Phaal and Calladine b.c. (mesh A) Phaal and Calladine b.c. (mesh B) Common b.c. (mesh A) Common b.c. (mesh B) Flores and Oñate b.c. (mesh A) Flores and Oñate b.c. (mesh B) Morley (mesh A) Morley (mesh B)

1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 10

1

2

10 Number of nodes

10

Normalised displacement for centre point

1.20 1.15 1.10 1.05 1.00 0.95 0.90 0.85 0.80 10

1

2

10 Number of nodes

10

3

Normalised displacement for centre point

Fig. 11. Simply supported square plate subject to a distributed load: the vertical displacement of the centre point is normalised with 4.06 mm.

2.6

Phaal and Calladine b.c. (mesh A) Phaal and Calladine b.c. (mesh B) Common b.c. (mesh A) Common b.c. (mesh B) Morley (mesh A) Morley (mesh B)

2.4 2.2

2.0 1.8 1.6 1.4 1.2 1.0 10

Phaal and Calladine b.c. (mesh A) Phaal and Calladine b.c. (mesh B) Common b.c. (mesh A) Common b.c. (mesh B) Flores and Oñate b.c. (mesh A) Flores and Oñate b.c. (mesh B) Morley (mesh A) Morley (mesh B)

1.25

Phaal and Calladine b.c. (mesh A) Phaal and Calladine b.c. (mesh B) Common b.c. (mesh A) Common b.c. (mesh B) Morley (mesh A) Morley (mesh B)

2.2

3

Fig. 10. Simply supported square plate subject to a point load: the vertical displacement of the centre point is normalised with 11.60 mm.

1.30

Normalised displacement for centre point

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

Normalised displacement for centre point

5010

1

2

10 Number of nodes

10

Fig. 13. Clamped square plate subject to a distributed load: the vertical displacement of the centre point is normalised with 1.26 mm.

3.4. Scordelis-Lo roof The Scordelis-Lo roof [29], Fig. 15, is used to determine the elements ability to solve complex states of membrane strain. The roof is supported at each end by a rigid diaphragm (v = w = hx = 0). Because of symmetry only one quarter was modelled. A uniform vertical gravity load of q = 90 Pa is applied, L = 50 m, R = 25 m, c = 2p/9, t = 0.25 m, E = 432 MPa and m = 0. The convergence is satisfactory for all the elements, Fig. 14. No difference between the RF elements is visible if the same type of elements are used on the boundary. Since membrane effects dominate the problem, only small differences are obtained between different b.c. The results for the Flores–On˜ate and Phaal–Calladine b.c. are almost identical to those for the element by Morley for both mesh orientations, and slightly better than for the common b.c. The reason is the zero moments at the free edges for their elements, which results in a more flexible solution. The total error is therefore smaller, since this error partly cancels out the overestimated stiffness of the Kirchhoff

2.0 1.8 1.6

L y

1.4

tr me

sym

me

try

ym

s

1.2

A

1.0 1

10

2

10 Number of nodes

10

3

diaphragm

ree

f

Fig. 12. Clamped square plate subject to a point load: the vertical displacement of the centre point is normalised with 5.60 mm.

γ the uniform load case, the Flores–On˜ate b.c. produces the best results. For the clamped case the common b.c. is identical to that of Flores–On˜ate.

3

R

Fig. 14. Scordelis-Lo roof, mesh A.

Normalised displacement for point A

Normalised displacement for point A

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

1.00 0.95 0.90 0.85

Phaal and Calladine b.c. (mesh A) Phaal and Calladine b.c. (mesh B) Common b.c. (mesh A) Common b.c. (mesh B) Flores and Oñate b.c. (mesh A) Flores and Oñate b.c. (mesh B) Morley (mesh A) Morley (mesh B)

0.80 0.75 0.70 10

2

5011

1.4 1.3 1.2 1.1 1.0 0.9

Phaal and Calladine b.c. (mesh A) Phaal and Calladine b.c. (mesh B) Common b.c. (mesh A) Common b.c. (mesh B) Flores and Oñate b.c. (mesh A) Flores and Oñate b.c. (mesh B) Morley (mesh A) Morley (mesh B)

0.8 0.7 0.6 0.5

3

10 Number of nodes

10

2

3

10 Number of nodes

Fig. 15. Scordelis-Lo roof: the vertical displacement at the midside of the free edge is normalised with 0.3024 m.

Fig. 17. Pinched cylinder: the radial displacement at the point load is normalised with 18.248 lm.

assumption and the FE discretisation. Note that even faster convergence can be achieved with a higher-order membrane part, e.g. a reduced integration LST element as shown in [5,18].

was explained in Fig. 5: the stiffness of the node at the point load is overestimated for mesh A, and vice versa for mesh B. 3.6. Hemispherical shell

3.5. Pinched cylinder The pinched cylinder with a diaphragm [29], Fig. 16, is used to test inextensional bending modes and complex membrane states. The cylinder is constrained at each end by a rigid diaphragm (v = w = hx = 0). Because of symmetry only one-eighth was modelled. Opposing radial loads of P = 1.0 N are applied, L ¼ 600 m, R = 300 m, t = 3 m, E = 3 MPa and m = 0.3. Again, the only difference between the RF elements is in the b.c., Fig. 17. The common b.c. are marginally more accurate than the ones by Flores and On˜ate, and significantly more accurate than the ones by Phaal and Calladine and Morley. The b.c. by Phaal and Calladine is also significantly mesh-dependent. The explanation to the stiffer solution of mesh A, and the more flexible solution of mesh B,

symmetry L

The hemispherical shell [29], Fig. 18, tests the ability of the elements to represent inextensional modes. Moreover, it challenges the ability of the elements to handle rigidbody rotations normal to the surface of the shell. Because of symmetry only one quarter was modelled. Opposing radial point loads of P = 2 N were applied, R = 10 m, v = p/10, t = 0.04 m, E = 68.25 MPa and m = 0.3. Since the loads were applied symmetrically, changing from mesh A to mesh B is the same as swapping points A and B. The bottom circumferential edge of the hemisphere is free. The mean of the absolute values of the displacements at points A and B are shown in Fig. 19. This is the only of the four investigated benchmark examples where the structured mesh has elements of different size. Therefore, the results differ slightly between the RF elements when equal boundary elements are used. The results are normalised with the approximate value

P A

free symmetry

symmetry

R

v diaphragm

symmetry

R

P B free P Fig. 16. Pinched cylinder, mesh B.

A P

Fig. 18. Hemispherical shell, mesh B.

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

1.05

Phaal and Calladine (inside and b.c.) Phaal-Calladine (inside) and common b.c. Sabourin and Brunet (inside and b.c.) Oñate-Cervera (inside) and common b.c. Oñate-Cervera (inside) and Flores-Oñate (b.c.) Morley

1.10

1.05

Normalised vertical tip displacement

Average normalised displacement for points A and B

5012

1.00

0.95

0.90 10

2

1.00

0.95

Phaal and Calladine b.c. Common b.c. Flores and Oñate b.c. Morley

0.90

3

10 Number of nodes

0.85 10

Fig. 19. Hemispherical shell: the averages of the absolute values of the radial displacements at the loaded points are normalised with 0.0940 m [28].

0.0940 m [28], which seems to be too large. It is therefore difficult to determine which of the formulations that is more accurate, but it seems as if the b.c. by Flores–On˜ate and Phaal–Calladine converge more quickly than the over-stiff common b.c. Since only the mean value is shown, the strong mesh dependency of the b.c. by Phaal–Calladine is hidden. 3.7. Raasch’s hook The Raasch’s hook problem of a curved strip with an inplane shear load at the tip, Fig. 20, is challenging because of the inherent coupling between modes of bending, extension and twist. The left end of the hook is clamped, the other edges are free, a shear load of 8.756 N/m is applied at the tip. The dimensions and material parameters are: r = 0.3556 m, R = 1.1684 m, L = 0.508 m, c = p/6, t = 50.8 mm, E = 22.7527 MPa and m = 0.35. It is found [30] that shell elements without transverse shear flexibility converge to a value that is 5% stiffer than the solution by solid elements, whereas shell elements with transverse shear flexibility do not converge. However, if the hook is made thin, then all shell elements would converge. Here, the only visible difference between the RF elements is between different elements along the boundary, Fig. 21. In this example, the free edge b.c. has great impact

2

3

10 Number of nodes

10

4

Fig. 21. Raasch’s hook: the normalised tip displacement for grids with 35 · 6, 69 · 11, 137 · 21 and 273 · 41 nodes. The results are normalised using 0.120 m, i.e. the best available Kirchhoff solution obtained for Raasch’s hook [30].

since three of the four edges are free and few elements are used along the width. The results obtained with the b.c. by Flores–On˜ate and Phaal–Calladine resemble the results of the element by Morley. The good accuracy even for crude meshes may be due to errors of equal magnitudes, i.e. the overestimate of the stiffness due to Kirchhoff theory and the FE discretisation, and the underestimate of the stiffness from the free b.c. 3.8. Unstructured meshes All examples above have also been investigated using unstructured meshes. For the first three shell examples, grids of 21 · 21 nodes were used and for Raasch’s hook 69 · 11 nodes was taken as the norm. The results obtained with the analogous structured grids have been considered as exact. Then the interior grid points have been perturbed randomly, with values equally distributed from minus to plus the maximum perturbation, in both directions of the surface. The maximum perturbations were chosen to be 5%, 10%, 15% respectively 20% of the distance between two adjacent grid points in the structured mesh. 100 such meshes, Fig. 22, have been considered for each example,

y L

r

R x γ

Fig. 20. Raasch’s hook.

Fig. 22. Examples of grids for a coplanar square surface, randomly perturbed using maximum perturbations of (a) 10% and (b) 20% of the distance between two adjacent grid points.

0 –1 –2 –3 Phaal and Calladine Oñate and Cervera Sabourin and Brunet Morley

–4 –5 –6

0

2

4

Relative displacement error for point A (%)

2 0

–2 –4 –6 Phaal and Calladine Oñate and Cervera Sabourin and Brunet Morley

–10 –12 0

2

4

6 8 10 12 14 16 18 20 Max nodal perturbation (%)

Fig. 24. Pinched cylinder, unstructured mesh. Mean values (thick black lines) and 95% confidence intervals (thin grey lines) of the relative displacement error from 100 perturbed meshes.

and the same meshes have been used for all elements. The same elements, i.e. Sabourin and Brunet, are used along the boundary for all the RF simulations. Mean values and 95% confidence intervals of the displacement error compared to the ‘‘exact’’ value for the structured mesh are shown in Figs. 23–26. The results for the square plates showed similar behaviour and are therefore not shown. From the simulations shown here it become clear that the RF element by Phaal and Calladine is not as sensitive to mesh distortion as the other RF elements. Unexpectedly, the element by On˜ate and Cervera is less sensitive than the element by Sabourin and Brunet. A possible explanation is that taking the average of the rotation angles smooths out the effects of the mesh distortion as the geometries of the triangles in the patch have less influence.

5013

2 0 –2 –4 –6 –8

Phaal and Calladine Oñate and Cervera Sabourin and Brunet Morley

–10 –12 0

2

4

6 8 10 12 14 16 18 20 Max nodal perturbation (%)

Fig. 23. Scordelis-Lo roof, unstructured mesh. Mean values (thick black lines) and 95% confidence intervals (thin grey lines) of the relative displacement error from 100 perturbed meshes.

–8

Relative displacement error for point A (%)

1

6 8 10 12 14 16 18 20 Max nodal perturbation (%)

Fig. 25. Hemispherical shell, unstructured mesh. Mean values (thick black lines) and 95% confidence intervals (thin grey lines) of the relative displacement error from 100 perturbed meshes.

Relative tip displacement error (%)

Relative displacement error for point A (%)

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

2 0 –2 –4 –6 –8

Phaal and Calladine Oñate and Cervera Sabourin and Brunet Morley

–10 –12

0

2

4

6 8 10 12 14 16 18 20 Max nodal perturbation (%)

Fig. 26. Raasch’s hook, unstructured mesh. Mean values (thick black lines) and 95% confidence intervals (thin grey lines) of the relative displacement error from 100 perturbed meshes.

4. Discussion All the investigated RF elements use a patch of four elements with only translational d.o.f., which is enough to describe constant curvatures. However, there is a choice between two-dimensional constant curvatures over the patch of elements or superposition of one-dimensional constant curvatures for the three pairs of triangles. The first approach satisfies the Kirchhoff assumption of constant curvatures exactly, whereas the second one only approximates it. The consequence is that theoretically constant curvatures are not approximated correctly on unstructured grids when the second approach is used. Therefore, the performance of the element by Phaal and Calladine becomes significantly better for unstructured meshes. For a structured mesh the only difference between the RF elements is in the formulation of the boundary conditions. In this case, the element by Phaal and

5014

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015

Calladine is not effective since it takes more time to compute than the other elements. The element by On˜ate and Cervera is slightly faster than the others, and is preferred then. The b.c. by Flores–On˜ate gave the best results for most of the linear benchmark examples in this study. The common b.c. is too stiff for simply supported and free edges. Although the b.c. by Phaal–Calladine gave accurate results for structures with free edges, its drawbacks are: introduction of fictitious nodes and severe mesh orientation dependency for the symmetry b.c. This study suggests that the element by Phaal and Calladine is used inside an unstructured mesh, and the element by On˜ate and Cervera on the boundary with Flores–On˜ate b.c. The combination of elements has convergence properties similar to the element by Morley when the same number of nodes are used, and it is more efficient since less d.o.f. are used, but slightly more prone to severe mesh distortions. It may seem cumbersome to use two different elements, but different interpretations must be used on the boundary for all RF elements. The combination of the two elements does not imply additional problems as the same d.o.f. and approximately the same curvatures assumptions are used. It is not possible to draw any conclusions about the membrane behaviour of the RF elements from this study. To provide a comparison of the RF element in the nonlinear range, a non-linear version of the element by Phaal and Calladine must be developed. An implementation using the total or updated Lagrangian formulations may become problematic as the element does not use traditional shape functions. A more straightforward alternative is to implement the element by Phaal and Calladine in a corotational framework [31]. This requires the development of a new corotational framework for RF elements. However, the use of the corotational framework is limited to problems with large displacements and small strains. Acknowledgements

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14] [15] [16]

[17]

[18]

[19] [20]

Financial support from the Swedish Research Council and KTH is gratefully acknowledged. The authors also thank one of the reviewers for constructive and valuable comments which have helped to improve the paper.

[21]

References

[23]

[1] F.G. Flores, E. On˜ate, K. Schweizerhof, Session on advances in rotation-free shell elements, in: 5th International Conference on Computation of Shell & Spatial Structures, Salzburg, Austria, 1–4 June, 2005. [2] F. Sabourin, M. Brunet, Detailed formulation of the rotation-free triangular element ‘‘S3’’ for general purpose shell analysis, Engrg. Comput. 23 (2006) 469–502. [3] R. Phaal, C.R. Calladine, A simple class of finite elements for plate and shell problems II: An element for thin shells with only

[22]

[24]

[25] [26]

translational degrees of freedom, Int. J. Numer. Methods Engrg. 35 (1992) 979–996. H. Laurent, G. Rio, Formulation of a thin shell finite element with continuity C0 and convected frame notion, Comput. Mech. 27 (2001) 218–232. F.G. Flores, E. On˜ate, Improvements in the membrane behaviour of the three node rotation-free BST shell triangle using an assumed strain approach, Comput. Methods Appl. Mech. Engrg. 194 (2005) 907–932. E. On˜ate, M. Cervera, Derivation of thin plate bending elements with one degree of freedom per node: A simple three node triangle, Engrg. Comput. 10 (1993) 543–561. M. Brunet, F. Sabourin, A simplified triangular shell element with a necking criterion for 3-D sheet-forming analysis, J. Mater. Process. Technol. 50 (1995) 238–251. Y.Q. Guo, W. Gati, H. Naceur, J.L. Batoz, An efficient DKT rotation free shell element for springback simulation in sheet metal forming, Comput. Struct. 80 (2002) 2299–2312. M.R. Barnes, Form-finding and analysis of tension space structures by dynamic relaxation, Ph.D. thesis, Department of Civil Engineering, The City University, London, 1977. A.S.L. Chan, G.A.O. Davies, A simplified finite element model for the impact of thin shells, in: P.S. Bulson (Ed.), Structures under Shock and Impact, Proc. 1st Int. Conf., Cambridge, MA, USA, 1989, pp. 365–380. J.K. Hampshire, B.H.V. Topping, H.C. Chan, Three node triangular bending elements with one degree of freedom per node, Engrg. Comput. 9 (1992) 49–62. R. Phaal, C.R. Calladine, A simple class of finite elements for plate and shell problems I: Elements for beams and thin flat plates, Int. J. Numer. Methods Engrg. 35 (1992) 955–977. M.J. Turner, R.W. Clough, H.C. Martin, L.J. Topp, Stiffness and deflection analysis of complex structures, J. Aeronaut. Sci. 23 (9) (1956) 805–823. M. Brunet, F. Sabourin, Analysis of a rotation-free 4-node shell element, Int. J. Numer. Methods Engrg. 66 (2006) 1483–1510. E. On˜ate, F. Za´rate, Rotation-free triangular plate and shell elements, Int. J. Numer. Methods Engrg. 47 (2000) 557–603. F.G. Flores, E. On˜ate, A basic thin shell triangle with only translational DOFs for large strain plasticity, Int. J. Numer. Methods Engrg. 51 (2001) 57–83. F.G. Flores, E. On˜ate, A rotation-free shell triangle for the analysis of kinked and branching shells, Int. J. Numer. Methods Engrg. 69 (2007) 1521–1551. E. On˜ate, F.G. Flores, Advances in the formulation of the rotationfree basic shell triangle, Comput. Methods Appl. Mech. Engrg. 194 (2005) 2406–2443. F. Sabourin, M. Brunet, Analysis of plates and shells with a simplified three node triangular element, Thin-Walled Struct. 21 (1995) 209–223. J.M. Roelandt, J.L. Batoz, Shell finite element for deep drawing problems: computational aspects and results, in: IUTAM Symposium on Finite Inelastic Deformations, 1992, pp. 423–430. L.S.D. Morley, The constant-moment plate-bending element, J. Strain Anal. 6 (1971) 20–24. F. Cirak, M. Ortiz, P. Schro¨der, Subdivision surfaces: a new paradigm for thin shell finite-element analysis, Int. J. Numer. Methods Engrg. 47 (2000) 2039–2072. F. Cirak, M. Ortiz, Fully C1-conforming subdivision elements for finite deformation thin-shell analysis, Int. J. Numer. Methods Engrg. 51 (2001) 813–833. R. Hauptmann, K. Schweizerhof, A systematic development of solidshell element formulations for linear and non-linear analysis employing only displacement degrees of freedom, Int. J. Numer. Methods Engrg. 42 (1998) 49–69. P. Jetteur, Development of a new shell element for very thin shell structures, Tech. Rep., Samtech, Liege, Belgium, Report 236, 2002. M. Ga¨rdsback, G. Tibert, Evaluation of triangular shell elements for thin membrane structures, in: Proc. 5th International Conference on

M. Ga¨rdsback, G. Tibert / Comput. Methods Appl. Mech. Engrg. 196 (2007) 5001–5015 Computation of Shell & Spatial Structures, Salzburg, Austria, 1–4 June, 2005. [27] F.G. Flores, E. On˜ate, Rotation-free element for the non-linear analysis of beams and axisymmetric shells, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5297–5315. [28] R.H. MacNeal, R.L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Elem. Anal. Des. 1 (1985) 3–22.

5015

[29] T. Belytschko, H. Stolarski, W.K. Liu, N. Carpenter, J.S.J. Ong, Stress projection for membrane and shear locking in shell finite elements, Comput. Methods Appl. Mech. Engrg. 51 (1985) 221–258. [30] N.F. Knight Jr., Raasch challenge for shell elements, AIAA J. 2 (1997) 375–381. [31] J.-M. Battini, A modified corotational framework for triangular shell elements, Comput. Methods Appl. Mech. Engrg. 196 (2007) 1905–1914.