Efficient mixed Timoshenko-Mindlin shell elements

G.M.Kulikov* and S.V.Plotnikova

Department of Applied Mathematics and Mechanics,

Tambov State Technical University,

Sovetskaya Street, 106, Tambov, 392000 Russia

SUMMARY

The precise representation of rigid body motions in the displacement patterns of curved

Timoshenko-Mindlin (TM) shell elements is considered. This consideration requires the

development of the strain-displacement relationships of the TM shell theory with regard

to their consistency with the rigid body motions. For this purpose a refined TM theory of

multilayered anisotropic shells is elaborated. The effects of transverse shear deformation

and bending-extension coupling are included. The fundamental unknowns consist of

five displacements and eight strains of the face surfaces of the shell, and eight stress re-

sultants. On the basis of this theory the simple and efficient mixed models are devel-

oped. The elemental arrays are derived using the Hu-Washizu mixed variational princi-

ple. Numerical results are presented to demonstrate the high accuracy and effectiveness

of the developed 4-node shell elements and to compare their performance with other fi-

nite elements reported in the literature.

KEY WORDS: Timoshenko-Mindlin element; mixed model; rigid body motion; multi-

layered shell

*
Correspondence to: G.M.Kulikov, Department of Applied Mathematics and Mechanics, Tambov State
Technical University, Sovetskaya Street, 106, Tambov, 392000 Russia. E-mail: kulikov@apmath.tstu.ru,
Fax: (7 075) 253 2017, http://apm.tstu.ru/kulikov
 2

1. INTRODUCTION

One of the main requirements of the modern shell theory that is intended for the general

finite element (FE) formulation is that it must lead to strain-free modes for rigid body

motions. The adequate representation of rigid body motions is a necessary condition if

an element is to have good accuracy and convergence properties. Therefore, when an

inconsistent theory is used to construct any finite element, erroneous straining modes

under rigid body motions may be appeared. This problem has been studied for the

Kirchhoff-Love shell theory by Cantin [1] and Dawe [2].

Herein, the more general study on the basis of the refined Timoshenko-Mindlin

(TM) theory of multilayered shells is considered [3, 4]. The direct use of the traditional

TM shell theory (the first-order shear deformation theory) [5-9] for solving a series of

important shell problems such as the contact problems is not always convenient. In these

problems it is more convenient to select as unknown functions the displacements of the

top and bottom surfaces of the shell, since with the help of these displacements the ki-

nematic requirement of no penetration of the contact bodies can be fulfilled. Further-

more, the proposed TM shell theory can also simplify a formulation of new FE models

[10].

Using the classical linear TM shell theory in a FE formulation for plates and shells

is well established and has been shown to give acceptable results [11-18]. This theory

has the advantage that independent displacement and rotation trial functions may be

used and these functions need only to be C o continuous. The developed FE formulation,

based on the refined TM shell theory [3, 4] has the advantage because only independent

trial functions of displacements of the face surfaces may be used. It should be mentioned

that in some works (e.g. [19, 20]) developing the degenerate solid shell concept [11, 21]

displacement vectors of the bottom and top surfaces are also used and resolved with re-

spect to some global Cartesian basis in order to exactly describe rigid body motions.
 3

This allows, in particular special boundary conditions at the face surfaces of the shell to

be accounted for. The same idea of selecting as unknowns the displacements of the bot-

tom and top surfaces to construct any curved TM shell element is very attractive, since

only in this case we can deduce linear or non-linear TM strain-displacement relation-

ships that are free for all small or large rigid body motions, respectively. Taking into ac-

count that herein the displacement vectors of the face surfaces are represented in the lo-

cal reference surface basis, the developed FE formulation has substantial computational

advantages compared to the conventional isoparametric FE formulations, since it elimi-

nates the costly numerical integration by deriving the stiffness matrices. Indeed, our

element matrix requires only direct substitutions, no inversion is needed if the element

is rectangular, and it is evaluated by using the full exact analytical integration.

Our FE formulation is based on a simple and efficient approximation of shells via

quadrilateral 4-node elements developed by Hughes and Tezduyar [22], and by Wemp-

ner et al. [23]. The fundamental unknowns consist of five displacements and eight

strains of the face surfaces of the shell, and eight stress resultants. The simplest admissi-

ble approximations of the two-dimensional fields are used, namely, bilinear approxima-

tions of the displacements, and assumed approximations of the strains and stress resul-

tants. In this connection the element characteristic arrays are obtained by applying the

Hu-Washizu mixed variational principle. It is noteworthy that the stiffness matrix has

six, and only six, zero eigenvalues as required for satisfaction of the general rigid body

motion requirements.

Numerical results are presented to demonstrate the high accuracy and effectiveness

of the FE models developed and to compare their performance with other FE models

reported in the literature. For this purpose four tests were employed. They were pinched

cylinder tests, a shell roof test, and an open cylindrical multilayered composite shell test.
 4

2. STRAIN-DISPLACEMENT EQUATIONS

Let us consider the shell of uniform thickness h. The shell may be defined as a three-

dimensional body of volume V bounded by two bounding surfaces S − and S + , located

