# Approximation of Linear Elastic Shells by Curved Triangular Finite Elements Based on Elastic Thick Shells Theory

Authors: Joseph Nkongho Anyi; Robert Nzengwa; Jean Chills Amba; Claude Valery Abbe Ngayihi
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We have developed a curved finite element for a cylindrical thick shell based on the thick shell equations established in 1999 by Nzengwa and Tagne (N-T). The displacement field of the shell is interpolated from nodal displacements only and strains assumption. Numerical results on a cylindrical thin shell are compared with those of other well-known benchmarks with satisfaction. Convergence is rapidly obtained with very few elements. A scaling was processed on the cylindrical thin shell by increasing the ratioχ
=
h
/
2
R (half the thickness over the smallest radius in absolute value) and comparing results with those obtained with the classical Kirchhoff-Love thin shell theory; it appears that results diverge at 2
χ
=
1
/
10
=
0.316 because of the significant energy contribution of the change of the third fundamental form found in N-T model. This limit value of the thickness ratio which characterizes the limit between thin and thick cylindrical shells differs from the ratio 0.4 proposed by Leissa and 0.5 proposed by Narita and Leissa.


---

## Body

## 1. Introduction

In the field of structural mechanics the word shell refers to a spatial, curved structural member. The very high structural and architectural potential of shell structures is used in various fields of civil, architectural, mechanical, aeronautical, and marine engineering. The strength of the double-curved structure is efficiently and economically used, for example, to cover large areas without supporting columns [1]. In addition to the mechanical advantages, the use of shell structures leads to aesthetic architectural appearance [2, 3]. Examples of shells used in civil and architectural engineering are shell roofs, liquid storage tanks, silos, cooling towers, containment shells of nuclear power plants, and arch dams. Piping systems, curved panels, pressure vessels, bottles, buckets, and parts of cars are examples of shells used in mechanical engineering [4]. In aeronautical and marine engineering, shells are used in aircrafts, space crafts, missiles, ships, and submarines. Because of the spatial shape of the structure the behavior of shell structures is different from the behavior of beam and plate structures [4].The considerable effort in the development of rigorous shell theories dates back to the early twentieth century [5–7]. These shell theories reduce a basically three-dimensional problem to a two-dimensional one. Nevertheless, the analysis of shells with the aid of such theories involves complicated differential equations, which either cannot be solved at all [8–10], or whose solution requires the use of high-level mathematics unfamiliar to structural engineers. Therefore, many approximate shell theories have been developed, mainly on the assumption that the shell is thin, and to obtain generic analysis tools obviously some accuracy had to be traded for convenience and simplicity [4, 11]. Hence, it is not surprising that the development of the numerical formulations since the 1950s has led to a gradual cessation of attempts to find closed-form solutions to rigorous formulations [4]. But, with today’s availability of greatly increased computing power (also since the mid-twentieth century), completeness rather than simplicity is given more emphasis.In computational mechanics, many numerical techniques are developed to solve governing equations of shells. We can identify for this purpose the generalized differential quadrature finite elements (DGQFEM) [11–13], discrete singular convolution method (DSC) [14], the finite difference (FD), the discrete volume finite element (DVFEM), and the finite elements method (FEM) under which solid-shells finite elements (SSFEM) and mixt finite element (MFEM) are found [15–18]. Our approach is found under the MFEM category in which three-dimensional shell parameters are reduced to the midsurface of the shell of thickness h. Some finite elements in this category describe well the behavior of shells under several loads. All these models of finite elements are link with either Kirchhoff-Love (K-L) or Reissner-Mindlin (R-M) hypothesis which neglected the effect of curvature in elastic stiffness [11, 19]. This is one of the main reasons why K-L and R-M theories of shells are incomplete and unsuitable to handle various types of shells [20]. Structure engineers emphasis more on simplifying assumptions than limit analysis. Because of the third fundamental form that has been neglected in their theoretical approach of shells governing equations [21] and severe assumption on shell thickness, some mathematical terms without evident mechanic signification are introduced in the stiffness matrix in order to improve their efficiency [6, 22].After a brief presentation of Nzengwa and Tagne (N-T) kinematic equations of elastic thick shells, we develop a curved cylindrical 3-node finite element based on both K-L and N-T’s shell models using strain tensor interpolation assumption [23]. Next, through numerical implementation on classical benchmark, convergence of N-T’s shell model is investigated. The deviation and the limit thickness ratio between the two theoretical approaches of thin shells and thick shells are established and discussed by applying successive scaling of a sample of thin shell.


## 2. Materials and Methods

### 2.1. N-T’s Two-Dimensional Elastic Thick Shells Model

The N-T’s two-dimensional model for linear elastic thick shells has been deduced from the three-dimensional problem without any ad hoc assumption whether of geometrical or mechanical nature. The two-dimensional equations are deduced by applying asymptotic analysis on a family of variational equations obtained from an abstract scaled shell through multiple scaling of the initial three-dimensional equations. The theoretical elastic thick shells governing equations and limits displacements obtained from this approach (at the limit analysis) are more general because they contain additional terms to those found in Kirchhoff-Love model. The model also completes the thick shells theory of Reissner-Mindlin [24].For a three-dimensional shell, the following basesg
1
,
g
2
,
g
3
,
g
1
,
g
2
,
g
3 are dual bases and a
1
,
a
2
,
a
3
,
a
1
,
a
2
,
a
3 defined on the midsurface also constitute dual bases of shell if χ
<
1. Let U
i
x
1
,
x
2
,
z be the three-dimensional displacement field vector; it is defined as follows in the local coordinate system: (1)U=Uαgα+U3g3=u-αaα+u-3a3,u-ρ=u-ρx1,x2,z,u-3=u-3x1,x2,z.The strain tensor for a three-dimensional shell reads(2)εαβU=12Uα/β+Uβ/α=εαβu-=μαν∇βu-ν-bαβu-3+μβν∇αu-ν-bνβu-32,(3)εα3U=12Uα/3+U3/α=εα3u-=μανu-ν,3+u-3,α+bανu-ν2,(4)ε33U=U3/3=U3,3=u-3,3=ε33u-.From the limit analysis N-T demonstrated that the displacements field satisfies the equation ε
i
3
(
U
)
=
0 and the unique solution is(5)u-ρx1,x2,z=μργuγx1,x2-z∂ρw,W=w=u3,Uαx1,x2,z=μαρu-ρ,μργ=δργ-zbργ,Uαx1,x2,z=uαx1,x2-zθαx1,x2+z2ψαx1,x2,where (
u
α
) are local membrane displacement components and w
=
u
3 the transvers displacement:(6)θα=2bατuτ+∇αw,ψα=bαγbγτuτ+bαγ∇γw.The general shell deformation tensor derived from the kinematics defined above is(7)εαβU=eαβu-zkαβu+z2Qαβu,(8)εi3U=0,where e
α
β is the membrane strain tensor, k
α
β is the bending strain tensor, and Q
α
β is the tensor of Gauss curvature.The strain tensor contains the change of the third fundamental form:(9)Qαβ=12cαβ-Cαβ,c=cαβ=bαμbμβ=bαμbβμ,Cαβ=BαμBμβ=BαμBβμ.In Kirchhoff-Love thin shells theory, the shell is assumed to be thin. Following that assumption, the three-dimensional displacements field vectorU
α
x
1
,
x
2
,
z is taken equal to the midsurface displacements field:(10)Uαx1,x2,z=uαx1,x2-z∇αwx1,x2,Wx1,x2,z=wx1,x2=u3.It should be noted that the terms ψ
α
=
b
α
τ
b
τ
ρ
u
ρ
+
b
α
τ
∇
τ
w in the current model of thick shells disappear in the classical theory of thin shells because it is proportional to χ
2 and considered small. The thin shell theory, with the above assumption, cannot perform good shells behavior because the thickness has a considerable influence on the shells behavior. The thick shell model used in this survey appears to be suitable to handle both thin and thick shells because of its completeness.


