Computers and Structures 75 (2000) 321334

www.elsevier.com/locate/compstruc

A unied approach for shear-locking-free triangular and

rectangular shell nite elements
Kai-Uwe Bletzinger a,*, Manfred Bischo b, Ekkehard Ramm b
a
Institute of Structural Analysis, University of Karlsruhe Kaiserstr 12, 76131, Karlsruhe, Germany
b
Institute of Structural Mechanics, University of Stuttgart, 70550 Stuttgart, Germany

Abstract

A new concept for the construction of locking-free nite elements for bending of shear deformable plates and
shells, called DSG (Discrete Shear Gap) method, is presented. The method is based on a pure displacement
formulation and utilizes only the usual displacement and rotational degrees of freedom (dof) at the nodes, without
additional internal parameters, bubble modes, edge rotations or whatever. One unique rule is derived which can be
applied to both triangular and rectangular elements of arbitrary polynomial order. Due to the nature of the method,
the order of numerical integration can be reduced, thus the elements are actually cheaper than displacement
elements with respect to computation time. The resulting triangular elements prove to perform particularly well in
comparison with existing elements. The rectangular elements have a certain relation to the Assumed Natural Strain
(ANS) or MITC-elements, in the case of a bilinear interpolation, they are even identical. # 2000 Elsevier Science
Ltd. All rights reserved.

Keywords: Finite element; Reissner/Mindlin shells and plates; Triangular and rectangular; Shear-locking; ANS; Linear elastic

1. Introduction are either the most ecient in the experience of the

authors, or have some similarity to the method itself.
The development of ecient shear deformable plate Plate and shell formulations of the KirchhoLove
and shell elements of the ReissnerMindlin type has a type, without consideration of transverse shear eects,
more than 30 year old history and it is completely im- are not addressed in the present study.
possible to give an overview over the innumerable con- Regarding the multitude of dierent concepts, it can
cepts, invented by both mathematicians and engineers be recognized that most element developers concen-
in the past. Therefore, the method described in the pre- trate their eorts either on triangular elements or on
sent paper is solely compared to those concepts that quadrilaterals. There are only few papers, where a suc-
cessful concept for quadrilaterals is transferred to tri-
angles or vice versa.
* Corresponding author. Tel.: +49-721-608-2283; fax: +49-
It is also apparent, that the problem of transverse
721-608-6015. shear locking is practically solved for structured
E-mail addresses: kub@bau-verm.uni-karlsruhe.de (K.-U. meshes with regular element shapes. Here, most of the
Bletzinger), bischo@statik.uni-stuttgart.de (M. Bischo), elements described in the past 25 years perform well.
eramm@statik.uni-stuttgart.de (E. Ramm). However, in the case of distorted meshes (e.g. for a

0045-7949/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved.
PII: S 0 0 4 5 - 7 9 4 9 ( 9 9 ) 0 0 1 4 0 - 6
322 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

complex geometry or when automatic meshers within The method is based on interpolation of the shear
adaptivity are used) there are still some problems. strains from particularly chosen sampling points and
The question, whether triangles or quadrilaterals are successfully eliminates their parasitic part. Until today,
the `better' choice, still seems to be not yet decided. the MITC4 element is probably the most ecient
While most of the quadrilaterals exhibit better per- bilinear element for the analysis of both thick and thin
formance concerning convergence rates, triangles are plates and shells. Recent developments concern the re-
denitely easier to apply for free-meshing algorithms, duction of distortion sensitivity of the MITC4 element
and therefore, preferred in adaptivity. through stabilization methods [17].
For triangular elements it is remarkable that many It seems, however, that the transfer of both the KM
formulations contain awkward procedures while deriv- and the ANS concept to triangles is not trivial. Es-
ing the element stiness matrix. Especially, when ad- pecially a proper choice of feasible sampling points in
ditional dof are introduced (e.g. rotations at the mid- the ANS formulation proves to be more problematic
points of the edges, bubble modes, etc.) and condensed than for rectangles. One of the rst successful linear
out later on to preserve the global number of dof, the triangular elements has been developed by Xu [27] on
question arises, if a similar result could be obtained by the basis of the KM concept, introducing additional
simpler procedures. dof in the element center (`bubble modes').
In the present paper, a methodology is described Although the method presented in this study could
which allows the formulation of ecient nite elements be classied as an ANS method, from the point of
of arbitrary polynomial order, either triangular or rec- view of the authors, there are some advantages. Due
tangular, with one unique, simple rule. The method is to the fact that the element formulation evolves in a
based on the explicit satisfaction of the kinematic natural way for any kind of element, regardless of
equation for the shear strains at discrete points and shape and polynomial order, there is no need to
eectively eliminates the parasitic shear strains. The choose an interpolation for the shear strains or to
essential step is the calculation of discrete shear gaps specify any sampling points. In the case of rectangles,
(DSG) at the nodes and their interpolation across the the present concept leads to the same stiness matrices
element domain, thus obtaining a shear strain distri- as the KM or ANS elements, respectively. For tri-
bution which is free of parasitic parts angles, no equivalent could be found in the literature.
The concept could be regarded as a B-bar method, This leads us to the opinion that the present DSG-el-
because only the dierential operator for the strain-dis- ements might be the missing link between rectangular
placement relation is aected. The formulation uses and triangular KM or ANS elements.
the standard dof of pure displacement elements and Due to the multitude of dierent concepts to elimin-
does not introduce extra nodes or internal parameters. ate shear locking in beams, plates and shells this short
The only modication with respect to displacement el- review is necessarily incomplete. To sum up, one can
ements is the dierent calculation of transverse shear say, that the basic idea of the DSG concept appears in
strains. This, in turn, makes it most easy to implement the literature in numerous dierent shapes. The main
the element into an existing code. contribution of the present paper is the systematic
The resulting elements are free of locking, pass the development of a class of ecient elements by straight-
patch-test, and show reduced sensitivity to mesh dis- forward realization of the basic concept, described in
tortions. The computation time for the construction of Section 2.
the element stiness matrix is less than for pure displa-
cement elements, which makes the method extremely
ecient.
2. The basic idea
The most apparent similarities to the present method
can be observed in the so-called Kirchho mode (KM)
The basic idea is most simply explained with the
concept, originally proposed by Hughes and Tezduyar
example of the Timoshenko beam theory. The defor-
[16] (see Hughes and Taylor [15] for a corresponding
mation of the beam continuum is described by the dis-
linear triangle). Here, conditions for the shear strains
placement vx of the beam center line and rotation
are formulated along the edges of the element and the
bx of the cross section, Fig. 1. The dierence of ro-
resulting discrete shear strains at the nodes are interp-
tation bx and the gradient of the displacement v 0 x
olated over the element domain with the standard
denes the shear deformation gx at any point x along
shape functions.
the beam:
Although initially not realized, the so-called ANS
(Assumed Natural Strain) or MITC (Mixed Interp- gx v 0 x bx 1
olation of Tensorial Components) approach of Bathe
and co-workers [4,8] leads to identical elements as the The special case of pure bending is reected by the so-
KM concept in the case of a bilinear interpolation. called Bernoulli condition (Kirchho for plates), i.e.
 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334 323