at distances δ − and δ + measured with respect to the reference surface S, and the edge

boundary surface Ω that is perpendicular to the reference surface (see Figure 1). Let the

reference surface S be referred to an orthogonal curvilinear coordinate system α1 and

α 2 which coincides with the lines of principal curvatures of its surface; e1 and e 2 de-

note the tangent unit vectors to the lines of principal curvatures. The α 3 axis is oriented

along the outward unit vector e 3 normal to the reference surface.

Fig. 1. Shell element.

The three-dimensional strain-displacement relationships for the general shell in a

vector form can be written as

1 ∂u 1 ∂u 1 ∂u
ε eii = ei , ε eij = ej + ei (i ≠ j)
H i ∂α i H i ∂αi H j ∂α j (1)
H α = A α (1 + k α α3 ) , H3 = 1
 5

where u = å u i e i is the displacement vector; u i (α1 , α 2 , α 3 ) are the components of this

i

vector; A α and k α are the Lamé coefficients and principal curvatures of the reference

surface; H α are the Lamé coefficients of any surface parallel to the reference surface.

Both here and in the following developments, unless otherwise specified, Greek indices

may take the values 1, 2 while Latin indices take the values 1, 2, 3.

The refined TM shell theory is based on the linear approximation of the displace-

ment vector in the thickness direction [3]:

u = N − (α3 )v − + N + (α3 )v + (2a)

v ± = å v α± e α + v 3e 3 (2b)
α

( )
N − (α3 ) = δ + − α 3 / h, ( )
N + (α3 ) = α3 − δ − / h (2c)

where v ± are the displacement vectors of the face surfaces S ± ; v α± (α1 , α 2 ) are the tan-

gential displacements of the face surfaces; v 3 (α1 , α 2 ) is the transverse displacement of

the reference surface; N − (α 3 ) and N + (α 3 ) are the linear shape functions. The linear

approximation (2) may be treated as a refined Timoshenko kinematic hypothesis (e.g.

works [7, 8], where as unknown functions the displacements and rotation components of

the reference surface are selected). The advantage of the proposed approach is obvious,

since with the help of the displacements v α± the special loading conditions at the face

surfaces of the shell can be formulated. Moreover, this simplifies a formulation of new

FE models and provides a convenient way to express the non-linear strain-displacement

relationships in terms of face surface strains [3, 10].

Substituting displacements (2a) into the strain-displacement relationships (1) and

taking into account formulas for the derivatives of the unit vectors e i along the coordi-
 6

nate lines α1 and α 2 [24], one can obtain the following strain-displacement equations

of the TM theory of the thick shells:

é 1 ∂v − 1 ∂v + ù
ε aγγ = ê N − (α 3 ) + N + (α 3 ) ú eγ, ε a33 = 0 (3a)
ëê H γ ∂α γ H γ ∂α γ ûú

é − 1 ∂v − 1 ∂v + ù é − 1 ∂v − 1 ∂v + ù
a
ε12 = ê N (α 3 ) + N (α 3 )
+
ú e 2 + ê N (α3 ) + N (α 3 )
+
ú e1
ë H1 ∂α1 H1 ∂α1 û ë H 2 ∂α 2 H 2 ∂α 2 û

ε aγ 3 =
Hγ
Hγ
βe γ +
1 ∂v
H γ ∂α γ
e3, β=
1 +
h
(
v − v− , ) v=
1 −
2
(
v + v+ ) (3b)

where H γ = A γ (1 + k γ δ ) are the Lamé coefficients of the middle surface;

( )
δ = δ − + δ + / 2 is the distance from the reference surface S to the middle surface of the

shell.

Replacing further the Lamé coefficients H γ by their values on the top and bottom

( )
surfaces H ±γ = A γ 1 + k γ δ ± in formulas (3a) for the tangential strains and by their val-

ues on the middle surface H γ in formulas (3b) for the transverse shear strains, the

strain-displacement equations of the refined TM theory of the moderately thick shells

are obtained

é − 1 ∂v − 1 ∂v + ù
ε bγγ = ê N (α 3 ) − + N (α 3 ) +
+
ú eγ ,
b
ε 33 =0 (4)
ëê H γ ∂α γ H γ ∂α γ ûú

é − 1 ∂v − 1 ∂v + ù é − 1 ∂v − 1 ∂v + ù
b
ε12 = ê N (α3 ) − + N (α3 ) +
+
ú e 2 + ê N (α3 ) − + N (α3 ) +
+
ú e1
ë H1 ∂α1 H1 ∂α1 û ë H 2 ∂α2 H 2 ∂α2 û

ε bγ 3 =βe γ +
1 ∂v
H γ ∂α γ
e3, β=
h
(
1 +
)
v − v− , v=
2
(
1 −
v + v+ )
The strain-displacement equations (4) are more attractive than equations (3) because

they are completely free for small rigid body motions. It will be shown in the next sec-
 7

b
tion. Besides, tangential strains ε αβ are dis-

tributed over the shell thickness according

to the linear law. As can be seen in Figure

2, it is an acceptable assumption for the

moderately thick shells because the cou-

pling conditions