### 2.2. Curved Shell Elements

Curved shell elements are suitable to model the midsurface geometry more accurately. In the case of certain surfaces such as cylindrical shells, it means the exact description of the original surface. For more complicated cases, similar 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 [25].


### 2.3. Triangular Curve Shell Element

The triangular shell element is described in Figure1; it can be double-curved (e.g., spherical structures) or single-curved element (cylindrical structures). In accordance with the N-T’s equations of the theory of elastic thick shells, the geometrical properties of the cylindrical shell are the following:(11)OM→=X=Rcos⁡φ,Y=Rsin⁡φ,Z=x,where x and φ are curvilinear coordinates, 0
≤
x
≤
L and 0
≤
φ
≤
θ, R
≡ cylinder radius, and X
,
Y
,
Z are global coordinates.Figure 1
Triangular shell element.


### 2.4. Displacement Field

LetU
α be the global displacement field vector defined as follows: (12)Uα=uα-zθα+z2ψα,W=w.With respect to the cylinder’s geometry, the angles of rotations and the Gaussian curvatures are developed in function of membrane displacements u
1 and u
2 and the derivatives of the transverse displacement w
,
φ
,
w
,
x including curvature components b
α
τ:(13)θ1=2u1R1-w,xA1,θ2=2u2R2-w,φA2,ψα=bατbτρuρ+bατ∇τw,which leads to(14)ψ1=-1R1u1R1-1A1w,x=1R1b1τuτ+θ1,ψ2=-1R2u2R2-1A2w,φ=1R2b2τuτ+θ2,where A
1
,
A
2
,
R
1, and R
2 are geometrical properties of the midsurface of a shell within the respective coordinate-lines x and φ.In the case of a cylindrical shell they are defined as follows [25]:(15)A1=a1→=∂M∂x→=1;R1=∞,A2=a2→=∂M∂φ→=R;R2=R.The third component of rotation is added for a smooth transfer of element stiffness matrix from one system of coordinates to another.


### 2.5. Compatibility

We recall the strain components in thick shell theory found in (3) as follows:(16)
e
α
β
=
1
2
∇
α
v
+
∇
β
u
-
b
α
β
w
,
k
α
β
=
∇
α
b
β
ρ
u
ρ
+
b
α
ρ
∇
β
u
ρ
+
b
β
ρ
∇
α
u
ρ
+
∇
α
∇
β
w
-
b
α
ρ
b
ρ
β
w
,
Q
α
β
=
0.5
b
α
μ
∇
β
b
μ
ρ
u
ρ
+
b
α
μ
b
μ
ρ
∇
β
u
ρ
+
b
β
μ
b
μ
ρ
∇
α
u
ρ
+
b
β
μ
∇
α
b
μ
ρ
u
ρ
+
b
α
μ
∇
β
∇
μ
w
+
b
β
μ
∇
α
∇
μ
w
.The membrane deformation tensore
α
β and bending tensor k
α
β are widely used in thin shell theory. But the tensor Q
α
β is neglected [24]. In classical thin shell and thick shell computation analysis, the contribution of this tensor (Gauss curvature tensor) to the deformation energy disappears in the stiffness matrix. As the shell becomes thicker, the contribution of Gauss curvature tensor in terms of energy is no more negligible as compared to that of the first two tensors used in classical thin shells. So this model is able to handle thick shells as well as the classical R-M models.We express the strain components using the parameters of the cylindrical shell’s midsurface. In detail, for the case of a cylindrical shell, we have(17)ex=u1,x,eφ=1Ru2,φ+w,2exφ=1Ru1,φ+u2,x,kx=w,xx,kφ=-1R22u2,φ-w,φφ+w,2kxφ=-2R-w,xφ+u2,x,Qx=0,Qφ=-1R3u2,φ-w,φφ,2Qxφ=-1R2u2,x-w,xφ.


### 2.6. Interpolation

The approach we are using here is based on strain interpolation. This method consists in solving a system of differential equations from strain assumption. We have to find the displacements functions which perfectly capture the rigid-like motion; then particular displacements are calculated from approximation of the deformations components. These components are approximated such that membrane, shearing, and bending behavior are decoupled. So a pure bending state or a pure shearing state can be well represented. In this section we consideru
1
=
u and u
2
=
v
.The rigid-body-like motion is defined forε
α
β
=
0. It means that e
α
β
=
0, k
α
β
=
0, and Q
α
β
=
0.(18)ux,φ=uφ.Let g
(
φ
) be an arbitrary function that can satisfy the following equations:(19)vx,φ=-1Ru′φ·x+gφ,wx,φ=+1Ru′′φ·x-g′φ,(20)u3φ+u′φ=0,(21)u4φ+u2φ=0,(22)g3φ+g1φ=0.Since our displacement field is of six components, three rotations and three translations, we need six independent parameters to define our rigid body-like motion [18, 25].(23)uφ=a1+a2Rcos⁡φ-cos⁡θ-a3Rsin⁡φ,(24)gφ=a4Rcos⁡φcos⁡θ-1-a5sin⁡φ+a6cos⁡φ.When (23) and (24) are taken into (18) and (19), it comes that the rigid body-like motion is well handled by the following displacement components:(25)u0=a1+a2Rcos⁡φ-cos⁡θ-a3Rsin⁡φ,v0=a2xsin⁡φ+a3xcos⁡φ+a4Rcos⁡φcos⁡θ-1-a5sin⁡φ+a6cos⁡φ,w0=-a2xcos⁡φ+a3xsin⁡φ+a4Rsin⁡φcos⁡θ+a5cos⁡φ+a6sin⁡φ.The displacements found in (25) are expressed in matrix form:(26)U0=P0A0,where U
0
t
=
u
0
v
0
w
0. (27)A·t0=a1a2a3a4a5a6.
A
0 is the vector of unknown coefficient which captures the rigid body-like motion:(28)P0=1Rcos⁡φ-cos⁡θ-Rsin⁡φ0000xsin⁡φxcos⁡φRcos⁡φcos⁡φ-1-sin⁡φcos⁡φ0-xcos⁡φxsin⁡φRsin⁡φcos⁡θcos⁡φsin⁡φ.
P
0 is the matrix of interpolation functions.The general solution displacements fieldU is the sum of the displacements field of rigid body-like motion and that of a particular solution of the deformation U
d:(29)U=U0+Ud=uvw·t,where U
·
t
d
=
u
d
v
d
w
d
.The strain components are interpolated in such a way that interference between shear strains and bending is not possible. The following displacement components satisfy strains defined above:(30)ud=a7Rx+a8Rφ+a9Rxφ,vd=a10φ+a11xφ,wd=a12x2+a13xφ+a14φ2+a15x3+a16x2φ+a17xφ2+a18φ3,Ud=PdAd,with(31)Ad·t=a7a8a9a10a11a12a13a14a15a16a17a18,Pd=RxRφRxφ000000000000φxφ000000000000x2xφφ2x3x2φxφ2φ3.
P
d is the interpolation functions matrix related to the particular solution of displacements field.From the total displacement vector field expression above, we have 18 unknowns. It then requires 18 displacement parameters to find all the above unknowns (see Figure2).Figure 2
Triangular single-curved element.So the total displacement vector over the triangular element is(32)U^e·t=u^1v^1w^1θ^x1θ^φ1γ^1u^2v^2w^3θ^x2θ^φ2γ^2u^3v^3w^3θ^x3θ^φ3γ^3.In the local (curvilinear) system of coordinates, each node has six (06) components of displacement.


