17.

MODELING OF CURVED AND DOUBLY-

CURVED SHELLS BY FINITE ELEMENT METHOD
BASED SOFTWARE SYSTEMS

Szerző: Dr. András Szekrényes

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

2 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems

17. MODELING OF CURVED AND DOUBLY-CURVED SHELLS BY

FINITE ELEMENT METHOD BASED SOFTWARE SYSTEMS

17.1 Curved shell elements

Curved shell elements are suitable to model the midsurface geometry more accurately. In
the case of certain surfaces – for example the cylindrical shell – it means the exact descrip-
tion of the original surface. For more complicated cases – similarly to the displacement
field - the curvatures of the surface are approximated by interpolation functions. In this
respect such elements belong to the parametric element types [1].

17.2 Thin-walled cylindrical shell element

The thin cylindrical shell element is presented in Fig.17.1. In accordance with the basic
equations of the technical theory of thin shells the geometrical properties of the cylindrical
shell are the followings [1,2]:
q1 = x , H 1 = 1 , R1 = ∞ , (17.1)
q 2 = ϕ , H 2 = R , R2 = R .

Fig.17.1. Parameters of the thin cylindrical shell element.

Apart from the displacement components u, v and w we can derive the angle of rotations
by applying the basic equations of the technical theory of thin shells based on Eqs.(14.69)
and (14.70):
u 1
β1 = β x = − w,1 = − w, x , (17.2)
R1 H 1
v 1 1
β 2 = βϕ = − w, 2 = (v − w,ϕ ) ,
R2 H 2 R

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 3

1
β3 = ( Rv, x − u ,ϕ ) .
2R
Next, we calculate the strain components using the parameters of the cylindrical shell sur-
face (see Eq.(14.67)):
1 H 1, 2 1
ε 11 = ε x = u ,1 + v+ w = u,x , (17.3)
H1 H1H 2 R1
1 H 2,1 1 1
ε 22 = ε ϕ = v, 2 + u+ w = (u ,ϕ + w) ,
H2 H1H 2 R2 R
H1  u  H  v  1
2γ 12 = 2γ xϕ =   + 2   = u ,ϕ + v, x ,
H2  H 1  , 2 H 1  H 2  ,1 R
1 H
κ 11 = κ x = β1,1 + 1, 2 β 2 = w, xx ,
H1 H1H 2
1 H 1
κ 22 = κϕ = β 2, 2 + 2,1 β1 = 2 (v,ϕ − w,ϕϕ ) ,
H2 H1H 2 R
H1  β1  H β   1 1 1
2κ12 = 2κ xϕ =   + 2  2  +  −  β 3 = (−2 w, xϕ + v, x + β 3 ) .
H 2  H 2 , 2 H1  H1 ,1  R2 R1  R
The rigid body-like motion of the element involves six degrees of freedom, which are giv-
en by the displacement vector field given below [1]:
u 0 = a1 + a 2 R(cos ϕ − cos θ ) − a3 R sin ϕ , (17.4)
v0 = −a 2 x sin ϕ + a3 x cos ϕ + a 4 R (cos ϕ cos θ − 1) − a5 sin ϕ + a 6 cos ϕ ,
w0 = a2 x cos ϕ + a3 x sin ϕ + a4 R sin ϕ cos ϕ + a5 cos ϕ + a6 sin ϕ .
In matrix form:
u 0 = Φ0α 0 , (17.5)
where:
u 0 = [u 0 w0 ] ,
T
v0 (17.6)
1 R (cos ϕ − cos θ ) − R sin ϕ 0 0 0 
Φ 0 = 0 − x sin ϕ x cos ϕ R (cos ϕ cos θ − 1) − sin ϕ cos ϕ  ,
0 x cos ϕ x sin ϕ R sin ϕ cos ϕ cos ϕ sin ϕ 
α T0 = [a1 a 2
a3 a 4 a5 a6 ] .
where Φ 0 is the matrix of interpolation functions, which capture the rigid body-like mo-
tions, α0 is the vector of unknown coefficients. The displacement field of rigid body-like
motion and that of the deformation together give the total displacement field, which is:
u 
u = u 0 + u 01 =  v  , (17.7)
 w
where:
u 01 = [u 01 v01 w01 ],
T
(17.8)
u 01 = a 7 x + a8ϕ + a9 xϕ ,

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

4 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems
v01 = a10ϕ + a11 xϕ ,
w01 = a12 x + a13 xϕ + a14ϕ 2 + a15 x 3 + a16 x 2ϕ + a17 xϕ 2 + a18ϕ 3 +
2

+ a19 x 3ϕ + a20 x 2ϕ 2 + a21 xϕ 3 + a22 x 3ϕ 2 + a23 x 2ϕ 3 + a24 x 3ϕ 3 .