( )
ε aαβ δ ± = ε αβ
b
( )
δ ± = E αβ
±

are valid.
Fig. 2. Distribution of tangential
More simple strain-displacement equa- strains over shell thickness.

tions can be obtained for the thin shells replacing the Lamé coefficients H γ and H γ by

the Lamé coefficients of the reference surface A γ in formulas (3a) and (3b). As a result

we have

é 1 ∂v − 1 ∂v + ù
ε cγγ = ê N − (α 3 ) + N + (α 3 ) ú eγ, ε c33 = 0 (5)
ëê A γ ∂α γ A γ ∂α γ ûú

é − 1 ∂v − 1 ∂v + ù é − 1 ∂v − 1 ∂v + ù
c
ε12 (
= ê N α3 ) +
(
+ N α3 ) (
ú e 2 + ê N α3 ) +
(
+ N α3 ) ú e1
ë A1 ∂α1 A1 ∂α1 û ë A 2 ∂α 2 A 2 ∂α 2 û

ε cγ 3 =βe γ +
1 ∂v
A γ ∂α γ
e3, β=
h
(
1 +
v − v− , ) v=
2
(
1 −
v + v+ )
As we shall see in the next section, the strain-displacement equations (5) can never be

free for small rigid body motions.

3. RIGID BODY MOTIONS

A small rigid body motion is defined as [24]

uR = ∆ + Φ × R (6)
 8

where ∆ = å ∆ i e i is the constant displacement (translation) vector; Φ = å Φ i e i is the

i i

constant rotation vector; R = r + α 3e 3 is the position vector of any point of the shell; r

is the position vector of any point of the reference surface (see Figure 1). In particular,

rigid body motions of the face surfaces are

v±R = ∆ + Φ × R ± (7)

where R ± = r + δ ± e 3 are the position vectors of points of the top and bottom surfaces.

The derivatives of the translation and rotation vectors with respect to the reference

surface coordinates are zero, i.e.,

∂∆ ∂Φ
= 0, =0 (8)
∂α γ ∂α γ

Taking into account the formulas for the derivatives of the unit vectors e i along the co-

ordinate lines [24] and using the equations (7) and (8), one can obtain the following ex-

pressions for the derivatives:

∂v ± R
= H ±γ Φ × e γ (9)
∂α γ

It can be verified by using equations (7) and (9) that the strains given by the equa-

tions (4) are all zero in a general rigid body motion, i.e.,
bR
ε αα = (Φ × e α ) e α = 0, ε ijbR = 0 (i ≠ j)
So, the TM theory of moderately thick shells is completely strain-free for all rigid body

motions.

Using again equations (7) and (9) in the strain-displacement equations (5), we get

the following results:

αα = (1 + k α α 3 )(Φ × e α ) e α = 0,
ε cR = (k1 − k 2 )α3 (Φ × e1 ) e 2
cR
ε12

α 3 = k α δ (Φ × e α ) e 3
ε cR
 9

showing that rigid body motions can never be completely strain-free for the TM theory

of thin shells. However, the transverse shear strains ε cR

α 3 = 0 when a reference surface is

selected to be the middle surface, since in this case δ = 0 . The tangential shear strain
cR
ε12 = 0 in the case of spherical shells or symmetrically loaded and supported isotropic

or orthotropic shells of revolution, and at the reference surface points of anisotropic

shells of arbitrary geometry.

4. HU-WASHIZU FUNCTIONAL

Let us consider the shell built up in the general case by the arbitrary superposition across

the wall thickness of N thin layers of uniform thickness hk. The kth layer may be defined

as a three-dimensional body of volume Vk bounded by two surfaces Sk-1 and Sk, located

at the distances δk-1 and δk measured with respect to the reference surface S, and the

edge boundary surface Ωk that is perpendicular to the reference surface (see Figure 3).

Fig. 3. Multilayered shell element.

 10

Here and in the following developments the index k = 1, N identifies the belonging of

any quantity to the kth layer. The full edge boundary surface Ω = Ω1 + Ω 2 + K + Ω N is

generated by the normals to the reference surface along the bounding curve Γ (with the

arc length s) of this surface. It is also assumed that the bounding surfaces Sk-1 and Sk are

continuous, sufficiently smooth and without any singularities. Let the reference surface

be referred to an orthogonal curvilinear coordinate system α1 and α2 which coincides

with the lines of principal curvatures of its surface. The α 3 axis is oriented along the

outward unit vector e 3 normal to the reference surface.

The constituent layers of the shell are supposed to be rigidly joined, so that no slip

on contact surfaces and no separation of layers can occur. The material of each constitu-

ent layer is assumed to be linearly elastic, anisotropic, homogeneous or fiber reinforced,

such that in each point there is a single surface of elastic symmetry parallel to the refer-

ence surface. Let p i− and p i+ be the intensities of the external loading acting on the bot-

tom surface S − = S 0 and top surface S + = S N in the αi coordinate directions, respec-

tively, while q ( k ) = q (νk ) ν + q (t k ) t + q (3k ) e 3 be the external loading vector acting on the

edge boundary surface Ω k . Here, q (νk ) , q (t k ) and q (3k ) are the components of its vector