### 2.7. Stiffness Matrix

In order to calculate the stiffness matrix, we have to formulate the variational problem over a domain. LetS be the border of the domain; then let ∂
S
=
γ
0
∪
γ
1 be partitioned in two and the border of the shell ∂
Ω
=
Γ
0
∪
Γ
1 with Γ
0
=
γ
0
×
-
h
/
2
,
h
/
2, and Γ
1
=
γ
1
×
-
h
/
2
,
h
/
2
∪
Γ
-
∪
Γ
+. We define Γ
-
=
S
×
-
h
/
2 and Γ
+
=
S
×
h
/
2. We suppose here that the shell is clamped on Γ
0 and loaded by volume and surface forces as stated above; the three-dimensional variational equation related to the equilibrium is(33)findU∈IHΓ01∫Ω·σ-~TU:ε-~VdΩ=∫Ω·f·VdΩ+∫∂Ω·g·VdГ=Lv,forV∈IHΓ01.
f is volume forces in the domain Ω and g is surface forces in the domain ∂
Ω.The constitutive law of the linear elastic homogenous material is(34)σ-~=c-̿ε-~,where c
-
̿
=
C
i
j
k
l
=
λ
g
i
j
g
k
l
+
μ
g
i
k
g
j
l
+
g
i
l
g
j
k
,
λ and μ are Lamé coefficients which depend on intrinsic properties of materials.Replacing (34) in (33), the problem is stated as follows:(35)findu∈IHΓ01∫Ω·λgijgkl+μgikgjl+gilgjkε-~u:ε-~vdΩ=∫Ω·f·vdΩ+∫∂Ω·g·vdГ=Lv,forv∈IHΓ01.The strain tensor is defined at (7), if substitution is done in the left side of the above equation, it can now be written:(36)Au,v=∫Ω·λgαβgγρ+μgαγgβρ+gαρgβγ·eαβu-zkαβu+z2Qαβu:eγρv-zkγρv+z2QγρvdΩ=Lv.Naturally ε
i
3
(
U
)
=
0 at the limit analysis. Then at the midsurface, g
α
β
=
(
μ
-
1
)
ρ
α
(
μ
-
1
)
λ
β
a
ρ
λ. By using Taylor’s expansion on (
μ
-
1
)
ρ
α, (37)μ-1ρα=δρα+zbρα+∑n≥2∞znbnραand truncating at n
=
1, we obtain the best first-order two-dimensional (38)A1u,v=Eh1-v2·∫S·1-veαβu+veλλaαβeαβvdS+Eh3121-v2·∫S.1-vkαβu+vkλλaαβkαβvdS+Eh3121-v2·∫S·1-veαβu+veλλaαβQαβvdS+Eh3121-v2·∫S·1-vQαβu+vQλλaαβeαβvdS+Eh5801-v2·∫S·1-vQαβu+vQλλaαβQαβvdS=Lv.Remark 1.
For Taylor’s expansion used above, if we truncate (37) for n
=
0, the bilinear continuous form A
u
,
v will take the value:(39)A0u,v=Eh1-v2·∫S·1-veαβu+veλλaαβeαβvdS+Eh3121-v2·∫S·1-vkαβu+vkλλaαβkαβvdS,which is found in the literature to be K-L shell deformation energy. It is included in A
1
u
,
v which is the N-T shell deformation energy. In what follows, we compute both K
G
N
-
T and K
G
R
-
M which are, respectively, stiffness matrix related to A
1
u
,
v and A
0
u
,
v.For the computing aim,A
1
u
,
v is transformed (40)A1u,v=∫SEvTCGEudS=∫SEvTCGEuJdξdη.For this purpose, components of strain tensor are expressed as a product of the matrix of interpolation functions’ gradient and that of the total displacement vector over the domain. From (39), the stress tensor components are deduced in function of the product above.(41)εm=eαβ=DmA=DmP-1U^=BmU^,Bm=DmP-1,σm=cεm=cBmU^.Also, the bending stress and strain are expressed:(42)εb=kαβ=DbA=DbP-1U^=BbU^,Bb=DbP-1,σb=cεb=cBbU^.Finally as previous, the Gauss curvature tensor is expressed as(43)εg=Qαβ=DgA=DgP-1U^=BgU^,Bg=DgP-1,σg=cεg=cBgU^.Total displacement is stated as the sum of particular solution and homogeneous solution of our system of strain differential equations.(44)U=U0+Ud=P0A0+PdAd=P0PdA0Ad=N·A,where N
=
P
0
P
d is the matrix of interpolation functions and(45)U=N·A=N·P-1U^.The left-hand side of the variational problem reads(46)
A
1
u
,
v
=
1
2
∫
Ω
ε
m
t
c
t
ε
m
-
z
ε
m
t
c
t
ε
b
+
z
2
ε
m
t
c
t
ε
g
-
z
ε
b
t
c
t
ε
m
+
z
2
ε
b
t
c
t
ε
b
-
z
3
ε
b
t
c
t
ε
g
+
z
2
ε
g
t
c
t
ε
m
-
z
3
ε
g
t
c
t
ε
b
+
z
4
ε
g
t
c
t
ε
g
d
Ω
=
L
v
.This is written as follows:(47)A1u,v=12V^t∫ShBmtCtBmRdφdxU^-12V^t∫S∇h28BmtCtBbRdφdxU^+12V^t∫Sh312BmtCtBgRdφdxU^-12V^t∫S∇h28BbtCtBmRdφdxU^+12V^t∫Sh312BbtCtBbRdφdxU^-12V^t∫S∇h432BbtCtBgRdφdxU^+12V^t∫Sh312BgtCtBmRdφdxU^-12V^t∫S∇h432BgtCtBbRdφdxU^+12V^t∫Sh580BgtCtBgRdφdxU^=Lv.One can deduct from (47) what follows: (48)A1u,v=V^tKGN-TU^,A0u,v=V^tKGR-MU^,where K
G
(
·
) is the global stiffness matrix. The midsurface is symmetric in shell thickness: that is, the midsurface is at +
h
/
2 from the top surface and at -
h
/
2 from the bottom one.(49)
K
G
N-T
=
∫
S
h
B
m
t
C
t
B
m
+
h
3
12
B
b
t
C
t
B
b
+
h
3
12
B
m
t
C
t
B
g
+
h
3
12
B
g
t
C
t
B
m
+
h
5
80
B
g
t
C
t
B
g
R
d
φ
d
x
,
K
G
R-M
=
∫
S
h
B
m
t
C
t
B
m
+
h
3
12
B
b
t
C
t
B
b
R
d
φ
d
x
.The right-hand side of the best first-order two-dimensional variational equation reads(50)Lv=V^t∫ΩP-1·tNtfdΩ+V^t∫∂ΩP-1·tNtgdГ=V^tFV+M=0,where f and g are defined as in [24] and F is the force vector from distributed load F
V
+
M.(51)FV=∫S·NtfeRdφdx,M=∫S·NtgRdφdx.By equating A
1
u
,
v to L
(
v
), we obtain the structural equilibrium equation over the whole area as follows:(52)KGU^=F.Over a single-curved triangular element of area s
e the elementary stiffness matrix K
e and elementary force vectors F
e can be locally implemented.(53)
K
e
=
∫
s
e
h
B
-
m
t
C
t
B
-
m
+
h
3
12
B
-
m
t
C
t
B
-
g
+
h
3
12
B
-
b
t
C
t
B
-
b
+
h
3
12
B
-
g
t
C
t
B
-
m
+
h
5
80
B
-
g
t
C
t
B
-
g
R
d
φ
d
x
.By regular assembling process of elementary stiffness matrix K
e and elementary force vectors F
e, we then expressed the global stiffness matrix and force vector over the whole area as follows:(54)KG=∑e=1NKe,F=∑e=1NFe.
N is total number of elements.