Fig. 1. Timoshenko beam.

that the shear deformation has to vanish: The discrete shear gap is dened at node i by inte-
gration of the discretized shear strains:
0  v 0 x bx 2
x i x i
Dvig x 1 vh jxx i1 bh dx gh dx 4
The total displacement of the beam is due to defor- x1 x1
mation with respect to bending and shear. The shear
related part is determined by integration of (1): where x i is the coordinate of node i. For the case of
pure bending, the discretized Bernoulli (or Kirchho)
x^ x^
condition means zero discrete shear gaps. This con-
Dvg x^ g dx vjxx^ 0 b dx dition can be fullled, leading to the correct zero shear
x0 x0
deformation. Although formulated for shear deform-
^
vx-vx 0 able beams the concept meets the discrete Kirchho
|{z} Dvb 3 idea for the case of pure bending.
Dv After discretization of displacement and rotation
eld, in general, (1) does not apply anymore to deter-
which describes the increase of displacement due to mine the discretized shear deformation gh :
shear between the positions x 0 and x: ^ Dvg can be
identied as the `shear gap', the dierence between gh 6 vh0 bh 5
the increase of the actual displacement Dv and the dis-
placement Dvb which corresponds to pure bending, This is obvious in the case of pure bending. In particu-
Fig. 2. lar, if displacement vh and rotation bh are interpolated

Fig. 2. Shear gap.

324 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

by the same functions as it is usually the case

the dierence of displacement gradient and rotation
cannot anymore vanish identically, i.e. the Bernoulli
(Kirchho) condition cannot be satised for pure
bending:

0 6 vh0 bh 6

If, however, (5) is used although it is not valid but

is usually done with standard displacement elements
any deformation, even pure bending, exhibits some
Fig. 3. Linear Timoshenko beam.
parasitic shear. The element behaves too sti; the eect
is known as `shear locking'. l
As a consequence, the discretized shear deformation Dv1g 0; Dv2g v2 v1 bh dx 10
has to be formulated alternatively, e.g. in an integral 0
sense by discretization of the shear gap (3). The idea is
The discretized rotations are given by
to represent the shear deformation by their equivalent
part of the nodal displacements, the discrete shear X2
1 1
gaps Dvig at the nite element nodes. bh x N i
bi l x b1 xb2 11
l l
The discretized shear gap eld is determined by in- i1
terpolation of the nodal shear gaps:
and with (10) we obtain
X
N
l l
Dvg N i Dvig 7
i1
Dv2g v2 v1 N 1 dx
b1 N 2 dx
b2
0 0

N is the number of element nodes and N i are the l 1

shape functions. The shear deformation is evaluated v2 v1 b b2 12
2
straightforward by dierentiation:
The distribution of the shear gap across the element is
X
N
dN i now calculated by interpolation from its nodal value
gh Dvig 8 with the standard shape functions N i ,
i1
dx

By that means the shear deformation is consistently X