acting in the ν, t and α 3 directions; ν and t are the normal and tangential unit vectors to

the bounding curve Γ.

The refined TM theory of multilayered shells is also based on the linear approxima-

tion of the displacement vector in the thickness direction (2), where we should set

δ − = δ0 and δ + = δ N . Substituting the displacements (2) and strains rewritten in the

more convenient form

ε αβ = N − (α 3 ) E αβ
−
+ N + (α3 ) E αβ
+
, ε α3 = E α3 , ε 33 = 0 (10)
 11

into the Hu-Washizu functional [3, 25] and taking into account the strain-displacement

relationships (4) for the moderately thick shell, one can obtain

ìï
J = òò í Π − å R αβ
− −
E αβ[ (
− e −αβ + R αβ
+ +
E αβ ) +
− e αβ ( )]
S ïî α≤β

[ ( ) ( ) ] ( ü
− å R α3 (E α3 − e α3 ) + Q α− − p −α v −α + Q α+ + p +α v +α − Q 3 − p 3− + p 3+ v 3 ý dS (11) )
α þ
( )
− ò R̂ −νν v −ν + R̂ +νν v +ν + R̂ −νt v −t + R̂ +νt v +t + R̂ ν 3 v 3 ds
Γ

where

1 æ 1 ∂v ±γ ö
e ±γγ = ± ç + Bδ v δ± + k γ v 3 ÷
ζγ ç A γ ∂α γ ÷
è ø

± 1 æ 1 ∂v 2± ±ö 1 æ 1 ∂v 1± ö 1
e12 = ç
ç A ∂α − B v ÷
2 1 ÷ + ç
± ç
− B1v 2± ÷÷ , e γ3 = β γ − θγ
ζ1± è 1 1 ø ζ 2 è A 2 ∂α 2 ø ζγ
(12)
θγ = k γ v γ −
1 ∂v 3
A γ ∂α γ
, βγ =
1 +
h
(
v γ − v −γ , ) vγ =
2
(
1 −
v γ + v +γ )
1 ∂A δ
ζ ±γ = 1 + k γ δ ± , ζ γ = 1 + k γ δ, Bγ = (δ ≠ γ )
A1A 2 ∂α γ

Here, Π is the strain energy density; v ±ν , v ±t and v 3 are the components of the dis-

placement vectors of the face surfaces in the coordinate system ν , t and α 3 (see Figure

±
3); E αβ are the tangential strains of the bottom and top surfaces; E α3 are the transverse

±
shear strains of the middle surface; R αβ and R α3 are the generalized and classical stress

resultants; Q α± and Q 3 are the generalized and classical body force resultants; R̂ ±νν , R̂ ±νt

and R̂ ν 3 are the generalized and classical external load resultants, which are defined as
 12

Π=
1
å å
2 α≤β γ ≤δ
[
A 00 − − 01 − + +
( − 11 + +
αβγδ E αβ E γδ + A αβγδ E αβ E γδ + E αβ E γδ + A αβγδ E αβ E γδ ) ]
δk
1
+ å A α3 γ 3 E α3 E γ 3 , A α3 γ 3 = å ò Cα3γ 3dα3
(k)
(13)
2 α ,γ k δ k −1

δk

ò Cαβγδ [N (α3 )] [N (α3 )] dα3

2− m − n m+n
mn
A αβγδ =å (k) − +
(m, n = 0, 1)
k δk −1