## 2.1. N-T’s Two-Dimensional Elastic Thick Shells Model

The N-T’s two-dimensional model for linear elastic thick shells has been deduced from the three-dimensional problem without any ad hoc assumption whether of geometrical or mechanical nature. The two-dimensional equations are deduced by applying asymptotic analysis on a family of variational equations obtained from an abstract scaled shell through multiple scaling of the initial three-dimensional equations. The theoretical elastic thick shells governing equations and limits displacements obtained from this approach (at the limit analysis) are more general because they contain additional terms to those found in Kirchhoff-Love model. The model also completes the thick shells theory of Reissner-Mindlin [24].For a three-dimensional shell, the following basesg
1
,
g
2
,
g
3
,
g
1
,
g
2
,
g
3 are dual bases and a
1
,
a
2
,
a
3
,
a
1
,
a
2
,
a
3 defined on the midsurface also constitute dual bases of shell if χ
<
1. Let U
i
x
1
,
x
2
,
z be the three-dimensional displacement field vector; it is defined as follows in the local coordinate system: (1)U=Uαgα+U3g3=u-αaα+u-3a3,u-ρ=u-ρx1,x2,z,u-3=u-3x1,x2,z.The strain tensor for a three-dimensional shell reads(2)εαβU=12Uα/β+Uβ/α=εαβu-=μαν∇βu-ν-bαβu-3+μβν∇αu-ν-bνβu-32,(3)εα3U=12Uα/3+U3/α=εα3u-=μανu-ν,3+u-3,α+bανu-ν2,(4)ε33U=U3/3=U3,3=u-3,3=ε33u-.From the limit analysis N-T demonstrated that the displacements field satisfies the equation ε
i
3
(
U
)
=
0 and the unique solution is(5)u-ρx1,x2,z=μργuγx1,x2-z∂ρw,W=w=u3,Uαx1,x2,z=μαρu-ρ,μργ=δργ-zbργ,Uαx1,x2,z=uαx1,x2-zθαx1,x2+z2ψαx1,x2,where (
u
α
) are local membrane displacement components and w
=
u
3 the transvers displacement:(6)θα=2bατuτ+∇αw,ψα=bαγbγτuτ+bαγ∇γw.The general shell deformation tensor derived from the kinematics defined above is(7)εαβU=eαβu-zkαβu+z2Qαβu,(8)εi3U=0,where e
α
β is the membrane strain tensor, k
α
β is the bending strain tensor, and Q
α
β is the tensor of Gauss curvature.The strain tensor contains the change of the third fundamental form:(9)Qαβ=12cαβ-Cαβ,c=cαβ=bαμbμβ=bαμbβμ,Cαβ=BαμBμβ=BαμBβμ.In Kirchhoff-Love thin shells theory, the shell is assumed to be thin. Following that assumption, the three-dimensional displacements field vectorU
α
x
1
,
x
2
,
z is taken equal to the midsurface displacements field:(10)Uαx1,x2,z=uαx1,x2-z∇αwx1,x2,Wx1,x2,z=wx1,x2=u3.It should be noted that the terms ψ
α
=
b
α
τ
b
τ
ρ
u
ρ
+
b
α
τ
∇
τ
w in the current model of thick shells disappear in the classical theory of thin shells because it is proportional to χ
2 and considered small. The thin shell theory, with the above assumption, cannot perform good shells behavior because the thickness has a considerable influence on the shells behavior. The thick shell model used in this survey appears to be suitable to handle both thin and thick shells because of its completeness.


## 2.2. Curved Shell Elements

Curved shell elements are suitable to model the midsurface geometry more accurately. In the case of certain surfaces such as cylindrical shells, it means the exact description of the original surface. For more complicated cases, similar 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 [25].


## 2.3. Triangular Curve Shell Element

The triangular shell element is described in Figure1; it can be double-curved (e.g., spherical structures) or single-curved element (cylindrical structures). In accordance with the N-T’s equations of the theory of elastic thick shells, the geometrical properties of the cylindrical shell are the following:(11)OM→=X=Rcos⁡φ,Y=Rsin⁡φ,Z=x,where x and φ are curvilinear coordinates, 0
≤
x
≤
L and 0
≤
φ
≤
θ, R
≡ cylinder radius, and X
,
Y
,
Z are global coordinates.Figure 1
Triangular shell element.


## 2.4. Displacement Field

LetU
α be the global displacement field vector defined as follows: (12)Uα=uα-zθα+z2ψα,W=w.With respect to the cylinder’s geometry, the angles of rotations and the Gaussian curvatures are developed in function of membrane displacements u
1 and u
2 and the derivatives of the transverse displacement w
,
φ
,
w
,
x including curvature components b
α
τ:(13)θ1=2u1R1-w,xA1,θ2=2u2R2-w,φA2,ψα=bατbτρuρ+bατ∇τw,which leads to(14)ψ1=-1R1u1R1-1A1w,x=1R1b1τuτ+θ1,ψ2=-1R2u2R2-1A2w,φ=1R2b2τuτ+θ2,where A
1
,
A
2
,
R
1, and R
2 are geometrical properties of the midsurface of a shell within the respective coordinate-lines x and φ.In the case of a cylindrical shell they are defined as follows [25]:(15)A1=a1→=∂M∂x→=1;R1=∞,A2=a2→=∂M∂φ→=R;R2=R.The third component of rotation is added for a smooth transfer of element stiffness matrix from one system of coordinates to another.