2
Dvg x N i
Dvig
separated from the bending deformation and, there- i1
fore, the procedure could also be understood as a de-  
composition of shear and bending modes. However, 1 l

x v2 v1 b1 b2 13
this is the result of the formulation and was not a pre- l 2
sumption as in other comparable approaches. In prin-
ciple, the idea has been brought up in earlier Finally, the `correct shear' is obtained by dierentiat-
contributions for beams as well as it is part of special ing Dvg with respect to x.
formulations for plate and shell elements, as indicated
dDvg x 1 1

in the introduction. The dierence to the present con- gx v2 v1 b1 b2 14
dx l 2
cept is that now the procedure can be unied and
equivalently applied to plates and shells. This will be Modication of the dierential operator for the
demonstrated in the following sections. strains according to (14) leads to a shear-locking-free
The idea is put in concrete form by the example of a beam element. In fact, the resulting stiness matrix is
simple linear Timoshenko beam element, Fig. 3. identical to that of a reduced integrated beam element.
According to Eq. (4), the shear gap is evaluated at For comparison, the shear strains in the pure displace-
node i with coordinate x i x 1 0; x 2 l : ment element are
x i x i
1 2 1 x
gh x dx vh jx0 i bh x dx Dvig 9 gd x v v1 l x b1
b2 15
0 0 l l l
Since the shear strains alone are regarded, only the Note, that only the part containing the rotations is
shear gap dierence is of interest, the gap at node 1 aected by the procedure and that the undesired linear
which can be identied as the integration constant is components in Eq. (15) do not show up in Eq. (14). In
set to zero. Hence the case of curved beams (or shells, see Section 3) an
 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334 325

additional term evolves, that also aects the part con- use of Einstein's summation convention. A bar denotes
taining the displacements vi : variables in the deformed conguration.
The geometry of the shell in the undeformed and the
deformed state is represented by
3. DSG-elements in curvilinear coordinates h a1  a2
x r y3
a3 ; a3
; aa r,a 18
2 ja1  a2 j
3.1. Shell formulation including thickness stretch

The basic idea described above is now used to derive a a r ,a aa v,a 19

a class of shear deformable plate and shell elements,
either triangular or rectangular, of arbitrary poly- a 3 a3 w 20
nomial order. The underlying shell formulation is a
seven-parameter model, including a thickness stretch of i.e. the position vector x to any arbitrary point of the
the shell. This model is feasible for implementation of shell can be expressed by the position vector r to a
arbitrary three-dimensional constitutive laws, and has reference point on the midsurface of the shell and the
been described extensively in Refs. [7,9]. The procedure so-called `director' a3 , Fig. 4. The position vector of a
to obtain the modied dierential operator is exactly point in the deformed conguration is given by
the same as in Section 2.
As the shell formulation at hand covers the fully x r y3
a 3 21
three-dimensional stress and strain state, the strains
can be expressed by the three-dimensional Green and thus, the deformation can be described as
Lagrange strain tensor E
u v y3
w 4 x x u; r r v;
i j 1
22
E Eij
g
g ; Eij g i g j gi gj 16 a 3 a3 w
2
Here, gi and gi are the co- and contravariant base vec- In contrast to a classical ReissnerMindlin type ve-
tors of the shell body, respectively. The covariant base parameter plate or shell formulation, the update of the
vectors in the undeformed and deformed conguration, director is formulated through the dierence vector w
respectively, are given by instead of a rotation tensor. The three independent
@x components of the dierence vector allow for a change
ga x,a ; g a x ,a 17 of direction of the director a3 as well as a change of its
@ ya
length.
Note that Greek indices run from 1 to 2, Latin indices To avoid a certain `thickness locking' phenomenon,
run from 1 to 3. Unless otherwise stated, we will make the formulation with six parameters has to be further

Fig. 4. Geometry and kinematics of the shell.

326 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

enhanced by a seventh parameter, introducing a linear shell surface is calculated:

varying normal strain in thickness direction. This can
be achieved either indirectly, by rising the polynomial h an1  an2
a n3 32
order of the displacements [18,22] or directly via a 2 jan1  an2 j
hybrid mixed method. A possibility to introduce this
Then, the dierent directors a n3 of common nodes of
parameter with the help of the Enhanced Assumed
adjacent elements are averaged to yield one director
Strain (EAS) method has been described by Buchter
per node:
and Ramm [9]. For details see also Buchter et al. [10]
and Bischo and Ramm [7].
1 X nel
Another thickness locking phenomenon, termed an3
a n ; nel no: of adjacent elements 33
`curvature thickness locking', can be overcome with nel n1 3
the help of an ANS approach [6,7,21]. It should be
noted, that both thickness locking eects are conse- The director eld within one element can now be
quences of the special shell formulation and do not expressed as:
appear in a `classical' ve-parameter shell formulation. X
N X
N
For the calculation of the strain, the covariant base a3 N n
an3 ; a 3 a3 N n
wn 34
vectors of the shell body are needed in terms of the n1 n1
mid-surface.
The constant part of the transverse shear strains is due
g i x ,i x,i u,i gi u,i 23 to (26)

1
aa3 aa a 3 aa a3
g a ga v,a y3 w,a ; g 3 g3 w 24 2

With 1 
aa v,a a3 w aa a3 35
2
3
ga x,a aa y a3, a ; g3 a3 25
For geometrically linear problems the terms, which are
we nally obtain expressions for the strain tensor com- quadratic in the displacements, are neglected, thus
ponents as Eij 1aij y3 bij , with 1
aa3 aa w a3 va 36
1
2
aij a i a j ai aj 26
2 The linear part of the transverse shear strain ba3 does
not contribute to the shear locking phenomenon and
1
remains, therefore, unchanged.
bab a a a 3, b a b a 3, a aa a3, b ab a3, a 27
2

1

ba3 a 3, a a 3 a3, a a3 ; b33 0 28 3.2. Modication of shear strains for DSG elements
2
Discretization of geometry and displacements Now, the same steps as for the derivation of the
beam element in Section 2 are performed. First, the
X
N X
N X
N
r N n xn ; v N n vn ; w N n wn 29 discrete shear gaps are evaluated by integrating the
n1 n1 n1 transverse shear strains (36) over the element domain.
xi Zi
leads to the covariant base vectors of the midsurface in
Dvig1 a13 dx; Dvig2 a23 dZ 37
discrete form x1 Z1

X
N
@Nn X
N
@Nn Here, xi ; Zi are the natural coordinates of node i, Dviga
a1 xn ; a2 xn 30
@x @Z is the discrete shear gap associated with the a-direc-
n1 n1
tion. Introducing (36) into (37) yields
xi
X
N
@Nn X
N
@Nn 1
a 1 a1 n
v ; a 2 a2 v n
31 Dvig1 a3 v,1 a1 w dx
@x @Z 2 x1
n1 n1
xi
The shell `director' is calculated as follows. First, at 1
a3 v,1 r,1 w dx Dviv1 Dviw1 38
each element node the normal vector to the discretized 2 x1
 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334 327