and
δk δk
±
R αβ =å ò
k) ±
σ (αβ N (α3 )dα3 , R α3 = å ò σ α3 dα3
(k)
(14)
k δ k −1 k δ k −1

δk δk
Q ±α = å ò f α N (α3 )dα3 , Q3 = å
(k) ±
ò f3
(k)
dα 3
k δ k −1 k δk −1

δk δk
R̂ ±νν =å ò q (νk ) N ± (α 3 )dα 3 , R̂ ±νt =å ò qt
(k)
N ± (α 3 )dα3
k δk −1 k δ k −1

δk
R̂ ν 3 = å ò q3
(k)
dα 3
k δ k −1

In formulas (14) fα( k ) and f3( k ) are the externally applied body forces of the kth layer,

while σ (αβ
k)
and σ (αk3) denote the tangential and transverse shear stresses of the kth layer

that can be found as

σ (αβ
k)
= å C (αβγδ
k)
ε γδ , σ (αk3) = å C α( k3)γ 3ε γ 3
γ ≤δ γ

where C (αβγδ
k)
and C (αk3)γ 3 are the stiffness coefficients of the kth layer.

5. FE FORMULATION

It is well-known that the Hu-Washizu variational principle provides the basis for the

derivation of various variational principles, and many different mixed and hybrid finite
 13

elements can be designed [26]. Herein, the Hu-Washizu functional (11)-(14) for the ele-

ment can be written in the following form:

1 1
J el =
é1
ò ò êë 2 E
T
( ) ù
AE − ET − v T B T R − v T (P + Q )úΛdξ1dξ 2 − ò v TΓ R
ˆ Γ ds (15)
−1 −1 û Γel

where ξ1 and ξ 2 are the local coordinates of the element that vary from -1 to +1;

Λ(ξ1 , ξ 2 ) is the function characterizing the metric of the element;

[
v = v1− v1+ v 2− v 2+ v 3 ]
T
[
is the displacement vector; v Γ = v −ν v +ν v −t v +t v 3 ]
T
is the

−
displacement vector of the element edge Γ el ; E = E11 +
E11 −
E 22 +
E 22 −
E12 +
E12 [
E13 E 23 ]T

−
is the strain vector; R = R 11 +
R 11 −
R 22 [+
R 22 −
R 12 +
R 12 R 13 R 23 ]
T
is the stress resultant

[
ˆ Γ = R̂ −νν R̂ +νν R̂ −νt R̂ +νt R ν3
vector; R ]T
is the loading resultant vector acting on the ele-

[
ment edge Γ el ; Q = Q1− Q1+ Q 2− Q 2+ Q 3 ]
T
is the body force resultant vector;

[
P = − p1− p1+ − p 2− p 2+ − p 3− + p 3+ ] T
is the surface traction vector; A is the constitutive

stiffness matrix; B is the strain-displacement matrix.

For the quadrilateral 4-node shell element the displacement field is approximated

according to the standard C o interpolation

v = å Nlvl (16)
l

[
where v l = v1−l v1+l v 2−l v 2+l v 3l ]T
are the displacement vectors of the element nodes;

N l (ξ1 , ξ 2 ) are the linear shape functions of the element; l = 1, 4 .

In accordance with [22, 23], the 20 modes of this element are the six rigid body mo-
−00 +00 00
tions, the eight homogeneous states of strains E αβ , E αβ , E α3 , and the six additional

−01 +01
modes representing higher approximations of tangential normal strains E11 , E11 ,
 14

−10 +10 01
E 22 , E 22 and transverse shear strains E13 , E10
23 , i.e., we have following strain interpo-

lations:

E = E 00 + E10 ξ1 + E 01ξ2

[ −00
E 00 = E11 + 00
E11 −00
E 22 +00
E 22 −00
E12 + 00
E12 00
E13 E 00
23 ]
T
(17)

[−01
E 01 = E11 + 01
E11 01
0 0 0 0 E13 ]T
[
−10
0 , E10 = 0 0 E 22 +10
E 22 0 0 0 E10
23 ]T

The interpolations of the stress resultants follow the forms of the conjugate strains

R = R 00 + R 10 ξ1 + R 01ξ2 (18)

where the vectors R 00 , R 01 and R 10 are defined from equations (17) by replacing a

letter E by a letter R.

The governing equations for the element are obtained by applying the Hu-Washizu

variational principle (15). Using the equations (16)-(18), and eliminating the strain and

stress resultant parameters on the element level, one can obtain

Ku = F

where K is the elemental stiffness matrix; F is the load vector; u is the vector of five

displacement components at nodal points of the element.

It should be noted that the formulation of the stiffness matrix K requires only direct

substitutions; no inversion is needed if the element is rectangular. The matrix is sym-

metric and positive definite, and has six, and only six, zero eigenvalues as required for

satisfaction of the general rigid body motion representation. Furthermore, the element

matrix is evaluated by using the full exact analytical integration and the element does

not contain any spurious zero energy modes. So, our FE formulation especially for the

rectangular elements is very economical and efficient.

6. NUMERICAL TESTS

Four tests were employed to assess the effectiveness of the developed TMS4 element

based on the strain-displacement equations (4) that are completely free for all rigid body
 15

motions. They were a pinched cylinder with rigid end-diaphragms, a pinched cylinder

with free edges, a cylindrical shell roof, and a multilayered angle-ply cylindrical shell.

6.1. Pinched cylinder with rigid diaphragms

To illustrate the capability of the developed TMS4 element to overcome membrane

and shear locking phenomenon and to compare it with the different 4-node quadrilateral

elements [15, 16, 18, 27, 28], we consider

one of the most demanding standard linear

test [29]. A thin cylinder supported by rigid

end-diaphragms is loaded by two opposite

concentrated forces in its middle section.

The geometrical and material data of the

problem are shown in Figure 4(a).

Three types of boundary conditions

Fig. 4. Pinched cylinder under
modeling the rigid end-diaphragm can be opposite radial forces with: (a) rigid
used diaphragms and (b) free edges.

v −2 = v 3 = 0 (19a)

v +2 = v 3 = 0 (19b)

v −2 = v 2+ = v 3 = 0 (19c)

Besides, due to symmetry of the problem, only one octant of the cylinder is modeled

with a regular mesh of TMS4 elements.

Table 1 lists a comparison of the normalized radial displacement under the applied

load between the TMS4 element and aforementioned 4-node quadrilateral elements

[30]. The displacements are normalized with respect to the analytical solution

− 1.8245 × 10 −5 [29]. Note that all the three variants of boundary conditions (19) lead to

the practically identical results and very well model the ‘real’ boundary conditions
 16

v2 = v3 = 0 (20)

used by Heppler and Hansen [29]. As can be seen from Table 1, our results show an ex-

cellent agreement even for coarse meshes.

Table 1. Pinched cylinder with rigid diaphragms. Normalized radial displacement under
the concentrated load.
Mesh MITC4 RSDS Mixed SRI QPH TMS4 TMS4
[15] [16] [18] [27] [28] (19a) (19c)
4x4 0.370 0.469 0.399 0.373 0.370 0.890 0.890

8x8 0.740 0.791 0.763 0.747 0.740 0.941 0.941

16x16 0.930 0.946 0.935 0.935 0.930 0.986 0.986

6.2. Pinched cylinder with free edges

The pinched cylinder with free edges has been also extensively treated for numeri-

cal testing of new FE models. The geometrical and material properties of the cylinder

are shown in Figure 4(b). Due to symmetry of the problem, only one octant of the cylin-

der is discretized. The radial deflection at the load is given in Figure 5. The curves

marked by — display the results obtained by using the TMS4 element. Additionally, in

Figure 5 the solutions, based on the traditional cubic Lagrange element [31] (see curves

marked by ®) and the cubic Lagrange element [31] that has been augmented in the man-

ner of Cantin [32] to provide the correct rigid body motions (see curves marked by p),

are presented. Both the Lagrange elements are based on the classical TM shell theory.

Note that the augmented Lagrange element solution converges to value − 0.1114 [31]

as do the developed TMS4 element and the cubic Lagrange element. However, our sim-

ple and efficient element solution converges more rapidly, since the cubic Lagrange

element possesses only two zero eigenvalues. At the same time its augmented formula-

tion has the six zero eigenvalues as required for satisfaction of the general rigid body

motion requirements.
 17

Fig. 5. Pinched cylinder with free edges. Radial displacement under load.

6.3. Cylindrical shell roof

Let us consider an open cylindrical shell segment subjected to gravitational self-

weight loading and supported at its curved edges by rigid diaphragms, while the straight

edges are free [33]. This problem has become a de facto standard test and has been fre-

quently used for numerical testing of FE approximations. The geometrical and material

characteristics of the shell are depicted in Figure 6. Due to symmetry of the problem,

only one quarter of the roof is modeled with a regular mesh of TMS4 elements. Figure 7

shows the distribution of the axial displacement at the diaphragm (around the curve AD)

and the vertical displacement at the middle span (around the curve BC), while Figure 8

( +
displays the distribution of the moment resultants M αα = h R αα −
− R αα)/ 2 at the middle

span. The solid curves show the results obtained by using the exact solution [33] as re-

ported in works [11, 12]. One can see that the moment resultant M 11 is calculated with

an excellent exactitude for both meshes. However, the results obtained for the moment

resultant M 22 are less satisfactory. This phenomenon is explained by using the constant
 18

±
approximation over the element for generalized stress resultants R 22 in ϕ direction in

accordance with equations (18).

Fig. 6. Cylindrical shell roof under self-weight loading.

(a) (b)

Fig. 7. Cylindrical shell roof. Displacements at: (a) diaphragm and (b) middle span.

Table 2 additionally presents a comparison of the normalized vertical displacement

at the point C of the middle span between the TMS4 element and aforementioned 4-

node quadrilateral elements [15, 18, 27, 28]. The displacements [30] are normalized to

the value − 0.3024 calculated by MacNeal and Harder [34], while our displacement is

normalized to the value − 0.3015 . Such a value is a computationally exact solution of

 19

this problem based on the developed TM

shell theory. The slightly small differ-

ence between the reference solutions is

explained by a simple fact: in this test

the shell response is more sensitive to the

types of boundary conditions modeling

the rigid end-diaphragm (see for the

comparison a pinched cylinder with rigid

diaphragms). One can observe that Fig. 8. Cylindrical shell roof. Moment
resultants at middle span.
boundary conditions (19c) used in our

FE formulation provide the more rigid end support than (20) used in works [33, 34].

Table 2. Cylindrical shell roof. Normalized vertical displacement at point C (see Figure 6).
Mesh MITC4 [15] Mixed [18] SRI [27] QPH [28] TMS4

4x4 0.940 1.083 0.964 0.940 0.864

8x8 0.970 1.015 0.984 0.980 0.962

16x16 1.000 1.000 0.999 1.010 0.989

6.4. Multilayered angle-ply cylindrical shell

It is apparent that using the simplified strain-displacement equations (5) can lead to

incorrect results for the moderately thick composite shells. To assess this statement, we

consider an open cylindrical two-layer angle-ply shell rigidly clamped at its curved

edges and supported at straight edges by rigid diaphragms. The shell is subjected to the

uniform stretching v 0 as shown in Figure 9. The material characteristics of each layer

were taken to be those typical of a high modulus composite [7] and are given in Figure

9, where subscripts L and T refer to the longitudinal and transverse directions of the in-

dividual ply. Let the ply thicknesses and ply orientations be [ h / 2, h / 2] and [− γ , + γ ] ,
 20

where γ is measured in the clockwise direction from x to the fiber direction. Due to the

anisotropic shell response, we did not adopt symmetry conditions and modeled the

whole shell by using regular meshes of TMS4 and TMS4c elements. A new developed

TMS4c element is based on the simplified strain-displacement equations (5).

Fig. 9. Two-layer angle-ply cylindrical shell subjected to uniform stretching.

In Table 3 the values of the dimensionless central transverse displacement of the

shell segment for both TMS4 and TMS4c elements and for various ply orientations and

meshes are presented. It can be seen that using the strain-displacement equations (4)

only insignificantly updates the results in a comparison with less general strain-

displacement equations (5). Thus, it is possible to recommend the geometrically linear

TM shell theory on the basis of strain-displacement equations (5) for using in engineer-

ing calculations. At the same time applying the TM shell theory, which is not com-

pletely strain-free for rigid body motions for solving the composite shells undergoing

large deflections and large rotations can lead to significant errors. Such a problem is cur-

rently under development. Finally, Figure 10 shows the distribution of the longitudinal

strains of the face surfaces at x = 45 mm in ϕ direction for various ply orientations. It is

 21

seen that a response of the shell is very unusual to the ply orientation − 30 o , 30 o , [ ]
where the longitudinal strain of the bottom surface is negative.

Table 3. Central transverse displacement − v 3 ( 0, 0) / v 0 of the two-layer angle-ply cy-

lindrical shell.
Mesh TMS4 TMS4c

[0o, 90o] [-15o, 15o] [-30o, 30o] [-60o, 60o] [0o, 90o] [-15o, 15o] [-30o, 30o] [-60o, 60o]

4x4 0.077 1.705 4.091 1.469 0.077 1.703 4.085 1.468

4x8 0.068 1.346 3.359 1.259 0.067 1.344 3.352 1.259

4x12 0.066 1.330 3.301 1.219 0.066 1.329 3.295 1.219

8x8 0.065 1.439 3.501 1.182 0.064 1.437 3.493 1.181

8x16 0.063 1.413 3.418 1.126 0.062 1.412 3.411 1.126

8x24 0.063 1.409 3.404 1.116 0.062 1.407 3.398 1.115

12x12 0.063 1.434 3.456 1.131 0.062 1.433 3.449 1.131

12x24 0.062 1.423 3.422 1.106 0.062 1.422 3.416 1.105

(a) (b)

−
Fig. 10. Distribution of longitudinal strains of bottom and top surfaces: (a) E11 and
+
(b) E11 at x = 45 mm in ϕ direction
 22

7. CONCLUSION

The simple and efficient mixed models have been developed for the analysis of multi-

layered anisotropic TM shells. The first FE formulation is based on the strain-

displacement equations of the moderately thick shell, which are completely free for rigid

body motions. The second FE formulation is based on the strain-displacement equations

of the thin shell that cannot be completely free for rigid body motions in the case of ani-

sotropic shells of arbitrary geometry. As fundamental unknowns five displacements and

eight strains of the face surfaces of the shell, and eight stress resultants have been cho-

sen. This allows, in particular special loading conditions at the bottom and top surfaces

of the shell to be accounted for.

The elemental stiffness matrices of our FE formulations are symmetric and positive

definite and have six zero eigenvalues as required for satisfaction of the general rigid

body motion representation. Besides, the elemental matrices require only direct substitu-

tions (no inversion is needed) if elements are rectangular and they are evaluated using

the full exact analytical integration. It is important that the developed TMS4 and TMS4c

elements do not contain any spurious zero energy modes.

To demonstrate the high accuracy and effectiveness of the developed elements four

tests were employed. They were pinch tests, an open cylindrical shell roof test and an

open cylindrical composite shell test.

The extension to finite deflections poses no additional difficulties but requires alge-

bra and computation efforts. For this purpose the Hu-Washizu mixed functional for an

analysis of geometrically non-linear multilayered anisotropic shells [3, 25] can be used.

The extension to initially stressed multilayered shells can be also done [35].

ACKNOWLEDGEMENTS

The present research was supported by Russian Fund of Basic Research (Grant No. 98-

01-04076) and by Research Programme of Tambov State Technical University (Grant

No. 1Γ/00-10).
 23

REFERENCES

1. Cantin G. Strain displacement relationships for cylindrical shells. AIAA Journal

1968; 6: 1787-1788.

2. Dawe DJ. Rigid-body motions and strain-displacement equations of curved shell

finite elements. International Journal of Mechanical Sciences 1972; 14: 569-578.

3. Kulikov GM, Plotnikova SV. Comparative analysis of two algorithms for numerical

solution of nonlinear static problems for multilayered anisotropic shells of revolu-

tion. 1. Account of transverse shear. Mechanics of Composite Materials 1999; 35:

241-248.

4. Kulikov GM. Refined global approximation theory of multilayered plates and

shells. Journal of Engineering Mechanics 2001; 127:119-125.

5. Timoshenko SP. On the correction for shear of the differential equation for trans-

verse vibrations of prismatic bars. Philosophical Magazine and Journal of Science.

Ser. 6 1921; 41: 744-746.

6. Mindlin RD. Influence of rotatory inertia and shear on flexural motions of isotropic

elastic plates. Trans. ASME, Journal of Applied Mechanics 1951; 18: 31-38.

7. Grigolyuk EI, Kulikov GM. Multilayered Reinforced Shells: Analysis of Pneumatic

Tires. Mashinostroyenie: Moscow, 1988 (in Russian).

8. Noor AK, Burton WS. Assessment of computational models for multilayered com-

posite shells. Applied Mechanics Reviews. 1990; 43: 67-97.

9. Reddy JN. Mechanics of Laminated Composite Plates: Theory and Analysis. CRC

Press: Boca Raton, 1997.

10. Kulikov GM, Plotnikova SV. Finite element formulation of straight composite

beams undergoing finite rotations. Trans. Tambov State Technical University 2001;

7: 617-633.
 24

11. Ahmad S, Irons BM, Zienkiewicz OC. Analysis of thick and thin shell structures by

curved finite elements. International Journal for Numerical Methods in Engineer-

ing 1970; 2: 419-451.

12. Zienkiewicz OC, Taylor RL, Too JM. Reduced integration technique in general

analysis of plates and shells. International Journal for Numerical Methods in Engi-

neering 1971; 3: 275-290.

13. Hughes TJR, Taylor RL, Kanok-Nukulchai W. A simple and efficient finite element

for plate bending. International Journal for Numerical Methods in Engineering

1977; 11: 1529-1543.

14. Malkus DS, Hughes TJR. Mixed finite element methods - reduced and selective in-

tegration techniques: a unification of concepts. Computer Methods in Applied Me-

chanics and Engineering 1978; 15: 63-81.

15. Bathe KJ, Dvorkin EN. A formulation of general shell elements - the use of mixed

interpolation of tensorial components. International Journal for Numerical Methods

in Engineering 1986; 22: 697-722.

16. Lam D, Liu WK, Law ES, Belytschko T. Resultant-stress degenerated-shell ele-

ment. Computer Methods in Applied Mechanics and Engineering 1986; 55: 259-

300.

17. Stolarski H, Belytschko T. On the equivalence of mode decomposition and mixed

finite elements based on the Hellinger-Reissner principle. Part I: Theory. Part II:

Applications. Computer Methods in Applied Mechanics and Engineering 1986; 58:

249-284.

18. Simo JC, Fox DD, Rifai MS. On a stress resultant geometrically exact shell model.

Part II: The linear theory; computational aspects. Computer Methods in Applied Me-

chanics and Engineering 1989; 73: 53-92.

 25

19. Park HC, Cho C, Lee SW. An efficient assumed strain element model with six DOF

per node for geometrically non-linear shell. International Journal for Numerical

Methods in Engineering 1995; 38: 4101-4122.

20. Klinkel S, Gruttmann F, Wagner W. A continuum based three-dimensional shell

element for laminated structures. Computers and Structures 1999; 71: 43-62.

21. Bathe KJ. Finite Element Procedures. Prentice Hall: New Jersey, 1996.

22. Hughes TJR, Tezduyar TE. Finite elements based upon Mindlin plate theory with

particular reference to the four-node bilinear isoparametric element. Trans. ASME,

Journal of Applied Mechanics 1981; 48: 587-596.

23. Wempner G, Talaslidis D, Hwang CM. A simple and efficient approximation of

shells via finite quadrilateral elements. Trans. ASME, Journal of Applied Mechanics

1982; 49: 115-120.

24. Gol'denveiser AL. Theory of Elastic Thin Shell. Pergamon Press: Oxford, 1961.

25. Kulikov GM. Variational equation for the nonlinear multilayered anisotropic shell

of variable stiffness. Trans. Tambov State Technical University 1997; 3: 119-125.

26. Pian THH, Tong P. Basis of finite element methods for solid continua. International

Journal for Numerical Methods in Engineering 1969; 1: 3-28.

27. Hughes TJR, Liu WK. Nonlinear finite element analysis of shells. Part II: Two-

dimensional shells. Computer Methods in Applied Mechanics and Engineering

1981; 27: 167-182.

28. Belytschko T, Leviathan I. Physical stabilization of the 4-node shell element with

one-point quadrature. Computer Methods in Applied Mechanics and Engineering

1994; 113: 321-350.

29. Heppler GR, Hansen JS. A Mindlin element for thick and deep shells. Computer

Methods in Applied Mechanics and Engineering 1986; 54: 21-47.

 26

30. Argyris JH, Papadrakakis M, Apostolopoulou C, Koutsourelakis S. The TRIC shell

element: theoretical and numerical investigation. Computer Methods in Applied Me-

chanics and Engineering 2000; 182: 217-245.

31. Hansen JS, Heppler GR. A Mindlin shell element that satisfies rigid-body require-

ments. AIAA Journal 1985; 23: 288-295.

32. Cantin G. Rigid body motions in curved finite elements. AIAA Journal 1970; 8:

1252-1255.

33. Scordelis AC, Lo KS. Computer analysis of cylindrical shells. American Concrete

Institute Journal 1964; 61: 539-560.

34. MacNeal RH, Harder RL. A proposed standard set of problems to test finite element

accuracy. Finite Elements in Analysis and Design 1985; 1: 3-20.

35. Kulikov GM. Analysis of initially stressed multilayered shells. International Jour-

nal of Solids and Structures 2001; 38: 4535-4555.