## 2.5. Compatibility

We recall the strain components in thick shell theory found in (3) as follows:(16)
e
α
β
=
1
2
∇
α
v
+
∇
β
u
-
b
α
β
w
,
k
α
β
=
∇
α
b
β
ρ
u
ρ
+
b
α
ρ
∇
β
u
ρ
+
b
β
ρ
∇
α
u
ρ
+
∇
α
∇
β
w
-
b
α
ρ
b
ρ
β
w
,
Q
α
β
=
0.5
b
α
μ
∇
β
b
μ
ρ
u
ρ
+
b
α
μ
b
μ
ρ
∇
β
u
ρ
+
b
β
μ
b
μ
ρ
∇
α
u
ρ
+
b
β
μ
∇
α
b
μ
ρ
u
ρ
+
b
α
μ
∇
β
∇
μ
w
+
b
β
μ
∇
α
∇
μ
w
.The membrane deformation tensore
α
β and bending tensor k
α
β are widely used in thin shell theory. But the tensor Q
α
β is neglected [24]. In classical thin shell and thick shell computation analysis, the contribution of this tensor (Gauss curvature tensor) to the deformation energy disappears in the stiffness matrix. As the shell becomes thicker, the contribution of Gauss curvature tensor in terms of energy is no more negligible as compared to that of the first two tensors used in classical thin shells. So this model is able to handle thick shells as well as the classical R-M models.We express the strain components using the parameters of the cylindrical shell’s midsurface. In detail, for the case of a cylindrical shell, we have(17)ex=u1,x,eφ=1Ru2,φ+w,2exφ=1Ru1,φ+u2,x,kx=w,xx,kφ=-1R22u2,φ-w,φφ+w,2kxφ=-2R-w,xφ+u2,x,Qx=0,Qφ=-1R3u2,φ-w,φφ,2Qxφ=-1R2u2,x-w,xφ.


## 2.6. Interpolation

The approach we are using here is based on strain interpolation. This method consists in solving a system of differential equations from strain assumption. We have to find the displacements functions which perfectly capture the rigid-like motion; then particular displacements are calculated from approximation of the deformations components. These components are approximated such that membrane, shearing, and bending behavior are decoupled. So a pure bending state or a pure shearing state can be well represented. In this section we consideru
1
=
u and u
2
=
v
.The rigid-body-like motion is defined forε
α
β
=
0. It means that e
α
β
=
0, k
α
β
=
0, and Q
α
β
=
0.(18)ux,φ=uφ.Let g
(
φ
) be an arbitrary function that can satisfy the following equations:(19)vx,φ=-1Ru′φ·x+gφ,wx,φ=+1Ru′′φ·x-g′φ,(20)u3φ+u′φ=0,(21)u4φ+u2φ=0,(22)g3φ+g1φ=0.Since our displacement field is of six components, three rotations and three translations, we need six independent parameters to define our rigid body-like motion [18, 25].(23)uφ=a1+a2Rcos⁡φ-cos⁡θ-a3Rsin⁡φ,(24)gφ=a4Rcos⁡φcos⁡θ-1-a5sin⁡φ+a6cos⁡φ.When (23) and (24) are taken into (18) and (19), it comes that the rigid body-like motion is well handled by the following displacement components:(25)u0=a1+a2Rcos⁡φ-cos⁡θ-a3Rsin⁡φ,v0=a2xsin⁡φ+a3xcos⁡φ+a4Rcos⁡φcos⁡θ-1-a5sin⁡φ+a6cos⁡φ,w0=-a2xcos⁡φ+a3xsin⁡φ+a4Rsin⁡φcos⁡θ+a5cos⁡φ+a6sin⁡φ.The displacements found in (25) are expressed in matrix form:(26)U0=P0A0,where U
0
t
=
u
0
v
0
w
0. (27)A·t0=a1a2a3a4a5a6.
A
0 is the vector of unknown coefficient which captures the rigid body-like motion:(28)P0=1Rcos⁡φ-cos⁡θ-Rsin⁡φ0000xsin⁡φxcos⁡φRcos⁡φcos⁡φ-1-sin⁡φcos⁡φ0-xcos⁡φxsin⁡φRsin⁡φcos⁡θcos⁡φsin⁡φ.
P
0 is the matrix of interpolation functions.The general solution displacements fieldU is the sum of the displacements field of rigid body-like motion and that of a particular solution of the deformation U
d:(29)U=U0+Ud=uvw·t,where U
·
t
d
=
u
d
v
d
w
d
.The strain components are interpolated in such a way that interference between shear strains and bending is not possible. The following displacement components satisfy strains defined above:(30)ud=a7Rx+a8Rφ+a9Rxφ,vd=a10φ+a11xφ,wd=a12x2+a13xφ+a14φ2+a15x3+a16x2φ+a17xφ2+a18φ3,Ud=PdAd,with(31)Ad·t=a7a8a9a10a11a12a13a14a15a16a17a18,Pd=RxRφRxφ000000000000φxφ000000000000x2xφφ2x3x2φxφ2φ3.
P
d is the interpolation functions matrix related to the particular solution of displacements field.From the total displacement vector field expression above, we have 18 unknowns. It then requires 18 displacement parameters to find all the above unknowns (see Figure2).Figure 2
Triangular single-curved element.So the total displacement vector over the triangular element is(32)U^e·t=u^1v^1w^1θ^x1θ^φ1γ^1u^2v^2w^3θ^x2θ^φ2γ^2u^3v^3w^3θ^x3θ^φ3γ^3.In the local (curvilinear) system of coordinates, each node has six (06) components of displacement.


## 2.7. Stiffness Matrix