Zi
1
Dvig2 a3 v,2 a2 w dZ triangles or rectangles of arbitrary polynomial order.
2 Z1 The only modication that has to be carried out to
Zi obtain a DSG-element in an existing code for pure dis-
1
a3 v,2 r,2 w dZ Dviv2 Dviw2 39 placement elements, is to replace the corresponding
2 Z1
part in the dierential operator matrix B by Eqs. (47)
The decisive advantage of the procedure is that the and (48).
integrals in Eqs. (38) and (39) can be determined ana- It should be mentioned that the formulae have been
lytically a priori, e.g. by the proper use of computer derived without considering conforming displacements
algebra packages (see Section 4.1 for a three-node el- at the element interfaces, thus violating the principle of
ement). This leads to a very compact and ecient pro- virtual work. However, in the case of quadrilateral el-
gram code. ements the formulation leads to well known and
Discretization of the displacement vectors v and w accepted results, as there is a variational justication
for ANS elements (see Ref. [26]). For triangular el-
X
N X
N ements the formulation still misses a rigid mathemat-
v N n x, Zvn ; w N n x, Zwn 40 ical justication, however, the elements pass the patch
n1 n1 test and are free of spurious kinematic modes.
Further, the idea of discrete shear gaps introduces
leads to the following expressions for the discrete shear
nodal indicators of shear deformation or, for thin el-
gap at node i (coordinates xi ; Zi ).
ements, discrete Kirchho constraints at the element
" N #
nodes. This means that for any specic element the
1 xi X @Nn
XN


Dviv1 x, Zi vn
N n x, Zi an3 dx 41 constraint count [14] has always the ideal number.
2 x1 n1 @ x n1

Zi " X
N n X
N
#
1
Dviv2 n
xi , Zv
n n
N xi , Za3 dZ 42
2 Z1 n1
@Z n1
4. Element matrices

xi " X #
1 N
@Nn
n X N

n 4.1. Three-node triangular element
Dviw1 x, Zi r
n
N x, Zi w dx 43
2 x1 n1
@x n1
4.1.1. Curvilinear coordinates
" N # As an example, a three-node DSG-element is de-
1 Zi X @Nn XN
rived.
Dviw2 n
x , Zr
n n
N xi , Zw dZ 44
2 Z1 n1 @ Z i n1
Geometry and shape functions, as well as their de-
rivatives are given for a three-node element in Fig. 5.
According to Eq. (13) the next step is the interpolation The discrete shear gaps at the nodes are calculated
of the discrete shear gaps across the element. according to Eqs. (41)(44)

X
N 1 2


Dv1v1 0; Dv2v1 v v1
a13 a23 ; Dv3v1 0
Dvg1 x, Z N n x, Z
Dvng1 45 4
n1

1 3

X
N Dv1v2 0; Dv3v2 v v1
a13 a33 ; Dv2v2 0
n
4
Dvg2 x, Z N x, Z
Dvng2 46
n1
1 2

The modied shear strains are nally obtained via par- Dv1w1 0; Dv2w1 r r1
w1 w2 ; Dv3w1 0
4
tial dierentiation of the shear gap distribution.
N   1 3