In matrix form we have:
u = u 0 + u 01 = Φ 0 α 0 + Φ 1 α 1 , (17.9)
where Φ 1 is the interpolation functions matrix related to the deformation displacement field:
x ϕ xϕ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
Φ1 = 0 0 0 ϕ xϕ 0 0 0 0 0 0 0 0 0 0 0 0 0 .
0 0 0 0 0 x2 xϕ ϕ 2 x3 x 2ϕ xϕ 2 ϕ 3 x 3ϕ x 2ϕ 2 xϕ 3 x 3ϕ 2 x 2ϕ 3 x 3ϕ 3 
(17.10)
α1 is the vector of unknown coefficients:
α 1T = [a7 a8 a9 a10 a11 a12 a13 a14 a15 .........
(17.11)
a16 a17 a18 a19 a20 a21 a22 a23 a24 ].
In the expression above there are 24 unknown coefficients. To determine all of them 24
displacement parameters are required, let us choose the followings:
T
[
u~ e = u~1 v~1 w ~ β~
1 x1
~
βϕ1 w ~
, xϕ 1 u~2 v~2 w~ β~
2 x2
~
βϕ 2 w~ .....
, xϕ 2

~ ~ ~ ~ ~ ~ ~ ~ ~
.........u3 v3 w3 β x 3 β ϕ 3 w, xϕ 3 u 4 v4 w4 β x 4 β ϕ 4 w
~ ~ ~ ,
, xϕ 4 ]
(17.12)
i.e. at each node there are displacements in the direction of the basis vectors, there are rota-
tions about e1 and e2, the sixth degrees of freedom is chosen to be the mixed derivative
w,xϕ. The degrees of freedom for the cylindrical shell element based on Eqs.(17.2) and
(17.7) are:
u = u 0 + u 01 , v = v0 + v01 , w = w0 + w01 , (17.13)
β x = − w, x = −a2 cos ϕ − a3 sin ϕ − 2a12 x − a13ϕ − 3a15 x 2 − 2a16 xϕ − a17ϕ 2 +
− 3a19 x 2ϕ − 2a20 xϕ 2 − a21ϕ 3 − 3a22 x 2ϕ 2 − 2a23 xϕ 3 − 3a24 x 2ϕ 3 ,
1 1
βϕ = (v − w,ϕ ) = (−a4 R + a10ϕ + a11 xϕ − a13 x − 2a14ϕ − a16 x 2 − 2a17 xϕ − 3a18ϕ 2 +
R R
− a19 x − 2a 20 x 2ϕ − 3a21 xϕ 2 − 2a22 x 3ϕ − 3a23 x 2ϕ 2 − 3a24 x 3ϕ 2 ),
3

w, xϕ = −a2 sin ϕ + a3 cos ϕ + a13 + 2a16 x + 2a17ϕ + 3a19 x 2 +

+ 4a20 xϕ + 3a21ϕ 2 + 6a22 x 2ϕ + 6a23 xϕ 2 + 9a24 x 2ϕ 2 .
The conditions for the determination of the parameters ai, i = 1...24 are:
u ( L / 2,−θ ) = u~1 , u ( L / 2, θ ) = u~2 , (17.14)
~ ~
u (− L / 2, θ ) = u 3 , u (− L / 2,−θ ) = u 4 ,
v( L / 2,−θ ) = v~1 , v( L / 2, θ ) = v~2 ,
v(− L / 2, θ ) = v~3 , v(− L / 2,−θ ) = v~4 ,
w( L / 2,−θ ) = w ~ , w( L / 2, θ ) = w ~ ,
1 2
~
w(− L / 2, θ ) = w3 , w(− L / 2,−θ ) = w ~ ,
4
~ ~
β x ( L / 2,−θ ) = β x1 , β x ( L / 2,θ ) = β x 2 ,
~ ~
β x (− L / 2,θ ) = β x 3 , β x (− L / 2,−θ ) = β x 4 ,

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 5

~ ~
βϕ ( L / 2,−θ ) = βϕ1 , βϕ ( L / 2,θ ) = βϕ 2 ,
~ ~
βϕ (−L / 2,θ ) = βϕ 3 , βϕ (−L / 2,−θ ) = βϕ 4 ,
~ , w ( L / 2, θ ) = w
w, xϕ ( L / 2,−θ ) = w ~ ,
, xϕ 1 , xϕ , xϕ 2
~ , w (− L / 2,−θ ) = w
w, xϕ (− L / 2, θ ) = w ~ .
, xϕ 3 , xϕ , xϕ 4

The vector of nodal displacements is formulated similarly to the plate element presented in
section 15, viz. [3]:
u~ e = M A , (17.15)
where vector A contains the elements of α 0 and α 1 in order, i.e.:
A = [a1
T
a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 .........
.......... a13 a14 a15 a16 a17 a18 a 23 a 24 ].
a19 a 20 a 21 a 22
(17.16)
The coefficients of the interpolation polynomials are determined by the inversion of M :
−1
A = M u~ . e (17.17)
Due to their large length, the coefficients are not detailed here. Assuming that the angle
rotation, β3 is approximately zero the strain components are formulated as follows:
1 1
ε x = u , x , ε ϕ = (u ,ϕ + w) , 2γ xϕ = u ,ϕ + v, x , (17.18)
R R
1 1
κ x = w, xx , κϕ = 2 (v,ϕ − w,ϕϕ ) , 2κ xϕ = (−2w, xϕ + v, x ) .
R R
The strain components can be classified into two parts: strains, shear strains and the curva-
tures, respectively. We can write that:
 εx   κx 
   
ε 0 =  ε y  = R 0 A , κ =  κ y  = R1 A , (17.19)
2γ xy  2κ xy 
   
where matrices R 0 and R 1 are calculated based on the derivatives of the displacement func-
tions:

0 0 0 0 0 0 1 0 ϕ 0 0 0
 − sin ϕ +  − cos ϕ + 
     cos ϕ sin ϕ 1 x x2
R 0 = 0  x   x  sin ϕ cosθ 0 0 0 ....
 + R cos ϕ  + R sin ϕ  R R R R R
 3 1 1
0 x
− sin ϕ cos ϕ 0 0 0 0 0 ϕ 0
 2 2 2R 2R
0 0 0 0 0 0 0 0 0 0 0 0 
xϕ ϕ 2 x 3 x 2ϕ xϕ 2 ϕ 3 x 3ϕ x 2ϕ 2 xϕ 3 x 3ϕ 2 x 2ϕ 3 x 3ϕ 3 
R 
.... ,
R R R R R R R R R R R
0 0 0 0 0 0 0 0 0 0 0 0 
(17.20)

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

6 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems

0 0 0 0 0 0 0 0 0 0 0 2 0 0 6x 2ϕ
 1 x 2
R 1 = 0 0 0 0 0 0 0 0 0 0 0 − 2 0 0 ....
 R2 R2 R
0 sin ϕ −
cos ϕ
0 0 0 0 0 0 0
ϕ
0 −
1
0 0 −
2x
 2R 2R 2R R R

0 0 6 xϕ 2ϕ 2 0 2ϕ 3 6 xϕ 2
6 xϕ 3 
2x 6ϕ 2x2 6 xϕ
6 x 2ϕ 6 x 3ϕ 
2x3
.... − 2 − 2 0 − 2 − 2
− 2 − 2 .
− 2
R R R R R R 
R
2ϕ 3x 2 4 xϕ 3ϕ 2
6 xϕ 2 6 x 2ϕ
9 x 2ϕ 2 
− 0 − − −− − −
R R R R R R 
R
(17.21)
For the calculation of the stiffness matrix and the force vector related to the distributed
load we formulate the total potential energy of a single element. We note that the expres-
sion below contains only the strain energy and the work of the distributed load:
1
Π e = ∫ σ ε dV − ∫ u pdA .
T T
(17.22)
2 Ve A pe

The constitutive law of the linear elastic material is:

σ = Cε , (17.23)
where, similarly to the plates we assume plane stress state, i.e. C = C . Using Eqs.(17.17)
str

and (17.19) we obtain:

ε 0 = R 0 A = R 0 M −1 u~ e = B 0 u~ e , (17.24)
−1
B0 = R0 M , σ 0 = C ε 0 = C B 0 u~ e ,
furthermore:
κ = R 1 A = R 1 M −1 u e = B 1 u~ e , (17.25)
−1
B 1 = R 1 M , σ 1 = − z Cκ = − z C B 1 u~ e
Based on Eq.(17.7) the displacement field becomes:
−1
u = λ A = λ M u~ e , (17.26)
where:
λ = [Φ 0 Φ1 ]. (17.27)
Eq.(17.24) is related to the stress resultants, while Eq.(17.25) is related to the stress cou-
ples. The total potential energy becomes:
1 T T 1
Π e = ∫ ε C ε dV + ∫ z 2 κ C κ dV − u~ e ∫ M λ pdA , (17.28)
T T T −1T T

2 Ve 2 Ve A pe

which is written as:

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 7

L/2 θ
1 ~T  ~ 1 ~T  L / 2 θ t 3 T T 
Π e = u e  ∫ ∫ t ⋅ B 0 C B 0 R ⋅ dϕ ⋅ dx u e + u e  ∫ ∫ B1 C B1R ⋅ dϕ ⋅ dx u~ e +
T T

2 − L / 2 −θ  2 − L / 2 −θ 12 
L/2 θ

∫ ∫θ M
−1T
− u~ e λ T p ⋅ R ⋅ dϕ ⋅ dx.
T

−L / 2 −
(17.29)
It is important to note that in Eq.(17.29) the term related to the concentrated loads is ex-
cluded, consequently the vector of concentrated loads should be produced additionally.
This is an easy task based on the nodal degrees of freedom:
~T
[
F ec = Fx1 Fϕ1 Fz1 M x1 M ϕ1 M xϕ1 Fx 2 Fϕ 2 Fz 2 M x 2 M ϕ 2 M , xϕ 2 .....
.........Fx 3 Fϕ 3 Fz 3 M x3 Mϕ3 M xϕ 3 Fx 4 Fϕ 4 Fz 4 M x4 Mϕ 4 ]
M xϕ 4 ,
(17.30)
and the completed total potential energy becomes:
1 T~ T ~ ~
Π e = u~ e K e u~ e − u~ e ( F ep + F ec ) . (17.31)
2
~
In Eq.(17.31) K e is the element stiffness matrix in the local coordinate system
L/2 θ
~  t3 T T 
∫ ∫ 0
Ke =  ⋅ + B C B 1  R ⋅ dϕ ⋅ dx .
T T
t B C B 0
(17.32)
− L / 2 −θ 
12 1 
The force vector from the distributed load is:
L/2 θ
~
F ep = ∫ ∫ M λ p ⋅ R ⋅ dϕ ⋅ dx .
−1T T
(17.33)
− L / 2 −θ