In order to calculate the stiffness matrix, we have to formulate the variational problem over a domain. LetS be the border of the domain; then let ∂
S
=
γ
0
∪
γ
1 be partitioned in two and the border of the shell ∂
Ω
=
Γ
0
∪
Γ
1 with Γ
0
=
γ
0
×
-
h
/
2
,
h
/
2, and Γ
1
=
γ
1
×
-
h
/
2
,
h
/
2
∪
Γ
-
∪
Γ
+. We define Γ
-
=
S
×
-
h
/
2 and Γ
+
=
S
×
h
/
2. We suppose here that the shell is clamped on Γ
0 and loaded by volume and surface forces as stated above; the three-dimensional variational equation related to the equilibrium is(33)findU∈IHΓ01∫Ω·σ-~TU:ε-~VdΩ=∫Ω·f·VdΩ+∫∂Ω·g·VdГ=Lv,forV∈IHΓ01.
f is volume forces in the domain Ω and g is surface forces in the domain ∂
Ω.The constitutive law of the linear elastic homogenous material is(34)σ-~=c-̿ε-~,where c
-
̿
=
C
i
j
k
l
=
λ
g
i
j
g
k
l
+
μ
g
i
k
g
j
l
+
g
i
l
g
j
k
,
λ and μ are Lamé coefficients which depend on intrinsic properties of materials.Replacing (34) in (33), the problem is stated as follows:(35)findu∈IHΓ01∫Ω·λgijgkl+μgikgjl+gilgjkε-~u:ε-~vdΩ=∫Ω·f·vdΩ+∫∂Ω·g·vdГ=Lv,forv∈IHΓ01.The strain tensor is defined at (7), if substitution is done in the left side of the above equation, it can now be written:(36)Au,v=∫Ω·λgαβgγρ+μgαγgβρ+gαρgβγ·eαβu-zkαβu+z2Qαβu:eγρv-zkγρv+z2QγρvdΩ=Lv.Naturally ε
i
3
(
U
)
=
0 at the limit analysis. Then at the midsurface, g
α
β
=
(
μ
-
1
)
ρ
α
(
μ
-
1
)
λ
β
a
ρ
λ. By using Taylor’s expansion on (
μ
-
1
)
ρ
α, (37)μ-1ρα=δρα+zbρα+∑n≥2∞znbnραand truncating at n
=
1, we obtain the best first-order two-dimensional (38)A1u,v=Eh1-v2·∫S·1-veαβu+veλλaαβeαβvdS+Eh3121-v2·∫S.1-vkαβu+vkλλaαβkαβvdS+Eh3121-v2·∫S·1-veαβu+veλλaαβQαβvdS+Eh3121-v2·∫S·1-vQαβu+vQλλaαβeαβvdS+Eh5801-v2·∫S·1-vQαβu+vQλλaαβQαβvdS=Lv.Remark 1.
For Taylor’s expansion used above, if we truncate (37) for n
=
0, the bilinear continuous form A
u
,
v will take the value:(39)A0u,v=Eh1-v2·∫S·1-veαβu+veλλaαβeαβvdS+Eh3121-v2·∫S·1-vkαβu+vkλλaαβkαβvdS,which is found in the literature to be K-L shell deformation energy. It is included in A
1
u
,
v which is the N-T shell deformation energy. In what follows, we compute both K
G
N
-
T and K
G
R
-
M which are, respectively, stiffness matrix related to A
1
u
,
v and A
0
u
,
v.For the computing aim,A
1
u
,
v is transformed (40)A1u,v=∫SEvTCGEudS=∫SEvTCGEuJdξdη.For this purpose, components of strain tensor are expressed as a product of the matrix of interpolation functions’ gradient and that of the total displacement vector over the domain. From (39), the stress tensor components are deduced in function of the product above.(41)εm=eαβ=DmA=DmP-1U^=BmU^,Bm=DmP-1,σm=cεm=cBmU^.Also, the bending stress and strain are expressed:(42)εb=kαβ=DbA=DbP-1U^=BbU^,Bb=DbP-1,σb=cεb=cBbU^.Finally as previous, the Gauss curvature tensor is expressed as(43)εg=Qαβ=DgA=DgP-1U^=BgU^,Bg=DgP-1,σg=cεg=cBgU^.Total displacement is stated as the sum of particular solution and homogeneous solution of our system of strain differential equations.(44)U=U0+Ud=P0A0+PdAd=P0PdA0Ad=N·A,where N
=
P
0
P
d is the matrix of interpolation functions and(45)U=N·A=N·P-1U^.The left-hand side of the variational problem reads(46)
A
1
u
,
v
=
1
2
∫
Ω
ε
m
t
c
t
ε
m
-
z
ε
m
t
c
t
ε
b
+
z
2
ε
m
t
c
t
ε
g
-
z
ε
b
t
c
t
ε
m
+
z
2
ε
b
t
c
t
ε
b
-
z
3
ε
b
t
c
t
ε
g
+
z
2
ε
g
t
c
t
ε
m
-
z
3
ε
g
t
c
t
ε
b
+
z
4
ε
g
t
c
t
ε
g
d
Ω
=
L
v
.This is written as follows:(47)A1u,v=12V^t∫ShBmtCtBmRdφdxU^-12V^t∫S∇h28BmtCtBbRdφdxU^+12V^t∫Sh312BmtCtBgRdφdxU^-12V^t∫S∇h28BbtCtBmRdφdxU^+12V^t∫Sh312BbtCtBbRdφdxU^-12V^t∫S∇h432BbtCtBgRdφdxU^+12V^t∫Sh312BgtCtBmRdφdxU^-12V^t∫S∇h432BgtCtBbRdφdxU^+12V^t∫Sh580BgtCtBgRdφdxU^=Lv.One can deduct from (47) what follows: (48)A1u,v=V^tKGN-TU^,A0u,v=V^tKGR-MU^,where K
G
(
·
) is the global stiffness matrix. The midsurface is symmetric in shell thickness: that is, the midsurface is at +
h
/
2 from the top surface and at -
h
/
2 from the bottom one.(49)
K
G
N-T
=
∫
S
h
B
m
t
C
t
B
m
+
h
3
12
B
b
t
C
t
B
b
+
h
3
12
B
m
t
C
t
B
g
+
h
3
12
B
g
t
C
t
B
m
+
h
5
80
B
g
t
C
t
B
g
R
d
φ
d
x
,
K
G
R-M
=
∫
S
h
B
m
t
C
t
B
m
+
h
3
12
B
b
t
C
t
B
b
R
d
φ
d
x
.The right-hand side of the best first-order two-dimensional variational equation reads(50)Lv=V^t∫ΩP-1·tNtfdΩ+V^t∫∂ΩP-1·tNtgdГ=V^tFV+M=0,where f and g are defined as in [24] and F is the force vector from distributed load F
V
+
M.(51)FV=∫S·NtfeRdφdx,M=∫S·NtgRdφdx.By equating A
1
u
,
v to L
(
v
), we obtain the structural equilibrium equation over the whole area as follows:(52)KGU^=F.Over a single-curved triangular element of area s
e the elementary stiffness matrix K
e and elementary force vectors F
e can be locally implemented.(53)
K
e
=
∫
s
e
h
B
-
m
t
C
t
B
-
m
+
h
3
12
B
-
m
t
C
t
B
-
g
+
h
3
12
B
-
b
t
C
t
B
-
b
+
h
3
12
B
-
g
t
C
t
B
-
m
+
h
5
80
B
-
g
t
C
t
B
-
g
R
d
φ
d
x
.By regular assembling process of elementary stiffness matrix K
e and elementary force vectors F
e, we then expressed the global stiffness matrix and force vector over the whole area as follows:(54)KG=∑e=1NKe,F=∑e=1NFe.
N is total number of elements.


## 3. Numerical Study

A numerical study which aims to quantitatively investigate the accuracy of different thick shell theories for the linear elastic analysis of cylindrical structures and to provide a basis for the long-term analysis is presented.One of the problems frequently used to evaluate the performances of a shell element is that of the cylindrical roof (of Scordelis-lo) [4, 26] subjected to its self-weight described in Figure 3. The flat rims are free and the curved edges rest on rigid diaphragms in their plans. The geometrical and mechanical characteristics are indicated for h
/
R = 0.01, L
/
h = 200. The length of the cylinder L = 600 cm; its radius is R = 300 cm; the angle subtended by the roof is 2
θ = 80 deg; the thickness is h
=
3 cm; Young’s modulus is E
=
3.10
e
+
10 Pa and Poisson’s ratio is ν
=
0.Figure 3
Definition of Scordelis-Lo problem. Cylindrical roof subjected to its self-weight.