@ Dvg1 x, Z X @Nn Dv1w2 0; Dv3w2 r r1
w1 w3 ; Dv2w2 0
gx
Dvng1 47 4
@x n1
@x
The modied shear strains for the three node elements
are consequently
N 
X 
@ Dvg2 x, Z @Nn h
gZ
Dvng2 48 1

i
@Z n1
@Z gx v2 v1
a13 a23 r2 r1
w1 w2
4
These formulae apply to any kind of element, either 49
328 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

Fig. 5. Three-node element.

h

i
1 3
gZ v v1
a13 a33 r3 r1
w1 w3 bilinear, reduced integrated quadrilateral with hourglas
4 stabilization. The shear strain is derived from a modi-
50 cation of the BatheDvorkin ANS-scheme with respect
to the hourglass modes of shear deformation.
In contrast to the pure displacement element, where
the shear strains are linearly varying across the el- 4.1.2. Plate element in Cartesian coordinates
ement, they are constant in this formulation. The spur- To demonstrate the simplicity and eciency of the
ious linear components which are responsible for shear presented method, a closed form B-operator matrix for
locking are eectively eliminated, independent of the the DSG3 plate element is given in this section. It
element shape. relies on a classical ve-parameter formulation with ro-
As for the linear beam element and the bi-linear tations bx and by instead of the dierence vector used
plate element, the procedure could be interpreted as an in Section 3.
averaging of the transverse shear strains along the The element is dened as shown in Fig. 6. Geome-
edges [16], this is obviously not true in the case of the try, deections and rotations are interpolated by the
linear triangle, where, only two edges are taken into same shape functions:
account. This is the decisive dierence to the KM-el- 8 9 8 9
ement, introduced by Hughes and Taylor [15], where >
> x >> >
> xi >
>
>
>y > >
> >
> yi >>
explicit satisfaction of the Kirchho condition is intro- >
< = X 3 >
< > = N1 1 x Z
duced along all three edges leading to an articial con- v Ni vi ; N2 x 51
straint. The resulting element is consequently not free >
> b >
> >
> i >
>
>
> >
x> i1 >
> bx>> N 3 Z
of shear locking. >
: by >; >
: bi >;
y
As a result, the element stiness matrices of triangu-
lar DSG-elements depend on the sequence of node
numbers (see Eq. (57)). The inuence on the solution,
however, diminishes with mesh renement as satisfac-
tion of the patch test is ensured anyhow. It should be
noted that this unusual loss of objectivity of the el-
ement formulation is in no way a result of any specic
construction for three node elements, but comes along
quite naturally with a consequent application of the
DSG-concept. In fact, it can be observed that the
attempt to treat all three edges of a triangle equally
turns out to be a major obstacle in the construction of
locking-free elements. Therefore, the DSG-concept
merely takes care of coordinate directions instead of el-
ement edges.
Additionally, an alternative formulation of a three
node element insensitive to node numbering should be
mentioned, which is given in the ABAQUS theory
manual [13] and was presented by Fox and Nagtegaal
[11]. The element is dened as the result of a collapsed Fig. 6. Three-node element.
 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334 329

The Jacobian matrix and its inverse are determined to as a piecewise plane `facet element' also for the analy-
be: sis of shells, provided no average director (see Section
3) is used.
   
x ,x y,x a b At the nodes the discrete shear gaps are evaluated to
J be
x ,Z y,Z d c
   ;
x Z,x 1 c b 52 Dv1g1 Dv3g1 Dv1g2 Dv2g2 0
J1 ,x
x,y Z,y det J d a
det J ac bd 2A 1   1  
Dv2g1 v2 v1 a b1x b2x b b1y b2y
2 2

Next, the `shear gap' distributions are calculated: 1   1  

Dv3g2 v3 v1 d b1x b3x c b1y b3y 55
2 2
xi


Dvg1 vx, Zjxxi1 bx a by b dx from which the modied interpolation scheme for the
x1
shear strain results:

@ N2 @ x 2 @ N3 @ Z 3
   xi gx Dv Dv
1 1 @ x @ x g1 @ Z @ x g1
vx, Z a x x2 xZ b1x x2 b2x xZb3x
2 2 x1
@ N2 @ x 2 @ N3 @ Z 3
gy Dv Dv 56
@ x @ y g2 @ Z @ y g2
   xi
1 1 Finally, the dierential operator B can be derived to
b x x2 xZ b1y x2 b2y xZb3y 53
2 2 x1 determine curvature and shear deformations from
nodal displacements and rotations:
and T
kx , ky , kxy , gx , gy Bu
Zi

 T
Dvg2 vx, ZjZZi1 bx d by c dx B w1 , b1x , b1y , w2 , b2x , b2y , w3 ,b3x , b3y
Z1

2   3
0 bc 0  0 c 0  0 b 0
 
6 0 0 da  0 0 d  0 0 a 7
6   7
6 0 da bc  0 d c  0 a b 7
1 6   7
6   7
B 6 1  ac bc  bd bc 7 57
det J 6 b c det J 0  c  b 7
6 2  2 2  2 2 7
6   7
4 1  ad bd  ad ac 5
da 0 det J  d  a