Finally the well-known finite element equilibrium equation in the local system is:
~ ~
K e u~ e = F e , (17.34)
which is applicable only for a single element. The global equation of developed by a prop-
er transformation. The structural equation is required when there are several elements con-
nected to each other, which is mathematically the same as Eq. (14.86). The advantage of
the thin cylindrical shell element is that the cylindrical surface is captured exactly; as a
consequence it provides accurate result even if the number of elements is relatively low.

17.3 Axisymmetric shell problems – conical shell element

The midsurface of axisymmetric shells is produced by the rotation of the meridian curve
about a straight axis [1]. An example is shown by Fig.17.2.

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

8 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems

Fig.17.2. Axisymmetric shell.

The meridian curves and the circular curves perpendicularly to the meridian curves are
principal curvature lines of the surface. If the load of the structure is axisymmetric, then in
this kind of problem the displacement field is the function of arc length along the meridian
curve only.
The meridian curve can be modeled by straight lines, and so we approximate the original
shell structure by conical shell elements. Referring to the basic equations of the technical
theory of thin shells, the parameters of the conical shell element shown in Fig.17.3 are:
q1 = s , H 1 = 1 , R1 = ∞ , (17.35)
r
q 2 = ϕ , H 2 = r , R2 = ,
cos θ
where s is the arc length, ϕ is the angle coordinate, r is the radius for a point P, θ is the half
angle of inclination. To calculate the strain components we need to determine the r(s) rela-
tionship, based on Fig.17.3 we have:
dr
r ( s ) = s sin θ + r1 and = sin θ . (17.36)
ds

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 9

Fig.17.3. Axisymmetric conical shell element and its nodal parameters.

Using Eqs.(14.67), (14.69) and (14.70) of the technical theory of thin shells we can calcu-
late the strain components as:
β1 = β s = − w, s , β 2 = 0 , β 3 = 0 , (17.37)
1
ε 11 = ε s = u , s , ε 22 = ε ϕ = (u sin θ + w cos θ ) , γ 12 = 0 ,
r
sin θ
κ 11 = κ s = − w, ss , κ 22 = κ ϕ = − w, s , κ 12 = 0 .
r
The displacement in the tangential direction at point P is v = 0 due to the axisymmetry. The
admissible rigid body-like motion of the element is a displacement given by d in direction
Z, for which the displacement components are u = -d⋅cosθ and w = d⋅sinθ (see Fig.17.3).
We consider three degrees of freedom at each node, these are: u (displacement along the
meridian direction), w (displacement perpendicularly to the meridian curve) and βs (angle
of rotation about the axis perpendicularly to the meridian curve in accordance with
Fig.17.3), therefore the element has six degrees of freedom. The displacement in the me-
ridian direction is interpolated by a linear function of the arc length. On the other hand we
apply third order interpolation with respect to the displacement in the normal direction:
 u  1 s 0 0 0 0 
 w = 0 0 1 s s 2 s 3 α = Φα , (17.38)
   
where α is the vector of unknown coefficients:
α T = [a1 a 2 a3 a 4 a5 a 6 ] . (17.39)
The vector of nodal displacements is:
[
u~ e = u~1 w
T ~ β~ u~ w
1 s1 2 2 ]
~ β~ .
s2 (17.40)
The conditions required for the determination of the coefficients are:
u ( s1 ) = a1 + a 2 s1 = u~1 , (17.41)
w( s1 ) = a3 + a 4 s1 + a5 s1 + a 6 s1 = w1 ,
2 3 ~
~
β s ( s1 ) = a 4 + 2a5 s1 + 3a 6 s12 = β s1 ,
u ( s 2 ) = a1 + a 2 s 2 = u~2 ,
w( s 2 ) = a3 + a 4 s 2 + a5 s 22 + a 6 s 23 = w ~ ,
2

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

10 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems
~
β s ( s 2 ) = a 4 + 2a5 s 2 + 3a 6 s 22 = β s 2 .
The solutions for the coefficients are moderately complicated, therefore they are not in-
cluded here. The displacement functions can be formulated also in the way presented be-
low:
u ( s ) = N 1u~1 + N 2 u~2 , (17.42)
~ + N β~ + N w ~ ~
w( s ) = N 3 w 1 4 s1 5 2 + N 6 β s2 ,

where Ni, i = 1...6 are the interpolation functions:

s −s s −s
N1 = 2 , N2 = − 1 , (17.43)
s 2 − s1 s 2 − s1
( s − s 2 ) 2 (3s1 − s 2 − 2 s ) ( s − s 2 ) 2 ( s1 − s )
N3 = , N = ,
( s 2 − s1 ) 3 ( s 2 − s1 ) 2
4

( s − s1 ) 2 (3s 2 − s1 − 2 s ) ( s − s1 ) 2 ( s 2 − s )
N5 = , N = − ,
( s 2 − s1 ) 3 ( s 2 − s1 ) 2
6

and:
u  N 0 0 N2 0 0 ~
u= = 1  u e = N u~ e . (17.44)
 w  0 N 3 N 4 0 N 5 N 6 
The strain components in matrix form are:
ε s 
ε =   = Bu e , (17.45)
ε
  ϕ

where the strain-displacement matrix is:

 ∂N 1 ∂N 2 
 ∂s 0 0 0 0 