### 3.1. Convergence

The quarter of the roof is discretized by considering triangular single-curved elements withN
=
4,6
,
… elements on edges A
B and A
D. The transverse displacements at B and C are plotted (Figure 4) and compared with other models of finite elements: see Figure 5.Figure 4
The undeformed and deformed configurations of one-quarter of the roof.Figure 5
Displacement-convergence plot for Scordelis-Lo Roof. (a) Convergence of transverse displacements at pointC (W
C). (b) Convergence of transverse displacements at point B (W
B).
(a)
(b)The computed deformed limit surface is shown in Figure4. Here, displacement-convergence results are shown in Figure 5 and Table 1. In the diaphragm case, we monitor the displacements under the self-weight at points B and C. The convergence properties of the method are evident from Figure 5 and Table 1. As may be seen, the CSFE3 converges as well as both the semifinite elements (SFE) and the DKT elements for the displacements under the loads.Table 1
Comparison of transverse displacements atB and C.
Mesh step
W
C (cm) reference: 0.541
W
B (cm) reference: −3.61
DKT18
DKT12
SFE3
CSFE3
DKT18
DKT12
SFE3
CSFE3
N
=
2
0.269
0.175
0.072
0.210
-
2.893
-
3.773
-
2.407
−0.470
N
=
4
0.330
0.330
0.292
0.060
-
2.436
-
2.499
-
2.541
−1.640
N
=
6
0.409
0.420
0.393
0.210
-
2.848
-
2.882
-
2.951
−2.460
N
=
8
0.455
0.465
0.445
0.360
-
3.106
-
3.129
-
3.206
−2.940
N
=
10
0.482
0.489
0.475
0.460
-
3.259
-
3.275
-
3.363
−3.300
N
=
12
0.498
0.504
0.493
0.510
-
3.354
-
3.365
-
3.462
−3.480
N
=
14
0.512
0.516
0.510
0.526
-
3.438
-
3.446
-
3.555
−3.580


### 3.2. Scaling and Variation of Displacements

A successive scaling is carried out on the thickness ratioh
/
R over the Scordelis-lo roof. We have used the following range of ratios: 0.1; 0.2; 0.3; 0.325; 0.4; 0.5, because it certainly belongs to both thin shells and thick shells ranges of thickness. The radius R is constant while the thickness h varies with the ratio. The scaled shell obtained here has been implemented using both K-L theory and that proposed by N-T. The results are plotted for each ratio as shown in Figure 6.Figure 6
(a) Variation of transverse displacements atB for h
/
R = 0.1. (b) Variation of transverse displacements at B for h
/
R = 0.2. (c) Variation of transverse displacements at B for h
/
R = 0.3. (d) Variation of transverse displacements at B for h
/
R = 0.325. (e) Variation of transverse displacements at B for h
/
R = 0.4. (f) Variation of transverse displacements at B for h
/
R = 0.5.
(a)
(b)
(c)
(d)
(e)
(f)We observe from the diagrams shown in Figure6 that for the thickness ratios h
/
R = 0.1, 0.2, and 0.3, the transverse displacements at point B are the same for single-curved triangular computation of both K-L and N-T models (see Figures 6(a), 6(b), and 6(c)). One can verify that, at any other points of the roof structure, these same observations are satisfied. For thickness ratios h
/
R greater or equal to 0.325, we observe that transverse displacements computed from the shells equations of K-L and N-T are not the same (see Figures 6(d), 6(e), and 6(f)). We can still investigate the exact value h
/
R in 3
/
10
13
/
40 from which displacements results issued from these two approaches of shells theories diverge in a cylindrical wall and probably other shells.


### 3.3. Deviations

N-T’s equations for elastic thick shell are more general [24] and with regard to that, we assume the displacements issued from this model of thick shell to be more accurate. We can then investigate the deviations encountered when using K-L model to evaluate a shell. The investigation consists in calculating the deviations (difference between K-L displacement and N-T displacement) for each mesh-step (number of elements per half side) with respect to the series of thickness ratios h
/
R = 0.1; 0.2; 0.3; 0.325; 0.4; and 0.5 at the same point. This investigation was carried out on two different points C and B (see Figure 3) and the results are plotted in Figures 7(a) and 7(b).Figure 7
(a) Evolution of displacements deviations at pointC. (b) Evolution of displacements deviations at point B.
(a)
(b)From diagrams shown in Figures7(a) and 7(b), we observe that (i) at points B and C of our cylindrical roof, deviations are encountered from thickness ratios h
/
R above 0.3; (ii) deviations increase with mesh-steps at point C and decrease at point B.


## 3.1. Convergence

The quarter of the roof is discretized by considering triangular single-curved elements withN
=
4,6
,
… elements on edges A
B and A
D. The transverse displacements at B and C are plotted (Figure 4) and compared with other models of finite elements: see Figure 5.Figure 4
The undeformed and deformed configurations of one-quarter of the roof.Figure 5
Displacement-convergence plot for Scordelis-Lo Roof. (a) Convergence of transverse displacements at pointC (W
C). (b) Convergence of transverse displacements at point B (W
B).
(a)
(b)The computed deformed limit surface is shown in Figure4. Here, displacement-convergence results are shown in Figure 5 and Table 1. In the diaphragm case, we monitor the displacements under the self-weight at points B and C. The convergence properties of the method are evident from Figure 5 and Table 1. As may be seen, the CSFE3 converges as well as both the semifinite elements (SFE) and the DKT elements for the displacements under the loads.Table 1
Comparison of transverse displacements atB and C.
Mesh step
W
C (cm) reference: 0.541
W
B (cm) reference: −3.61
DKT18
DKT12
SFE3
CSFE3
DKT18
DKT12
SFE3
CSFE3
N
=
2
0.269
0.175
0.072
0.210
-
2.893
-
3.773
-
2.407
−0.470
N
=
4
0.330
0.330
0.292
0.060
-
2.436
-
2.499
-
2.541
−1.640
N
=
6
0.409
0.420
0.393
0.210
-
2.848
-
2.882
-
2.951
−2.460
N
=
8
0.455
0.465
0.445
0.360
-
3.106
-
3.129
-
3.206
−2.940
N
=
10
0.482
0.489
0.475
0.460
-
3.259
-
3.275
-
3.363
−3.300
N
=
12
0.498
0.504
0.493
0.510
-
3.354
-
3.365
-
3.462
−3.480
N
=
14
0.512
0.516
0.510
0.526
-
3.438
-
3.446
-
3.555
−3.580


## 3.2. Scaling and Variation of Displacements