2 2 2 2 2

   Zi The stiness matrix K and the stresses s are deter-

1 1 mined by the standard nite element relations:
vx,Z d Z xZ Z2 b1x xZb2x Z2 b3x
2 2 Z1
K BT DB dV; s DBu 58
   Zi V
1 2 1 2 1 2 3
b Z xZ Z by xZby Z by 54
2 2 Z1 where D is the constitutive matrix of plate bending and
u is the vector of generalized nodal displacements. Note,
It has been assumed that the axes of rotation for bx that the proposed formulation results in a simple modi-
and by do not change within the element. Thus, the cation of the B-operator compared to the original dis-
formulation is, strictly speaking, not valid for shells, placement element. In the special case of the three node
although the described element accounts for membrane element, the B-operator is even a constant matrix, i.e.
action. However, the present element could be applied the volume integration in (58) can be performed analyti-
330 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

cally or numerically by a one point Gauss quadrature. several reasons speaking for the use of higher order el-
Again, compared to the original displacement formu- ements. Given a certain number of dof, the accuracy
lation, the order of integration can be reduced without of the results is usually signicantly better when using
activating zero energy modes. This result can be general- higher order elements. In addition, higher order el-
ized, i.e. the order of integration can be reduced by one ements are not as sensitive with respect to linear mesh
for any triangle or selectively for the integration of the distortions. Quadratically distorted meshes can be
shear deformation parts of quadrilateral elements which avoided by a sub-parametric interpolation of the geo-
leads to considerable enhancements of eciency. metry in the geometrical linear case.
It can be seen from the numerical investigations in
4.2. Four-node quadrilateral Section 5.1, that the DSG6-element exhibits a tremen-
dous rate of convergence, although the elements are
Analogous to the development of a three-node el- quadratically distorted due to the circular shape of the
ement, a bilinear four-node element can be derived, structure.
Fig. 7. Rectangular elements are usually superior to
triangles with respect to the rate of convergence. 4.4. Nine-node quadrilateral
Application of Eqs. (47) and (48) leads to the follow-
ing shear strain distributions for the DSG4 element. The nine node DSG9-element has again some simi-
h i larities to the corresponding ANS-element (e.g. ref
1


gx 1 Z v2 a23 v1 a13 w1 w2 r2 r1 [20]). For rectangular and linearly distorted elements,
4 the stiness matrices are again identical (cf. Section
h

i 4.3), for quadratic distortions the stiness matrices are
1
1 Z v3 a33 v4 a43 w3 w4 r3 r4 slightly dierent. However, these dierences are practi-
4
(59) cally not signicant.
The main merit of the nine-node element is robust-
h

i ness, rather than eciency. Even in the case of quadra-
1
gZ 1 x v4 a43 v1 a13 w1 w4 r4 r1 tically distorted meshes, performance is excellent. For
4
h i the application to shells, a modication of the mem-
1