B= ∂s
N 1 sin θ N 3 cos θ N 4 cos θ N 2 sin θ N 5 cos θ N 6 cos θ 
.
 
 r r r r r r 
(17.46)
We collect also the curvatures in matrix form:
κ s 
κ =   = H u~ e , (17.47)
κ ϕ 
where:
 ∂ 2 N3 ∂2 N4 ∂ 2 N5 ∂2 N6 
0 − − 0 − − 
H = ∂s 2 ∂s 2 ∂s 2 ∂s 2  .
0 − sin θ ∂N 3 − sin θ ∂N 4 0 − sin θ ∂N 5 − sin θ ∂N 6 
 r ∂s r ∂s r ∂s r ∂s 
(17.48)
Based on the constitutive law the vector of stress components is:
σ 0 = C ε = C Bu~ e , (17.49)
σ = − z Cκ = − z C H u~ .
1 e

The vectors of strain components and curvatures contain only two elements, therefore the
constitutive matrix reduces to:
E 1 ν 
C=  . (17.50)
1 − ν 2 ν 1 

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 11

Taking the former back into the total potential energy (similarly to the cylindrical shell
element) we can calculate the element stiffness matrix in the local coordinate system:
s2 l
~ t3 T T t3 T T
K e = ∫ (t ⋅ B C B + H C H ) ⋅ 2πrds = ∫ (t ⋅ B C B + H C H ) ⋅ 2π ( s sin θ + r1 )ds .
T T T T

s1
12 0
12
(17.51)
In the above expression it was considered that s1 = 0 and s2 = l and so r = s⋅sinθ. The exact
computation of the stiffness matrix is quite complicated, and consequently the finite ele-
ment codes implement numerical methods, e.g. the Gauss rule presented in section 12 is
suitable to calculate the matrix components. The force vector is composed by two terms.
The vector of concentrated forces can be constructed based on the nodal degrees of free-
dom:
F ec = [Fs1 Fn1 M 1 Fs 2 Fn 2 M 2 ] ,
~T
(17.52)
where F refers to the concentrated force, M is a concentrated moment about the same di-
rection tan that of βs. The force vector from the distributed load is calculated based on the
work of the load:
T p 
l l
T ~
We = ∫ ( p s u + p n w)2πrds = u~ e ∫ N  s  ⋅ 2πrds = u~ e F ep , (17.53)
T

0 0  pn 
accordingly:
T p 
l
~
F ep = ∫ N  s  ⋅ 2π ( s sin θ + r1 )ds . (17.54)
0  p n 
Considering that l⋅sinθ = r2-r1 and assuming that both ps and pn are constants, we obtain:
 1 
 3 p s πl (2r1 + r2 ) 
 1 
 p nπl (7 r1 + 3r2 ) 
 10 
 1 p πl 3 (3r + 2r ) 
~  n 1 2 
F ep =  30 . (17.55)
1
 p s πl (r1 + 2r2 ) 
 3 
 1 p πl (3r + 7 r ) 
 10 n 1 2

 1 
− p nπl 3 (2r1 + 3r2 )
 30 
In the local coordinate system the nodal displacement and reactions are calculated form the
usual:
~ ~
K e u~ e = F e (17.56)
equation, where:
~ ~ ~
F e = F ec + F ep . (17.57)
For a finite element structure we need the structural equation given by Eq.(14.86). Since
the elements are connected under a given angle, the local displacement coordinates should
be transformed into the global cylindrical coordinate system with longitudinal axis given
by Z. The transformation can be performed based on Fig.17.3:

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

12 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems

 u  cos θ − sin θ 0  u~ 
 w  =  sin θ cos θ 0  w ~ T ~
    ~  = L u . (17.58)
 β s   0 0 1  β s 
Based on the former the transformation of the stiffness matrix becomes:
T ~
Ke = λ Keλ , (17.59)
where:
LT 0 
T = T , (17.60)
 0 L 
is an orthogonal transformation matrix. The transformed force vector is:
~
Fe =TFe. (17.61)
For a single element the finite element equation in the global system is:
K eue = F e , (17.62)
Moreover, for the whole structure we have:
KU = F . (17.63)
In the finite element literature there are more element types, e.g. curved axisymmetric shell
element [4,5,6], which operates similarly to the conical shell element.

17.4 Thick-walled shell elements

For the solution of three-dimensional problems we can apply the spatial (SOLID type) el-
ements. Fig.17.4 shows a 20 node isoparametric element. Isoparametric representation
means that the geometry and the displacement field is described by the same set of interpo-
lation functions [1,4,5]:
 x  20  xi   u  20  ui 
 y =  
N i (ξ ,η , ζ )  y i  ,  v  = ∑ N i (ξ ,η , ζ )  vi  . (17.64)
 
  ∑
 z  i =1  z i   w i =1  wi 

Fig.17.4. Quadratic two and three dimensional elements.

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 13

The thick-walled shell elements are constructed in accordance with isoparametric formula-
tion, in this respect we point out that the sides perpendicularly to the shell midsurface are
straight, i.e. the interpolation in the thickness direction is linear. The element is determined
by the 8 nodes of the ζ = 0 midsurface. As it can be seen in Fig.17.4 the direction of the
unit basis vectors changes from point to point, therefore the nodal number is indicated by
subscript „i”. The coordinates of the points on the midsurface of the thick-walled shell
element are given by:
 x  xi 
 y = ζ