A successive scaling is carried out on the thickness ratioh
/
R over the Scordelis-lo roof. We have used the following range of ratios: 0.1; 0.2; 0.3; 0.325; 0.4; 0.5, because it certainly belongs to both thin shells and thick shells ranges of thickness. The radius R is constant while the thickness h varies with the ratio. The scaled shell obtained here has been implemented using both K-L theory and that proposed by N-T. The results are plotted for each ratio as shown in Figure 6.Figure 6
(a) Variation of transverse displacements atB for h
/
R = 0.1. (b) Variation of transverse displacements at B for h
/
R = 0.2. (c) Variation of transverse displacements at B for h
/
R = 0.3. (d) Variation of transverse displacements at B for h
/
R = 0.325. (e) Variation of transverse displacements at B for h
/
R = 0.4. (f) Variation of transverse displacements at B for h
/
R = 0.5.
(a)
(b)
(c)
(d)
(e)
(f)We observe from the diagrams shown in Figure6 that for the thickness ratios h
/
R = 0.1, 0.2, and 0.3, the transverse displacements at point B are the same for single-curved triangular computation of both K-L and N-T models (see Figures 6(a), 6(b), and 6(c)). One can verify that, at any other points of the roof structure, these same observations are satisfied. For thickness ratios h
/
R greater or equal to 0.325, we observe that transverse displacements computed from the shells equations of K-L and N-T are not the same (see Figures 6(d), 6(e), and 6(f)). We can still investigate the exact value h
/
R in 3
/
10
13
/
40 from which displacements results issued from these two approaches of shells theories diverge in a cylindrical wall and probably other shells.


## 3.3. Deviations

N-T’s equations for elastic thick shell are more general [24] and with regard to that, we assume the displacements issued from this model of thick shell to be more accurate. We can then investigate the deviations encountered when using K-L model to evaluate a shell. The investigation consists in calculating the deviations (difference between K-L displacement and N-T displacement) for each mesh-step (number of elements per half side) with respect to the series of thickness ratios h
/
R = 0.1; 0.2; 0.3; 0.325; 0.4; and 0.5 at the same point. This investigation was carried out on two different points C and B (see Figure 3) and the results are plotted in Figures 7(a) and 7(b).Figure 7
(a) Evolution of displacements deviations at pointC. (b) Evolution of displacements deviations at point B.
(a)
(b)From diagrams shown in Figures7(a) and 7(b), we observe that (i) at points B and C of our cylindrical roof, deviations are encountered from thickness ratios h
/
R above 0.3; (ii) deviations increase with mesh-steps at point C and decrease at point B.


## 4. Discussion

We have investigated in this framework the influence of the third fundamental form proportional toχ
2 (with χ
2
=
h
/
2
R
2) included in N-T’s model of elastic thick shells. This model was implemented using cylindrical shell triangular finite elements (CSFE3). Based on the results obtained, the following comments can be made.The convergence studies of CSFE3 show that this curved finite element of three nodes only converges as well as SFE, DKT12, and DKT18 which are robust and greedy (memory) [27]. The last three finite elements models are based on theoretical approach of Kirchhoff-Love (thin shells) or Reissner-Mindlin (Thick shells) which neglect the third fundamental form in their shell kinematic equations. In addition to that, these classical shell theories are based on two powerful assumptions of K-L or R-M [28, 29].The theoretical equations of thick and thin shells issued from these assumptions can be questioned because the transverse fiber (normal to the reference surface) has been forced to follow a specific behavior without any scientific or technical background that justifies their kinematic equations [20, 24, 26]. Consequently some quantities without evident mechanical significations are added in the formulation of finite elements models (like DKT, SFE …) in order to improve computational accuracy [7, 26].In the contrary, our finite element (CSFE3) is based on N-T’s theoretical model of thick shells which is deduced fromε
α
3
(
u
) = 0 obtained at limits analysis of 3-dimensional shell equations [24]. Resultant kinematic equations of shell here are more general since they contain additional terms to those found in classical equations of shells. They are mathematically and technically justified without ad hoc assumption on transverse fiber’s behavior. This can explain why a simple 3-node curved triangular finite element obtained efficient convergent result. The investigation with a 6-node curved triangle will surely improve the convergence rate [7, 26].The successive scaling over the benchmark structure has presented divergence of displacements results ash
/
R increases for K-L kinematic model of shells and that of N-T. This is due to the influence of the additional terms found in N-T’s shell theory in the deformation energy of the structure. For homogeneous and isotropic elastic shells with a midsurface between -
h
/
2 and +
h
/
2, the total deformation energy contains membrane stiffness K
m, bending stiffness K
b, coupled stiffnesses K
m
g and K
g
m (membrane and Gaussian bending), and finally Gaussian stiffness K
g.For Reissner-Mindlin and K-L classical model of shells, the global stiffness matrixK
G
R-M contains only K
m and K
b while in N-T’s model, the global stiffness matrix K
G
N-T contains additional terms (K
m
g
+
K
g
m
+
K
g) to those found in K
G
R-M. The proportion of energy contributed by each stiffness matrix to the global deformation energy is as follows: K
m
=
α
m
E
h
/
R; K
b
=
α
b
E
h
/
R
h
/
R
2; K
m
g
=
α
m
g
E
h
/
R
10
-
1
h
/
R
2; and K
g
=
α
g
E
h
/
R
10
-
1
h
/
R
4. We monitored that for thickness ratio h
/
R less or equal to 0.3, that is, (
2
χ
)
2
<
1
/
10, the couple stiffness K
m
g and Gaussian stiffness K
g disappear in the global deformation energy because they are inversely proportional, respectively, to 1000 and 10.000. We then obtain K
G
N-T
≈
K
G
R-M. This explains why deviations are zero for thickness ratio less than 0.3 (see Figure 7). For thickness ratio greater or equal to 1
/
10 displacements results encountered for the two models are different because the coupled stiffness K
m
g and the Gaussian stiffness K
g which are now inversely proportional, respectively, to 100 and 1000 cannot disappear anymore in the global deformation energy of the structure. This additive proportion of energy into the global deformation energy found in this model improves total stiffness and increases the rigidity of the structure.


## 5. Conclusion

A finite element model for predicting the elastic thick shell behavior of cylindrical structure has been developed in this work. It has been quantitatively shown that the thickness ratio (thickness on radius) plays important role in their elastic behavior and structural safety. Hence, a reliable control of this parameter is required. N-T’s theoretical approach for the modelling of the displacement and stains in cylindrical structures have been examined and compared with that of K-L for the case of self-weight loading cylindrical roof.The comparison has revealed that N-T’s shell theory is suitable for the analysis of cylindrical thin and thick structures from the perspectives of both accuracy and simplicity. A 3-node cylindrical shell finite element has been developed for the simulation of static behaviour of shells in general, which requires the ability of well handling thickness ratios from the vicinity of zero to one.The ability of the CSFE3 model to deal with a wide range of ratios has been demonstrated through numerical examples.It has been found that the use of K-L theoretical approach based on those assumptions leads to inaccurate results as the thickness becomes greater than the square root of 0.1.It has also been found that above the value square root 0.1, the deviations linearly increase with the increase of the thickness ratio.The issues addressed in this work highlight the need for N-T’s approach to be implemented for the elastic shells analysis of cylindrical structures.The analytical model and the numerical results presented here contribute to the establishment of a foundation of theoretical knowledge required for reliable analysis, effective design, and safe use of cylindrical walled structures.Finally, the paper provides a platform for double-curved finite elements model to be developed for studying the influence of thickness ratio on all shell geometries.



---

*Source: 8936075-2016-08-29.xml*