1 x v3 a33 v2 a23 w3 w2 r3 r2 brane part is recommended to avoid in-plane shear
4 locking and membrane locking. To achieve this, one
(60) possibility is the use of the Enhanced Assumed Strain
Note, that these linear-constant distributions of shear (EAS) method, introduced by Simo and Rifai [23,24],
strains are exactly the ones obtained by application of and rst applied to four-node, linear shell elements by
the classical ANS method [4,16]. In fact, the resulting Andelnger and Ramm [1]. An application of the
stiness matrix is identical to the one of the `Bathe method, leading to extremely robust and reliable nine-
Dvorkin' element, even for arbitrarily distorted el- node shell elements for geometrical non-linearity has
ements. been described by Bischo and Ramm [7].

4.3. Six-node triangle

5. Numerical investigations
In practical engineering applications low-order el-
ements are usually preferred. Nevertheless, there are In order to examine the presented elements with

Fig. 7. Four-node element.

 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334 331

respect to eciency, reference is made to some of the 5.1. Circular plate

most successful triangular and rectangular elements
known to the authors. This simple example has two advantages. Firstly, an
The triangular element of Xu [27] was the rst lock- analytical solution for the Kirchho theory is avail-
ing-free triangular element with nine dof. The ad- able, and second, although the geometry is regular, the
ditional dof of the nodes at the midpoints and the two mesh is `automatically' distorted, and there is no need
bubble modes for the rotations can be condensed on for arbitrary, articial mesh distortions to test the sen-
the element level. Thus, the total number of 9 dof per sitivity of the elements (Fig. 8). The analytical solution
element is preserved. is [25]:
The DST element (Discrete Shear Triangle) of Batoz
and Katili [5], is an extension of the successful DKT pr4
element (`Discrete Kirchho Triangle'), introducing the w ;
64D
eect of transverse shear deformations. The derivation 61
of the element matrices is very elaborate; the basic fea- Et3 1
D ) w10:010731 3 10:731
ture is the calculation of the shear strains from the 121 n2 t
bending moments. Thus, shear locking eects are com-
pletely avoided, but convergence of the shear strains is The elements are tested for an extremely thin plate
poor. (slenderness 1:500) and compared to some of the
As a member of the ANS- or MITC-family, the six- most popular and ecient elements known to the
node MITC7 element of Bathe et al. [3] is also exam- authors. Certainly, the given slenderness is outside
ined. The basic idea for the element formulation fol- the range of practical signicance and applicability
lows the concept of the four-node element of Bathe of a linear plate theory. This geometry is merely
and Dvorkin [4]. In addition to the mixed interpolation chosen to have a strong tendency to shear locking,
of shear strains, a bubble mode is introduced. It thus demonstrating the benet of the proposed
should be mentioned that the computational cost for method. Using symmetry, only one quarter of the
the calculation of the element stiness matrix is very system has been analyzed. The results are compiled
high, furthermore, a six-point quadrature is necessary in Table 1.
for its numerical integration (eigenvalue analyses, how- It can be seen that the present elements are among
ever, show that even a four point quadrature does not the best available triangular plate bending elements
produce any kinematic modes). In the present study, a available to date. Especially the six-node element
six-point integration rule is applied. shows a rapid convergence to the nal solution. Note
The only rectangular element we refer to is the well that the Kirchho solution neglects the inuence of
known MITC4 element of Bathe and Dvorkin [4] transverse shear deformations, thus the analytical sol-
which is based on the ANS method (see also Ref. [16]), ution for a Reissner/Mindlin kinematics leads to
already discussed in Section 1. slightly larger deformations.

Fig. 8. Circular plate under uniform load-system and discretization.

332 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

Table 1
Circular plate results

Number of nodes DSG-3 (3-node) Xu (3-node) DST (3-node) DSG-4 (4-node) DSG-6 (6-node) MITC7 (6-node)

19 10.324 9.844 10.949 10.653 10.828 8.479

61 10.632 10.514 10.793 10.715 10.739 10.162
127 10.687 10.636 10.759 10.724 10.733 10.481
217 10.707 10.678 10.747 10.728 10.732 10.591
331 10.716 10.698 10.742 10.729 10.732 10.642
469 10.721 10.708 10.739 10.730 10.732 10.670
631 10.724 10.714 10.737 10.730 10.732 10.686

In Fig. 9 the center displacement of the plate is angular DST and Xu elements are added. Due to the
plotted versus the number of nodes. On the left dia- fact that these are based on a classical ve-parameter
gram it can be seen that all curves approach the formulation, there are slight deviations in the nal
Kirchho solution with increasing number of nodes. results for the rather thick shell as can be seen on the
The dierent scale on the right diagram makes the right-hand side in Fig. 10. Nevertheless, it can be con-
excellent performance of the DSG-elements even more cluded that also in this example the DSG-elements are
obvious. Here, also the results of the DSG4-element competitive with well established formulations.
(identical to those of the ANS element of Bathe and While comparing the results, one should be aware
Dvorkin [4]) are added, demonstrating the superiority that for a xed number of nodes, the Xu and DST el-
of rectangular to triangular elements in the rate of con- ements are signicantly more expensive with respect to
vergence. computation time. The reason for that is the fact that
more quadrature points are needed for proper inte-
5.2. Cylindrical shell (`Scordelis-Lo Roof') gration of the element matrices and additional dof are
involved, which have to be condensed out on the el-
This cylindrical shell under dead load is an often ement level.
used benchmark problem for linear and non-linear The reader might recognize that the overall results
shell analysis. One advantage is, that in contrast to of all elements tested in this example are relatively
many other benchmarks there are no singularities poor. This, however, is due to the fact that for the
involved. sake of comparability no additional eorts have
The shell is analyzed using dierent discretizations been undertaken to improve the membrane behavior
with DSG3- and DSG4-elements in the setting of the of the elements and to remove membrane locking. For
three-dimensional shell formulation described in Sec- quadrilateral elements this could be achieved by appli-
tion 3. For comparison, results obtained by the tri- cation of the EAS method. For triangles, a combi-

Fig. 9. Circular plate results.

 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334 333

Fig. 10. Circular plate results.

nation of EAS and the so-called `free formulation' is 6. the order of numerical stiness integration can be
successful (see Ref. [12] for an overview). reduced without activating zero energy modes; these
together allow to generate an ecient element code
which can be derived from existing code for displa-
6. Conclusions cement elements by some simple modications of
the B-matrices.
A unied formulation of shear-locking-free, Reiss- So far the presented formulation is restricted to linear
ner/Mindlin plate and shell, rectangular and triangular problems. Further developments towards fully geome-
nite elements has been presented. The method is trically nonlinear behavior are in progress.
based on the decomposition of bending and shear de-
formation. The presence of shear is identied by a
References
nodal criterion. It compares the actual nodal displace-
ments and displacements which are related to a pure
[1] Andelnger U, Ramm E. EAS-elements for two-dimen-
bending mode. The dierence in displacements are the
sional, three-dimensional, plate and shell structures and
so called `shear gaps'. In turn, the shear strains are
their equivalence to HR-elements. International Journal
determined from the interpolated shear gaps using the for Numerical Methods in Engineering 1993;36:131137.
standard element shape functions. Compared to stan- [3] Bathe KJ, Brezzi F, Cho SW. The MITC7 and MITC9
dard displacement formulations the method results in plate bending elements. Computers and Structures
a modication of the shear related part of the B-oper- 1989;32:797814.
ator matrix. DSG-elements can be classied as ANS- [4] Bathe KJ, Dvorkin EN. A four-node plate bending el-
type or B-bar-elements. ement based on Mindlin/Reissner theory and a mixed in-
DSG-elements have several superior properties: terpolation. International Journal for Numerical
Methods in Engineering 1985;21:36783.
1. the formulation is the same for any triangular or [5] Batoz JL, Katili I. On a simple triangular Reissner/
rectangular element, Mindlin plate element based on incompatible modes and
2. the elements pass the patch test, discrete constraints. International Journal for Numerical
3. applications of the method to three-, four-, six-, and Methods in Engineering 1992;35:160332.
nine-node elements show that the behavior is either [6] Betsch P, Gruttmann F, Stein E. A 4-node nite shell el-
identical or even better than the best existing formu- ement for the implementation of general hyperelastic 3D-
lations available, elasticity at nite strains. Computer Methods in Applied
Mechanics and Engineering 1996;130:5779.
4. as a consequence of the nodal shear criterion, the
[7] Bischo M, Ramm E. Shear deformable shell elements
constraint count [14] has the ideal number of 1.5, for large strains and rotations. International Journal for
for any specic element type, Numerical Methods in Engineering 1997;40:442749.
5. all necessary mathematical operations, in particular [8] Bucalem ML, Bathe KJ. Higher-order MITC general
integration, are independent of the actual element shell elements. International Journal for Numerical
shape and can be performed a priori, Methods in Engineering 1993;36:372954.
334 K.-U. Bletzinger et al. / Computers and Structures 75 (2000) 321334

[9] Buchter N, Ramm E. 3D-extension of nonlinear shell [18] Parisch H. A continuum-based shell theory for nonlinear
equations based on the enhanced assumed strain concept. applications. International Journal for Numerical
In: Hirsch Ch, editor. Computational methods in applied Methods in Engineering 1993;38:185583.
sciences. Amsterdam: Elsevier, 1992. p. 5562. [20] Pinsky PM, Jang J. A C0-elastoplastic shell element
[10] Buchter N, Ramm E, Roehl D. Three-dimensional exten- based on assumed covariant strain interpolations. In:
sion of nonlinear shell formulation based on the Pande GN, Middleton J, editors. Proceedings of the
enhanced assumed strain concept. International Journal International Conference NUMETA 1987, Swansea,
for Numerical Methods in Engineering 1994;37:255168. 1987.
[11] Fox DD, Nagtegaal JC. A unied transverse shear treat- [21] Ramm E, Bischo M, Braun M. Higher order nonlinear
ment for 3-node and 4-node stress resultant shell el- shell formulations a step back into three dimensions.
ements. In: The 1997 US National Conference on In: Bell K, editor. From nite elements to the troll plat-
Computational Mechanics, San Francisco. 1997. form. Trondheim, Norway: Department of Structural
[12] Hauer C, Ramm E. Ecient shear deformable 3-node Engineering, Norwegian Institute of Technology, 1994.
plate/shell elements an almost hopeless undertaking. p. 6588.
In: Proceedings of the Third International Conference in [22] Sansour C. A theory and nite element formulation of
Computational Structures Technology CST '96, shells at nite deformations including thickness change:
Budapest, August 2123. 1996. circumventing the use of a rotation tensor. Archive of
[13] Hibbit, Karlsson, and Sorensen, Inc. ABAQUS Theory Applied Mechanics 1995;10:194216.
Manual. 1996. [23] Simo JC, Rifai S. A class of mixed assumed strain
[14] Hughes TJR. The nite element method linear static methods and the method of incompatible modes.
and dynamic nite element analysis. Englewood Clis, International Journal for Numerical Methods in
NJ: Prentice-Hall, 1987. Engineering 1990;29:1595638.
[15] Hughes TJR, Taylor RL. The linear triangular bending [24] Simo JC, Armero F. Geometrically non-linear enhanced
element. In: Whiteman JR, editor. The mathematics of strain mixed methods and the method of incompatible
nite elements and applications IV. New York: modes. International Journal for Numerical Methods in
Academic Press, 1981. Engineering 1992;33:141349.
[16] Hughes TJR, Tezduyar TE. Finite elements based upon [25] Szilard R. Theory and analysis of plates. Englewood
Mindlin plate theory with particular reference to the Clis, NJ: Prentice-Hall, 1974.
four-node isoparametric element. Journal of Applied [26] Simo JC, Hughes TJR. On the variational foundations of
Mechanics 1981;48:58796. assumed strain methods. Journal of Applied Mechanics
[17] Lyly M, Stenberg R, Vihinen T. A stable bilinear element 1986;53:514.
for the ReissnerMindlin plate model. Computer [27] Xu Z. A simple and ecient triangular nite element
Methods in Applied Mechanics and Engineering for plate bending. Acta Mechanica Sinica 1986;2:185
1993;110:34357. 92.