N i (ξ ,η )  yi  +
8 8

  ∑ ∑ N (ξ ,η ) 2 t n ,
i i i (17.65)
 z  i =1
 zi  i =1

where ni are the column vector of normal vectors at the midsurface nodes, ti is the thick-
ness in the actual node, Ni are the interpolation functions, respectively. The interpolation
functions are the same as those of the quadratic isoparametric plane membrane element
(see section 12). The compact form of the interpolation function is:
1
N i = (1 + ξξ i )(1 + ηη i )(ξξ i + ηη i − 1) , i = 1, 3, 5, 7, (17.66)
4
1 2 1
N i = ξ i (1 + ξξ i )(1 − η 2 ) + η i2 (1 + ηη i )(1 − ξ 2 ) , i = 2, 4, 6, 8,
2 2
where ξi and ηi are the local nodal coordinates. On the midsurface the ξ and η coordinate
lines are orthogonal, therefore the basis vectors are calculated as:
R × R i ,η R
n i = i ,ξ , e 2i = i ,η , e1i = e 2i × ni . (17.67)
R i ,ξ × R i ,η R i ,η
The nodal displacement parameters are the ui, vi, wi displacements and the β1i and β2i angle
of rotations. In the case of eight nodes it means that the element has 40 degrees of freedom.
Vector ni can be formulated by using the rotations and the basis vectors e1i, e2i:
~ ~
ni = β 2 e1i − β1 e 2i , (17.68)
which is the term capturing the transverse shear deformation, it causes an increment in u
and v. According to the isoparametric representation the displacement field becomes:
u  8  u~i  8
 v  = N (ξ ,η )  v~  + N (ξ ,η ) ζ t ( β~ e − β~ e ) . (17.69)
  ∑ i =1
i  i ∑ i =1
i
2
i 2 i 1i 1i 2 i

 w ~
 wi 
To calculate the stiffness matrix we have to establish the strain-displacement relationship.
The derivatives of the displacement parameters with respect to the local coordinates are:
 ∂u   ∂N i 1 1 x 
 ∂ξ [1 − 2 tiζe2i 2 tiζe1i ] ~
x
 ∂ξ 
  8    ui 
∂ ∂
tiζe1i ]  β1i  .
u
  = ∑ N 1 1 ~
i
[1 − tiζe2i
x x
(17.70)
 ∂η  i =1  ∂η 2 2  ~
 ∂u   1 x 1 x   β 2i 
   N i [0 − ti e2i ti e1i ] 
 ∂ζ   2 2 
For the other two components we obtain similar equations. The further computations re-
quire the Jacobi matrix and determinant [1,4,5]:

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

14 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems
 ∂  ∂ ∂  ∂ 
 ∂ξ   ∂x   ∂x   ∂ξ 
       
 ∂  = J  ∂  , and:  ∂  = J −1  ∂  . (17.71)
 ∂η   ∂y   ∂y   ∂η 
 ∂  ∂ ∂  ∂ 
       
 ∂ζ   ∂z   ∂z   ∂ζ 
The elements of the Jacobi matrix can be obtained using Eq.(17.65). Also, the derivatives
of the displacement components can be determined in the global coordinate system. For
example, the derivatives of the component u in matrix form are:
 ∂u   ∂u 
 ∂x   ( −1)  
 ∂u   J 11 J 12( −1) J 13( −1)   ∂ξ 
( −1)   ∂u 
  =  J 21( −1) ( −1)
J 22 J 23   , (17.72)
 ∂y   ( −1) ( −1) ( −1)  ∂η
 ∂u   J 31 J 32 J 33   ∂u 
   
 ∂z   ∂ζ 
(-1)
where Jij are the elements of the inverse Jacobi matrix. Based on Eq.(17.70) we obtain
the following:
∂u 8
∂N ∂N 1  ∂N ∂N ~
= ∑ ( J11( −1) i + J12( −1) i )u~i − ti e2xi ζ ( J11( −1) i + J12( −1) i ) + J13( −1) N i β1i +
∂x i =1 ∂ξ ∂η 2  ∂ξ ∂η 
1  ∂N ∂N ~ 8
∂N ~ ~
+ ti e1xi ζ ( J11( −1) i + J12( −1) i ) + J13( −1) N i β 2i = ∑ i u~i + Gix ( g1xi β1i + g 2xi β 2i ),
2  ∂ξ ∂η  i =1 ∂x

(17.73)
where:
∂N ∂N
Gix = ζ ( J11( −1) i + J12( −1) i ) + J13( −1) N i , (17.74)
∂ξ ∂η
and:
1 1
g 1 = − t i e 2i , g 2 = t i e1i .
i i
(17.75)
2 2
The derivatives with respect to the other two coordinates are:
∂u 8
( −1) ∂N i ( −1) ∂N i ~ 1  ( −1) ∂N i ( −1) ∂N i ~
= ∑ ( J 21 + J 22 )ui − ti e2xi ζ ( J 21 + J 22 ( −1)
) + J 23 N i β1i +
∂y i =1 ∂ξ ∂η 2  ∂ξ ∂η 
1  ( −1) ∂N i ( −1) ∂N i ~ 8
∂N ~ ~
+ ti e1xi ζ ( J 21 + J 22 ( −1)
) + J 23 N i β 2i = ∑ i u~i + Giy ( g1xi β1i + g 2xi β 2i ),
2  ∂ξ ∂η  i =1 ∂y

(17.76)
∂u 8
∂N ∂N 1  ∂N ∂N ~
= ∑ ( J 31( −1) i + J 32( −1) i )u~i − ti e2xi ζ ( J 31( −1) i + J 32( −1) i ) + J 33( −1) N i β1i +
∂z i =1 ∂ξ ∂η 2  ∂ξ ∂η 
1  ∂N ∂N ~ 8
∂N ~ ~
+ ti e1xi ζ ( J 31( −1) i + J 32( −1) i ) + J 33( −1) N i β 2i = ∑ i u~i + Giz ( g1xi β1i + g 2xi β 2i ),
2  ∂ξ ∂η  i =1 ∂z

where:
( −1) ∂N i ( −1) ∂N i
Giy = ζ ( J 21 + J 22 ( −1)
) + J 23 Ni , (17.77)
∂ξ ∂η

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 15

∂N i ∂N
Giz = ζ ( J 31( −1) + J 32( −1) i ) + J 33( −1) N i .
∂ξ ∂η
Written in matrix form we have:
 ∂u   ∂N i x x 
 ∂u 
 ∂x Gi g1i Gi g 2i   ∂ξ 
x x
 ∂x 
 ∂u  8  ∂N   ∂u 
  = ∑  i Giy g1xi Giy g 2xi    . (17.78)
 ∂y  i =1  ∂y   ∂η 
 ∂u   ∂N i  
z x  ∂u
i 2i   
z x
   G g G g
 ∂z   ∂z   ∂ζ 
i 1i

The derivatives of the other two components can be provided similarly. Using the deriva-
tives we can calculate matrix B , which is the relationship between the strain components
and the nodal displacement parameters:
ε~ = Bu~ e , (17.79)
where u~ e is the vector of nodal parameters in the local coordinate system. The vectors of
strain and stress components in the global system are:
[
ε T = ε x ε y ε z γ xy γ yz γ xz , ] (17.80)
σT = [σ x σ y σ z τ xy τ yz τ xz].
Hooke’s law in the local system can be written as:
σ~ = C ε~ , (17.81)
where C is the constitutive matrix:
1 ν 0 0 0 0 
ν 1 0
 0 0 0 
0 0 0 0 0 0 

E 0 0 0 1 − ν 
C=  0 0 . (17.82)
1 −ν 2  2 
1 −ν
0 0 0 0 k 0 
 2 
 1 −ν 
0 0 0 0 0 k
2 
The matrix above differs from the general three dimensional case in accordance with the
followings. The stress normal to the shell surface is zero (3rd row, 3rd column). Since the
element is thick-walled it considers also the effect of transverse shear deformation, but
only in the form of an average stress. The constant in the elements of the 5th row, 5th col-
umn, and the 6th row, 6th column is a shear correction factor, k = 5/6 [1,4,5]. The reason for
that is the real distribution of the shear stresses is assumed to be parabolic over the thick-
ness, and it is not constant as considered in the shell model. The correction factor k is the
ratio of the strain energies from the two different distributions. Based on the transformation
of local stress and strain components we can write the followings:
σ~ = T σ , σ = T T σ~ , (17.83)
ε~ = T ε , ε = T T ε~ ,

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

16 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems
where T is the transformation matrix for general spatial stress and strain states. The calcu-
lation of T is possible using the definitions given by Eq.(11.62). Taking back Eq.(17.83)
into Hooke’s law we have:
T σ = CT ε . (17.84)
−1
The premultiplication with T leads to:
−1 −1
T T σ = T CT ε . (17.85)
−1 −1
SinceT is an orthogonal matrix we can write that: T T = E , and T = T , viz.:
T

σ = T T CT ε . (17.86)
The transformation matrix is [4]:
 l12i m12i n12i l1i m1i m1i n1i n1il1i 
 2 
 l2 i m22i n22i l2i m2i m2i n2i n2il2i 
 l 2 2
m3i 2
n3i l3i m3i m3i n3i n3il3i 
T =  3i ,
 2l1il2i 2m1i m2i 2n1i n2i l1i m2i + l2i m1i m1i n2i + m2i n1i n1il2i + n2il1i 
2l l 2m2i m3i 2n2i n3i l2i m3i + l3i m2i m2i n3i + m3i n2i n2il3i + n3il2i 
 2 i 3i 
 2l3il1i 2m3i m1i 2n3i n1i l3i m1i + l1i m3i m3i n1i + m1i n3i n3il1i + n1il3i 
(17.87)
where li, mi and ni are the direction cosines of the unit basis vectors at the actual point
[4,7]:
l1i = cos(i, e1i ) , m1i = cos( j , e1i ) , n1i = cos(k , e1i ) , (17.88)
l2i = cos(i, e 2i ) , m2i = cos( j , e 2i ) , n2i = cos(k , e 2i ) ,
l3i = cos(i, e3i ) , m3i = cos( j , e3i ) , n3i = cos(k , e3i ) .
The transformation matrix should be evaluated in the nodes, moreover due to the numerical
integration even in the integration points. The stiffness matrix in the global coordinate sys-
tem can be calculated using Eq.(15.17):
K e = ∫ B T C T BdV = ∫ B T C T BJdξdηdζ ,
T T T T T T
(17.89)
Ve Ve

where J is the Jacobi determinant, which can be calculated using Eqs.(17.65) and (17.71).
For the determination of the force vector we recall the displacement vector field in the usu-
al form:
u (ξ ,η , ζ ) = N u~ e , (17.90)
where N is the matrix of interpolation polynomials. As a result, the vectors of body, surface
and line forces in the global coordinate system are:
~
F eb = ∫ N p b dV = ∫ N p b Jdξdηdζ ,
T T
(17.91)
Ve Ve
~
F ep = ∫ N p p dA = ∫ N p p Jdξdη ,
T T

Ae Ae
~
F el = ∫ N p l dS ,
T

which can be determined by transformation into the global system in a similar way to that
presented in section 16. In the nodes concentrated forces may act, the relevant vector can

 www.tankonyvtar.hu  Dr. András Szekrényes, BME

Alfejezetcím 17

be obtained in the same way as that shown in plate elements. Because of he high number of
nodes it is not detailed here. The finite element equilibrium equation is formed in the usual
way, for a single element it is:
K eue = F e , (17.92)
where Fe is the sum of the vectors of body, surface, line and concentrated forces. Finally,
the structural equation is:
KU = F . (17.93)

17.5 A shell-solid transition element

In complex structures sometimes there is the necessity of the simultaneous application of
solid and thick-walled shell elements. These elements can not be connected directly, be-
cause the nodal degrees of freedom are not identical. In these cases it is reasonable to use a
transition element between the solid and shell elements [1,4,5]. A quadratic transition ele-
ment is shown in Fig.17.5, where the nodes 1-8 are located in the solid side, nodes 10-12
care located in the shell side of the element.

Fig.17.5. A shell-solid transition element.

The geometry of the transition element is captured by the function below:

 x 8  xi  13  xi  13
 y  = N (ξ ,η , ζ )  y  + N (ξ ,η )  y  + N (ξ ,η ) ζ t n .
  ∑ i =1
i  i ∑ i =9
i  i ∑ i =9
i
2
i i

 z   zi   zi 
(17.94)
The indices i = 1...8 refer to the interpolation function of the solid element given by
Eq.(17.64), if i = 9...13 then the actual interpolation functions of the thick-walled shell el-
ements are referred to in accordance with Eq.(17.65). The composed system of functions
satisfies the following conditions [1]:
8 13

∑ N (ξ ,η , ζ ) + ∑ N (ξ ,η ) = 1 ,
i =1
i
i =9
i (17.95)

0, if i ≠ j
N i (ξ j ,η j , ζ j ) =  ,
1, if i = j

 Dr. András Szekrényes, BME  www.tankonyvtar.hu

18 Modeling of curved and doubly-curved shells by finite element method based soft-
ware systems
0, if i ≠ j
N i (ξ j ,η j ) =  ,
1, if i = j
where ξj, ηj and ζj are the nodal coordinates in the local coordinate system. Similarly to the
thick-walled shell elements the displacement field is expressed by:
u  8  u~i  13  u~i  13
 v  = N (ξ ,η , ζ )  v~  + N (ξ ,η )  v~  + N (ξ ,η ) ζ t ( β~ e − β~ e ) .
  ∑ i =1
i  i ∑ i =8
i  i ∑ i =8
i
2
i 2 i 1i 1i 2 i

 w ~
 wi  ~
 wi 
(17.96)
The degrees of freedom in nodes 1-8 are equal to three, in nodes 9-13 there are five de-
grees of freedom. Consequently the transition element has 49 degrees of freedom. The fur-
ther calculations can be performed in similar fashion to that presented in the thick-walled
shell element.

17.6 Bibliography
[1] Imre Bojtár, Gábor Vörös, Application of the finite element method to plate and
shell structures. Technical Book Publisher, 1986, Budapest (in Hungarian).
[2] Gábor Vörös, Lectures and practices of the subject Applied mechanics, manuscript,
Budapest University of Technology and Economics, Faculty of Mechanical Engi-
neering, Department of Applied Mechanics, 1978, I. semester, Budapest (in Hun-
garian).
[3] Singiresu S. Rao, The finite element method in engineering – fourth edition, Else-
vier Science & Technology Books, 2004,
[4] Klaus-Jürgen Bathe, Finite element procedures. Prentice Hall, Upper Saddle River,
1996, New Jersey 17458.
[5] Thomas J. Hughes, The finite element method – Linear static and dynamic finite
element analysis. Prentice Hall Inc., A division of Simon & Schuster, Englewood
Cliffs, 1987, New Jersey 07632.
[6] O.C. Zienkiewicz, R.L. Taylor. The finite element method – fifth edition, Volume 2:
Solid mechanics. Butterworth-Heinemann, 2000, Oxford, Auckland, Boston, Jo-
hannesburg, Melbourne, New Delhi.
[7] Pei Chi Chou, Nicholas J. Pagano, Elasticity – Tensor, dyadic and engineering ap-
proaches. D. Van Nostrand Company, Inc., 1967, Princeton, New Jersey, Toronto,
London.

 www.tankonyvtar.hu  Dr. András Szekrényes, BME