Asymptotically accurate and locking-free finite element

implementation of the refined shell theory

Khanh Chau Lea,b1 and Hoang-Giang Buic
arXiv:2503.23369v1 [math.NA] 30 Mar 2025

a
Mechanics of Advanced Materials and Structures, Institute for Advanced Study in
Technology, Ton Duc Thang University, Ho Chi Minh City, Vietnam
b
Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Vietnam
c
Institute of Material Systems Modeling, Helmholtz-Zentrum Hereon, Geesthacht,
Germany

Abstract
A formulation of the 2D refined shell theory incorporating transverse shear
in the rescaled coordinates and angles of rotation is considered. This novel
approach provides the first asymptotically accurate and inherently locking-
free finite element implementation. Numerical simulations of semi-cylindrical
shells demonstrate excellent agreement between the analytical solution, the
2D refined shell theory, and three-dimensional elasticity theory, validating
the effectiveness and accuracy of the method.
Keywords: refined theory, shells, finite element, asymptotic accuracy,
locking-free.

1. Introduction
Shells, thin-walled structures with curved surfaces, are ubiquitous in en-
gineering, from civil and environmental to mechanical and aerospace appli-
cations. Understanding their structural behavior is crucial in designing safe,
efficient and resilient structures. When a shell’s thickness is significantly
smaller than its characteristic radius of curvature and longitudinal size, its
deformation can be effectively approximated using functions defined on 2D
surface coordinates. Asymptotic analysis, using shell thickness as a model
parameter, provides a hypothesis-free and systematic way to derive such 2D

1
Corresponding author. Phone: +84 93 4152458, email: [email protected]

Preprint submitted to Elsevier April 1, 2025

shell theories from 3D elasticity theory. This rigorous approach, based on
the asymptotic analysis of 3D elasticity’s strong [1] and weak formulation [2],
yields not only the classical shell theory [3, 4] but also various refined shell
theories (see [1, 2, 5]).
Accurate modeling of shells, especially concerning transverse shear ef-
fects, has driven the development of refined theories. Berdichevsky’s pioneer-
ing work constructed an asymptotically exact refined shell theory using the
variational-asymptotic method (VAM) [6, 7, 2]. This theory maintains the
accuracy up to O(h/R), O(h/l), and O(h2 /l2 ), where h, R, and l represent
the shell thickness, curvature radius, and longitudinal deformation scale, re-
spectively. While sharing asymptotic equivalence with Reissner’s first-order
shear deformation theory (FSDT) for plates [8], Berdichevsky’s formulation
uniquely provides pointwise accuracy for both displacements and stresses,
contrasting with Reissner’s integral-based accuracy. The necessity of refined
theories is evident from the potential inaccuracies in displacement predic-
tions by classical shell theory [5]. Subsequent efforts have expanded VAM
applications to laminated structures [9, 10, 11, 12], often optimizing FSDT
parameters for near-asymptotic correctness. Beyond linear elasticity, VAM
has been applied to various nonlinear materials [13, 14, 15, 16]. Recent inves-
tigations have explored asymptotically exact FSDT for functionally graded
plates [17], alongside VAM’s utility in dimension reduction, homogenization,
and dynamic analyses [18, 19, 20, 21, 22, 23, 24, 25].
The mathematical complexity of 2D refined shell theories often precludes
analytical solutions, making finite element implementation crucial for practi-
cal problems. However, membrane and shear locking can significantly reduce
the reliability [26, 27, 28]. Shear locking arises from the disparity between
shear and bending stiffnesses, coupled with much smaller rotations from pure
shear versus bending-induced curvature changes. This discrepancy, evident
in the variational-asymptotic analysis of the energy functional [2, 29], leads
to numerical instability as shell thickness vanishes (h → 0) with standard
low-order elements. Membrane locking stems from the contrast between
extension and bending stiffnesses and measures, causing instability due to
multiplication of small and large quantities as h → 0.
While sophisticated methods exist to alleviate locking, they often compro-
mise computational efficiency. Reduced and selective integration [30, 31, 32],
using different rules for different energy components, is a popular technique
but can introduce instability due to rank deficiency and spurious energy
modes [33]. Thus, alternative formulations and techniques have been devel-

2
oped for improved accuracy and stability [34, 35, 36, 37, 28, 38, 39, 40, 41, 42].
Recently, a novel rescaled formulation of FSDT for plates, inherently
free from shear locking regardless of discretization or integration, was pre-
sented [43]. Numerical simulations demonstrated the asymptotic accuracy
of Berdichevsky’s theory for plates using isogeometric elements ensuring C 1 -
continuity of primary variables. Building on this foundation, this work ex-
tends and pursuits two main goals. On one hand, we extend the rescaled for-
mulation to shells, providing a thickness-independent, locking-free approach
that boosts efficiency without high-order interpolation or complex integra-
tion. On the other hand, and most importantly, we aim for asymptotic
accuracy in the first finite element implementation of inherently locking-
free Berdichevsky’s refined shell theory. Following [6, 17], this requires C 1 -
continuity for both displacements and rotation angles. Numerical simulations
of semi-cylindrical shells, compared against analytical and 3D elasticity so-
lutions, validate the achievement of asymptotic accuracy using suitable iso-
geometric elements.
The paper is structured as follows: Section 2 outlines the 2D refined shell
theory and its variational principle. Section 3 details the rescaled, locking-
free variational formulation. Section 4 describes the finite element imple-
mentation. Section 5 presents numerical examples, including semi-cylindrical
shells under internal pressure and a composite structure, and Section 6 con-
cludes the paper.

2. 2D refined shell theory: Kinematics and variational principle

Consider a smooth, two-dimensional surface S within three-dimensional
Euclidean space, bounded by a continuous closed curve ∂S. We define a
shell volume V by constructing line segments of length h, orthogonal to S at
each point, with their midpoints on the surface. This definition is valid for
sufficiently small shell thickness h, ensuring no segment intersections. Here,
S represents the shell’s middle surface, and h its thickness (Fig. 1).
Mathematically, S is represented by the vector equation:

z = r(x1 , x2 ), (1)

or, in component form:

z i = ri (xα ), i = 1, 2, 3; α = 1, 2, (2)

3
 h

Figure 1: Schematic diagram of a shell segment

where z denotes the position vector in a 3D Cartesian system, and r(x1 , x2 )

is a smooth vector function. Curvilinear coordinates x1 and x2 are chosen
on S, with units of length.
The shell’s geometric characteristics are described by the first and second
fundamental forms:
aαβ = tα · tβ ,
(3)
bαβ = −tα · n,β = tα,β · n,
where tα = r,α are tangent vectors, and n is the unit normal vector to S,
given by:
t1 × t2
n= . (4)
|t1 ||t2 |
We utilize Latin indices (1 to 3) for Cartesian coordinates and Greek indices
(1 to 2) for surface coordinates. A comma preceding an index signifies partial
differentiation.
Within 2D linear kinematics, shell deformation is defined by the displace-
ment vector u(x1 , x2 ) = ui ei , with ei as Cartesian basis vectors. In the basis
{tα , n}, displacement components are:
uα = tα · u, u = n · u. (5)
The 2D refined shell theory (2D-RST), considering transverse shear, in-
troduces two additional degrees of freedom, φα (x1 , x2 ), representing shear-
induced rotation angles. The shell’s deformation is characterized by: (i)
extension measures:
i
γαβ = r,(α ui,β) = u(α;β) − bαβ u, (6)

4
(ii) bending measures:

ραβ = (ni ui,(α );β) + bλ(α ϖβ)λ − φ(α;β) = u;αβ + (uλ bλ(α );β) + bλ(α ϖβ)λ − φ(α;β) , (7)

where
1 i i 1
ϖαβ = (r,β ui,α − r,α ui,β ) = (uβ,α − uα,β ), (8)
2 2
and (iii) rotation angles φα . When φα = 0, ραβ reduces to classical shell
theory bending measures [3, 4]. Parentheses around indices denote sym-
metrization, and repeated indices imply summation. The metric tensor aαβ
and its inverse aαβ are used to lower or raise surface vector and tensor com-
ponent indices. Covariant derivatives, indicated by a semicolon in indices,
account for basis vector variations in surface tensors. For instance:
∇β uα ≡ uα;β = uα,β − Γλαβ uλ ,
(9)
∇β uα ≡ uα;β = uα,β + Γαβλ uλ ,

where Γλαβ are Christoffel’s symbols:

1
Γλαβ = aλµ (aµα,β + aµβ,α − aαβ,µ ). (10)
2
Note the difference in coordinate expressions for covariant derivatives of co-
and contravariant vector components in (9). These definitions extend to
tensor components [29].
The 2D-RST for linearly elastic, homogeneous shells is based on a varia-
tional principle [2, 6]. It states that the actual displacement field and rotation
angles minimize the two-dimensional average energy functional
Z
J[u, φα ] = Φ(γαβ , ραβ , φα ) dω − Aext (11)
S

among all admissible displacement fields and √ rotation angles that satisfy the
kinematic boundary conditions. Here, dω = adx1 dx2 is the area element,
and a = det aαβ . The two-dimensional energy density Φ(γαβ , ραβ , φα ) com-
prises three contributions, namely, (i) Φcl (γαβ , ραβ ): energy density of classi-
cal shell theory, (ii) Φgc (γαβ , ραβ ): geometric correction energy, (iii) Φsc (φα ):
shear correction energy. Thus:

Φ(γαβ , ραβ , φα ) = Φcl (γαβ , ραβ ) + Φgc (γαβ , ραβ ) + Φsc (φα ). (12)

5
These are explicitly given by:
h i µh3 h i
Φcl (γαβ , ραβ ) = µh σ(γαα )2 + γαβ γ αβ + σ(ραα )2 + ραβ ραβ ,
12
µh3 h αβ ′λ 3
Φgc (γαβ , ραβ ) = − ρ bα γβλ + σραβ bαβ γλλ + σρλλ bαβ γαβ
3 5 (13)
6  i
+ σ σ − 1 ρλλ Hγµµ ,
5
5
Φsc (φα ) = µhaαβ φα φβ ,
12
where λ and µ are Lamé’s constants, σ = λ+2µ λ ν
= 1−ν , and b′λ λ λ
α = bα − Hδα
is the deviator of the second quadratic form, with H = bαα /2 being the mean
curvature. The work done by external loads is expressed as
Z h
h σh 1 α β
Aext = (fi − hHgi )ui + gα (u,α + bλα uλ ) − (g + hf;α )γβ
S 2 2 6
1 1 α β 1 i
− σh2 (f + hg;α )ρβ − g α hφα dω , (14)
10 12 12
where
i
fi = τi |x3 =h/2 + τi |x3 =−h/2 , fα = fi r,α , f = f i ni ,
i
(15)
gi = τi |x3 =h/2 − τi |x3 =−h/2 , gα = gi r,α , g = gi ni ,

with τi representing the Cartesian components of the traction vector.

This 2D refined shell theory, by maintaining asymptotic accuracy up to
O(h/R) and O(h2 /l2 ), extends reliable analysis to shells of moderate thick-
ness, where R and l denote the characteristic radius of curvature and longi-
tudinal deformation length-scale, respectively [6, 2].

3. Rescaled formulation
The finite element (FE) implementation of the 2D variational problem
(Eq. (11)) presents challenges due to the differing orders of magnitude be-
tween bending and shear stiffnesses. Shear stiffness, significantly exceeding
bending stiffness by two orders relative to shell thickness h, can cause nu-
merical instability (shear locking) when multiplied by small rotation angles
as h decreases. A similar issue, membrane locking, arises from the contrast
between extension and bending stiffnesses, and their respective measures.

6
 To address these challenges, we introduce new unknown functions,

ψα = −ni ui,α + φα = −u,α − uλ bλα + φα , (16)

representing the total rotation angles of transverse fibers resulting from both
bending and shear, such that

ραβ = −ψ(α;β) + bλ(α ϖβ)λ ,

(17)
φα = u,α + uλ bλα + ψα .

Consequently, the functional in Eq. (11) is reformulated with these unknowns

as Z
J[u, ψα ] = Φ(γαβ , ραβ , u,α + uλ bλα + ψα ) dω − Aext , (18)
S
where function Φ of three arguments remains the same as in (12), while
ραβ should be taken from (17)1 . To avoid the membrane- and shear-locking
effects, we aim to make the problem independent of h, ensuring that exten-
sion, bending, and shear stiffnesses have comparable orders of magnitude.
We introduce rescaled coordinates:
zi xα
z̄ i = , x̄α = , (19)
h h
and express Eq. (2) as
z̄ i = r̄i (x̄α ), (20)
where r̄i (x̄α ) = ri (xα )/h. This rescaling renders the left-hand side of (19)
dimensionless. Importantly, the basis vectors and unit normal vector are
unaffected, preserving the metric tensor and its inverse. Specifically,
∂r̄ ∂(r/h)
r̄,ᾱ = α
= = r,α , (21)
∂ x̄ ∂(xα /h)
thus n̄ = n.
Displacements are unchanged, but to reflect their dependence on the new
base {r̄,ᾱ , n̄} and argument x̄α , they are denoted as:

ūᾱ (x̄α ) = uα (xα ), ū(x̄α ) = u(xα ). (22)

The second fundamental form of the surface transforms as

b̄ᾱβ̄ = r̄,ᾱβ̄ · n̄ = hr,αβ · n = hbαβ , (23)

7
resulting in scaled principal radii of curvature R̄1 = R1 /h and R̄2 = R2 /h.
Partial derivatives with respect to x̄α are related to those with respect to xα
by
∂ ∂
α
= h α , (.),ᾱ = h(.),α . (24)
∂ x̄ ∂x
From Eq. (10), the scaled Christoffel symbols are related to the original ones
by
1 1
Γ̄λ̄ᾱβ̄ = āλ̄µ̄ (āµ̄ᾱ,β̄ +āµ̄β̄,ᾱ −āᾱβ̄,µ̄ ) = haλµ (aµα,β +aµβ,α −aαβ,µ ) = hΓλαβ . (25)
2 2
Therefore, covariant derivatives with respect to x̄α are also related to those
with respect to xα by
¯ ᾱ = h∇α ,
∇ (.);ᾱ = h(.);α . (26)
For rotation angles, we introduce
ψ̄ᾱ = hψα , (27)
with dimensions of length, representing longitudinal displacements of the
shell’s positive face [6].
Applying the scaling rules (19), (22), and (27), we obtain the following
relationships between original and rescaled quantities:
1 1
γαβ = γ̄ , ραβ = ρ̄ ,
h ᾱβ̄ h2 ᾱβ̄
1 (28)
u,α + uλ bλα + ψα = (ū,ᾱ + ūλ̄ b̄λ̄ᾱ + ψ̄ᾱ ),
√ h √
dω = a dx dx = h2 ā dx̄1 dx̄2 = h2 dω̄ ,
1 2

where
γ̄ᾱβ̄ = ū(ᾱ;β̄) − b̄ᾱβ̄ ū,
(29)
ρ̄ᾱβ̄ = −ψ̄(ᾱ;β̄) + b̄λ̄(ᾱ ϖ̄β̄)λ̄ ,
and
1
ϖ̄ᾱβ̄ = (ūβ̄,ᾱ − ūᾱ,β̄ ). (30)
2
Substituting Eqs. (28) into the energy functional (18) gives a form with
rescaled quantities:
Z
J[ū, ψ̄ᾱ ] = µh Φ̄(γ̄ᾱβ̄ , ρ̄ᾱβ̄ , ū,ᾱ + ūλ̄ b̄λ̄ᾱ + ψ̄ᾱ ) dω̄ − Aext . (31)
S̄

8
Here S̄ = {(x̄1 , x̄2 ) | (x1 , x2 ) ∈ S} denotes the rescaled 2D domain, while
the rescaled energy density, Φ̄(γ̄ᾱβ̄ , ρ̄ᾱβ̄ , ū,ᾱ + ūλ̄ b̄λ̄ᾱ + ψ̄ᾱ ), is the sum of three
contributions
1 h ᾱ 2 i
Φ̄cl (γ̄ᾱβ̄ , ρ̄ᾱβ̄ ) = σ(γ̄ᾱᾱ )2 + γ̄ᾱβ̄ γ̄ ᾱβ̄ + σ(ρ̄ᾱ ) + ρ̄ᾱβ̄ ρ̄ᾱβ̄ ,
12
1 h ᾱβ̄ ′λ̄ 3
Φ̄gc (γ̄ᾱβ̄ , ρ̄ᾱβ̄ ) = − ρ̄ b̄ᾱ γ̄β̄ λ̄ + σ ρ̄ᾱβ̄ b̄ᾱβ̄ γ̄λ̄λ̄ + σ ρ̄λ̄λ̄ b̄ᾱβ̄ γ̄ᾱβ̄
3 5 (32)
6  i
λ̄ µ̄
+ σ σ − 1 ρ̄λ̄ H̄ γ̄µ̄ ,
5
5
Φ̄sc (ū,ᾱ + ūλ̄ b̄λ̄ᾱ + ψ̄ᾱ ) = āᾱβ̄ (ū,ᾱ + ūλ̄ b̄λ̄ᾱ + ψ̄ᾱ )(ū,β̄ + ūµ̄ b̄µ̄β̄ + ψ̄β̄ )
12
similar to (12). The work of external forces becomes
Z h
1 σ 1
Aext = h 2
(fi − H̄gi )ūi + g ᾱ (ū,ᾱ + ūλ̄ b̄λ̄ᾱ ) − (g + f;ᾱᾱ )γ̄β̄β̄
S 2 2 6
1 1 ᾱ β̄ 1 ᾱ λ̄
i
+ σ(f + g;ᾱ )ψ̄;β̄ − g (ū,ᾱ + ūλ̄ b̄ᾱ + ψ̄ᾱ ) dω̄ (33)
10 12 12
Since scaling a functional by a constant does not alter its minimizer, we
further simplify (31) by dividing by µh. The minimization problem becomes
Z
¯ ψ̄ᾱ ] =
J[ū, Φ̄(γ̄ᾱβ̄ , ρ̄ᾱβ̄ , ū,ᾱ + ūλ̄ b̄λ̄ᾱ + ψ̄ᾱ ) dω̄ − Āext → min , (34)
S̄ ūᾱ ,ū,ψ̄ᾱ

where
Z h
1 σ 1
Āext = (f¯i − H̄ ḡi )ūi + ḡ ᾱ (ū,ᾱ + ūλ̄ b̄λ̄ᾱ ) − (ḡ + f¯;ᾱᾱ )γ̄β̄β̄
S 2 2 6
1 ¯ 1 ᾱ β̄ 1 ᾱ λ̄
i
+ σ(f + ḡ;ᾱ )ψ̄;β̄ − ḡ (ū,ᾱ + ūλ̄ b̄ᾱ + ψ̄ᾱ ) dω̄ , (35)
10 12 12
and
hfi hgi
f¯i = , ḡi = . (36)
µ µ
Consequently, both f¯i and ḡi scale with shell thickness h and strain ε =
max{fi /µ, gi /µ}. Since all three stiffnesses (extension, bending, and shear)
in the rescaled functional (34) are O(1), the minimizer is O(f¯i ,ḡi ) times a
function depending on the characteristic size of S̄. Returning to original
functions, extension measures γαβ and rotation angles ψα are O(ε) (with

9
φα smaller [29, 2]), while bending measures ραβ are O(ε/h). The rescaled
problem (34) avoids membrane and shear locking.
To solve problem (34), boundary conditions must be specified. Typi-
cally, we assume a portion of the shell’s edge, ∂¯k , is clamped, requiring the
admissible functions to satisfy the following kinematic conditions:

ūᾱ = 0, ū = 0, ψ̄ᾱ = 0 at ∂¯k . (37)

For a simply supported remaining edge ∂¯s , kinematical conditions ūᾱ = 0,

ū = 0, and ψ̄ᾱ τ̄ ᾱ = 0 are enforced, while ψ̄ᾱ ν̄ ᾱ is allowed to vary freely,
where τ̄ ᾱ and ν̄ ᾱ are the components of the surface tangential and normal
unit vectors on ∂¯s , respectively. On a free edge ∂¯s , no constraints are imposed
on ūᾱ , ū and ψ̄ᾱ . Beyond these typical cases, Section 5 explores additional
boundary conditions to validate our FE implementation.
After obtaining the solution, the theory’s asymptotic accuracy is assessed
by comparing the true average displacement with 3D elasticity results. Ac-
cording to [6], ūᾱ represent the true average tangential displacements. The
true average normal displacement is obtained by correcting the normal dis-
placement:
h2 σ αβ σ
ǔ = u + a ραβ = ū + āᾱβ̄ (−ψ̄(ᾱ;β̄) + b̄λ̄(ᾱ ϖ̄β̄)λ̄ ). (38)
60 60
Since φα also depends on first derivatives of u, we seek solutions where ūᾱ , ū,
and ψ̄ᾱ are C 1 -functions for accurate computations and comparisons. Asymp-
totic accuracy is confirmed by agreement up to O(h/R) and O(h2 /l2 ) between
functions ūᾱ , ū, and ψ̄ᾱ found by the 2D refined shell theory with
i
r,α ⟨wi (xα , x3 )⟩, ni ⟨wi (xα , x3 )⟩, i
and r,α ⟨wi (xα , x3 )x3 ⟩/(h2 /12), (39)
R h/2
where ⟨.⟩ ≡ h1 −h/2 . dx3 denotes the averaging over the shell’s thickness
and wi (xα , x3 ) are the displacements computed by the 3D exact theory of
elasticity.

4. Finite element implementation

4.1. Weak and strong formulations
In this Section and the theoretical portion of the subsequent Section, we
will use rescaled coordinates and quantities exclusively. For brevity, overbars
will be omitted.

10
 Taking the first variation of the functional (Eq. (34)), a necessary con-
dition for its minimizer is that the virtual work of internal forces equals the
virtual work of external forces:
Z h i
nαβ δγ αβ + mαβ δραβ + q α (δu,α + bλα δuλ + δψ α ) dω = δAext , (40)
S

where the membrane forces nαβ , bending moments mαβ , and shear forces q α
are dual to γαβ , ραβ and φα , respectively. These are defined as:

∂Φ 1h ′β)

nαβ = = 2(σγλλ aαβ + γ αβ ) − ρ(αλ bλ + σaαβ bµν ρµν
∂γαβ 3
6   3 i
+ σ − 1 Hρλλ + σρλλ bαβ ,
5 5
∂Φ 1 1h ′β)
3
mαβ = = (σρλλ aαβ + ραβ ) − γ (αλ bλ + σaαβ bµν γ µν (41)
∂ραβ 6 3 5
6   i
+ σ − 1 Hγλλ + σγλλ bαβ ,
5
∂Φ 5
qα = = aαβ (u,β + bλβ uλ + ψβ ).
∂φα 6
The virtual work of external forces is:
Z h
1 σ 1 λ αβ
δAext = (f α −Hg α ) δuα +(f −Hg) δu+ g α (δu,α +bλα δuλ )− (g+ f;λ )a δγ αβ
S 2 2 6
σ 1 λ αβ 1 α λ
i
+ (f + g;λ )a δψ α;β − g (δu,α + bα δuλ + δψ α ) dω . (42)
10 12 12
Introducing the notations
Z Z
m αβ
n 1h ′β)

δW = n δγ αβ dω = dω 2(σγλλ aαβ +γ αβ )− ρ(αλ bλ +σaαβ bµν ρµν
S S 3
6   3 io
+ σ − 1 Hρλλ + σρλλ bαβ (δu(α;β) − bαβ δu) dω , (43)
5 5

Z Z n1
b αβ 1h ′β)
3
δW = m δραβ dω = dω (σρλλ aαβ +ραβ )− γ (αλ bλ +σaαβ bµν γ µν
S S 6 3 5
6   io
+ σ − 1 Hγλλ + σγλλ bαβ (− δψ (α;β) + bλ(α δϖβ)λ ) dω , (44)
5

11
and
Z Z
s α
δW = q δφα dω = aαβ (u,β +bλβ uλ +ψβ )(δu,α +bµα δuµ +δψ α ) dω , (45)
S S

for the virtual work of membrane forces, bending moments, and shear forces,
respectively, Eq. (40) becomes:

δW m + δW b + δW s = δAext . (46)

Let K = {(v, vα , χα ) | (v, vα , χα )|∂k = 0} be the space of kinematically ad-

missible functions. The weak formulation of the refined shell theory states
that for given f, f α , g, g α find (u, uα , ψα ) ∈ K such that Eq. (46) is satisfied
for all (δu , δuα , δψ α ) ∈ K. Looking at the integrals (43)-(45) we see that,
for this weak formulation to make sense, the admissible functions should be-
long at least to the Sobolev’s space of square integrable functions with the
square integrable first derivatives, H 1 (S). However, since the desired asymp-
totic accuracy of FSDT may require higher smoothness of u, uα , and ψα , the
continuity assumption in K still remains unspecified.
For completeness, we also present the strong formulation of the problem.
Assuming u, uα , and ψα are twice differentiable, we integrate Eq. (40) by
parts, yielding:
Z h
λ[α β]
−nαβ
;β δuα − bαβ n
αβ
δu + mαβ
;β δψ α − (m
α
bλ );β δuα − q;α δu + bαβ q β δuα
S Z h
i i
α β]
+ q δψ α ) dω + nαβ νβ δuα − mαβ νβ δψ α + mλ[α bλ νβ δuα + q α να δu ds
∂s
Z h
5 α 1 σ 1 λ
= (f α −Hg α ) δuα +(f −Hg) δu− g;α δu+ g β bαβ δuα + (g+ f;λ );β aαβ δuα
S 12 2 2 6
1 λ σ 1 λ 1 1 i
+ σH(g + f;λ ) δu − (f + g;λ );β aαβ δψ α − g β bαβ δuα − g α δψ α dω
Z h 6 10 12 12 12
5 α σ 1 λ σ 1 i
+ g να δu − (g + f;λ )νβ aαβ δuα + (f + g;λ λ
)νβ aαβ δψ α ds .
∂s 12 2 6 10 12
(47)

Here, square brackets enclosing indices denote anti-symmetrization: t[αβ] =

1 αβ
2
(t − tβα ). Vector να represents the outward surface normal to ∂S, and ds
is the length element. We assume the remaining portion of the shell edge,
denoted by ∂s , is free. Due to the arbitrariness of the variations δuα , δψ α , and

12
δu in S and on ∂s , we obtain the following second-order partial differential
equations from Eq. (47):
1 β α σ 1 λ 1
− tαβ α β α α
;β + bβ q = f − Hg + g bβ + (g + f;λ );β aαβ − g β bαβ
2 2 6 12
σ 1 1
mαβ α
;β + q = − (f + g;λ λ
);β aαβ − g α , (48)
10 12 12
α 5 1 λ
− q;α − bαβ nαβ = f − Hg − g;α α
+ σH(g + f;λ ),
12 6
β]
where tαβ = nαβ + mλ[α bλ . These equations are subject to the kinematic
boundary conditions (Eq. (37)) on ∂k and the following natural boundary
conditions on ∂s :
σ 1 λ
tαβ νβ = − (g + f;λ )νβ aαβ ,
2 6
σ 1 λ
− mαβ νβ = (f + g;λ )νβ aαβ , (49)
10 12
5
q α να = g α να
12
The derivation of boundary conditions for a simply supported edge ∂s is
similar. Equations (48), (37), and (49) constitute the strong formulation of
the problem.

4.2. Discretization
The weak formulation (Eqs. (42), (44), (43), and (45)) involves covariant
derivatives of various co- and contravariant vector and tensor components,
complicating direct computation of residual forces and stiffness matrices.
Therefore, unlike the first-order shear deformation theory for plates presented
in [43], this refined shell theory employs a different approach. Leveraging
the inherent smoothness of the weak formulation, Automatic Differentiation
(AD) is used for symbolic representation. Specifically, the Trilinos/Sacado
package provides the necessary AD capabilities.
The primal variables, u, uα , ψα , and their corresponding variations δu,
δuα , δψα , are assigned degrees of freedom (d.o.f.) with indices iu , iuα , iψα , iδu ,
iδuα , and iδψα , respectively. Defining δJ = δW m + δW b + δW s − δAext as the
first variation of the energy functional, the residual force and stiffness matrix
corresponding to primal d.o.f. δa and b (a, b = u, uα , ψα ) are computed as:


∂ δW m + δW b + δW s − δAext
rδa = = δJ .dx (iδa ) (50)
∂δa
13


∂rδa ∂ δW m + δW b + δW s − δAext
Kδa,b = = = δJ .dx (iδa ) .dx (ib ) (51)
∂b ∂δa∂b
We note that, to compute the double derivatives in Eq. (51), the symbolic
variables δa and b must be of type Sacado::Fad::DFad<Sacado::Fad::DFad<double>>.

4.3. Isogeometric analysis

As mentioned at the end of Section 3, accurate determination of the
true average normal displacement ǔ and rotation angle φα necessitates C 1
continuity. Isogeometric analysis (IGA) is well-suited for shell geometry dis-
cretization for two primary reasons: (i) Higher-order continuity: IGA readily
provides the required C 1 continuity through the use of non-uniform rational
B-splines (NURBS) shape functions, (ii) Geometric flexibility: Complex ge-
ometries can be accurately represented by multiple surface patches connected
at interfaces, with continuity preserved using techniques like the bending strip
method [44]. Furthermore, the order of the NURBS basis functions can be
conveniently increased via hpk-refinement.
NURBS shape functions used in IGA read:
m X
n
X N p (ξ1 ) Njq (ξ2 )
S (ξ1 , ξ2 ) = Pm Pni p q Pij . (52)
i=1 j=1 k=1 l=1 wkl Nk (ξ1 ) Nl (ξ2 )

Here, Nip and Njq are univariate B-spline basis functions of order p and q,
respectively, computed using the recursive Cox-de-Boor formula:
(
ξ − ξi ξi+p+1 − ξi 1 ξi ≤ ξ ≤ ξi+1
Nip (ξ) = Nip−1 (ξ)+ p−1
Ni+1 (ξ), Ni0 (ξ) = .
ξi+p − ξi ξi+p+1 − ξi+1 0 otherwise
(53)
{Pij }0≤i≤m,0≤j≤n represents the control point grid, and {wij }0≤i≤m,0≤j≤n
are the corresponding control weights. The Cox-de-Boor formula (Eq. (53))
requires a global knot vector. To maintain the local nature of finite elements
and enable parallel computation during assembly, Bézier decomposition [45]
is employed.

4.4. Structure of the finite element code

The proposed FSDT shell element is implemented within the Plate-
AndShellApplication extension of the Kratos Multiphysics framework ker-
nel [46]. The IsogeometricPlateAndShellApplication extension enables

14
 Kratos
(Variables, Database, Python interface, . . . )

IsogeometricApplication
PlateAndShellApplication
(IGA modeller, bending strip
(FSDT shell element)
patch, Bézier finite element)

IsogeometricPlateAndShellApplication
(Interlinking extension)

Figure 2: Structure of the computational code.

analysis using IGA. NURBS multipatch structures and Bézier elements are
supported in the IsogeometricApplication. The relationships between
these Kratos extensions are shown in Fig. 2.

5. Numerical examples
5.1. Semi-cylindrical shell under internal pressure with freely sliding side
edges
We analyze a semi-cylindrical shell occupying the region defined by
0 ≤ x ≤ L, 0 ≤ θ ≤ π, R − h/2 ≤ r ≤ R + h/2 (54)
in cylindrical coordinates {x, θ, r} (see Fig. 3). The shell is subjected to
internal pressure p applied to its inner surface at r = R − h/2. For the
2D shell theory, we introduce rescaled surface coordinates x̄1 = x/h varying
from 0 to L̄ = L/h and x̄2 = Rθ/h varying from 0 to W̄ = π R̄ = πR/h,
where W = πR is the semi-circular circumference. For brevity, we drop the
overbars on rescaled coordinates and quantities in the following theoretical
development.
The middle surface in rescaled coordinates is described by:
x2 x2
z = −x1 e1 + R cos e2 + R sin e3 . (55)
R R
15
 freely sliding edge

freely sliding edge L

R h

Figure 3: Schematic diagram of the semi-cylindrical shell

The tangent vectors to the coordinate lines on the middle surface are:
x2 x2
t1 = z,1 = −e1 , t2 = z,2 = − sin e2 + cos e3 . (56)
R R
The unit normal vector is:
t1 × t2 x2 x2
n= = cos e2 + sin e3 . (57)
|t1 × t2 | R R
In the 2D curvilinear coordinate system {x1 , x2 }, the components of the 2D
metric tensor are δαβ , the Kronecker delta. Consequently, the Christoffel
symbols vanish, and covariant derivatives coincide with partial derivatives.
Raising or lowering indices does not change tensor and vector components,
so we represent them in covariant form with lower indices. The non-zero
component of the second fundamental form is:
1
b11 = 0, b12 = 0, b22 = − . (58)
R
The mean curvature and b′αβ are:
1 1 1
H=− , b′11 = , b′12 = 0, b′22 = − . (59)
2R 2R 2R
16
 Using Eq. (28), the extension and bending measures, and the rotation
angles due to shear deformation are:
1 u
γ11 = u1,1 , γ12 = γ21 = (u1,2 + u2,1 ), γ22 = u2,2 + ,
2 R
ρ11 = −ψ1,1 , ρ22 = −ψ2,2 ,
1 1 (60)
ρ12 = ρ21 = − (ψ1,2 + ψ2,1 ) + (u1,2 − u2,1 ),
2 4R
u2
φ1 = u,1 + ψ1 , φ2 = u,2 − + ψ2 .
R
The energy density contributions (with overbars dropped) are:
u u 1
Φcl = σ(u1,1 + u2,2 + )2 + (u1,1 )2 + (u2,2 + )2 + (u1,2 + u2,1 )2
R R 2
1h 2 2 2
+ σ(ψ1,1 + ψ2,2 ) + (ψ1,1 ) + (ψ2,2 ) (61)
12
1  1 2 i
+ −ψ1,2 − ψ2,1 + (u1,2 − u2,1 ) ,
2 2R
1 h 6  6
Φgc = σ( σ − 1) − 1 ρ11 γ11 + (1 + σ)(1 + σ)ρ22 γ22
6R 5 5 (62)
6 1 6 i
+ σ(1 + σ)ρ22 γ11 + σ( + σ)ρ11 γ22 ,
5 5 5
5h u2 i
Φsc = (u,1 + ψ1 )2 + (u,2 − + ψ2 )2 . (63)
12 R
With pressure acting on the inner face,
f = p, g = −p, fα = gα = 0, (64)
and the external work is:
Z h
1 σ u σ i
Aext = p(1 − )u + p(u1,1 + u2,2 + ) + p(ψ1,1 + ψ2,2 ) dω . (65)
S 2R 2 R 10
We seek analytical solutions as benchmark solutions to validate our finite
element implementation of 2D problems. For this purpose, we assume fric-
tionless sliding of the side edges between rigid planes at x1 = 0 and x1 = L.
This implies the following boundary conditions:
u1 = ψ1 = 0, u2 , u, ψ2 -arbitrary at x1 = 0, L. (66)
For the bottom edges of the shell, located at x2 = 0, W , we consider three
boundary condition cases:

17
 1. Freely sliding edges:
u1 = u2 = 0, u -arbitrary, ψ1 = ψ2 = 0 at x2 = 0, W . (67)
2. Clamped edges:
u1 = u2 = u = ψ1 = ψ2 = 0 at x2 = 0, W . (68)
3. Simply supported edges:
u1 = u2 = u = ψ1 = 0, ψ2 -arbitrary at x2 = 0, W . (69)
These boundary conditions lead to a state of plane strain, where:
u1 ≡ ψ1 ≡ 0, u2 , u, ψ2 are function of x2 only. (70)
Consequently, Eq. (60) simplifies to:
1
γ11 = 0, γ12 = 0, γ22 ≡ γ = u2,2 + u,
R
ρ11 = 0, ρ12 = 0, ρ22 ≡ ρ = −ψ2,2 , (71)
1
φ1 = 0, φ2 ≡ φ = u,2 − u2 + ψ2 ,
R
where the comma before the index 2 denotes the ordinary derivative with
respect to x2 . The minimization problem (Eq. (34)) reduces to a 1-D problem
(arc-like model, 1D-RST): Minimize
Z Wh
1 1+σ 1 6 1 1
J= (1+σ)(u2,2 + u)2 + (ψ2,2 )2 − (1+σ)(1+ σ) (u2,2 + u)ψ2,2
0 R 12 6 5 R R
5 1 1 σ 1 σ i
+ (u,2 − u2 + ψ2 )2 − p(1 − )u − p(u2,2 + u) − pψ2,2 dx2
12 R 2R 2 R 10
(72)
with respect to u2 , u, ψ2 subject to the kinematic boundary conditions. The
integration over x1 yields a constant factor L, which is omitted in Eq. (72).
The Euler-Lagrange equations, obtained from the vanishing first variation
of functional (72), are:
1
−n,2 − q = 0,
R
m,2 + q = 0, (73)
1 1−σ
−q,2 + n = p(1 − ),
R 2R
18
where
∂Φ 1 6 1
n= = 2(1 + σ)γ + (1 + σ)(1 + σ) ρ,
∂γ 6 5 R
∂Φ 1 1 6 1
m= = (1 + σ)ρ + (1 + σ)(1 + σ) γ, (74)
∂ρ 6 6 5 R
∂Φ 5
q= = φ.
∂φ 6
Substituting Eq. (74) with γ, ρ, φ from Eq. (71) into Eq. (73) yields three
second-order ordinary differential equations for u2 , ψ2 , u.
For comparison, we also analyze the shell using Sanders-Koiter classical
shell theory (CST). Assuming plane strain, the problem reduces to minimiz-
ing the 1-D functional (1D-CST)
Z Wh
1 1+σ 1 i
Jcl = (1 + σ)(u2,2 + u)2 + (u,22 − u2,2 )2 − pu dx2 (75)
0 R 12 R
with respect to u2 and u subject to the kinematic boundary conditions. The
Euler-Lagrange equations are:
1
−n,2 − m,2 = 0,
R (76)
1
m,22 + n = p,
R
where
1
n = 2(1 + σ)γ, γ = u2,2 + u,
R (77)
1 1
m = (1 + σ)ρ, ρ = u,22 − u2,2 .
6 R
For Case 1 (Eq. (67)), the boundary conditions derived from the vanishing
first variation of J and the arbitrariness of u at at x2 = 0, W are:

u2 = ψ2 = 0, q = 0 at x2 = 0, W . (78)

The solution is
pR2 1−σ
u= (1 − ), u2 = ψ2 = 0, (79)
2(1 + σ) 2R

19
 10−5
2

L2 Error
10−7 1

2D-RST
10−9
101 102 103 104 105 106
Degrees of freedom

Figure 4: Convergence of 2D-RST displacement in the first case against the analytical
solution of 1D-RST in terms of L2 error.

yielding
pR 1−σ
γ= (1 − ), ρ = φ = 0. (80)
2(1 + σ) 2R
This solution satisfies Eqs. (73), (74), and the boundary conditions (Eq. (69)).
This problem has a known 3D elasticity solution. The radial displacement is
[47]:
pR2 1 2h 1 2 1 i
wr = (1 − ) (1 − 2ν)(1 + ζ) + (1 + ) , (81)
4 2R 2R 1 + ζ
where ζ = ξ/R and ζ ∈ (−1/(2R), 1/(2R)). The average displacement from
Eq. (81) agrees with Eq. (79) up to h/R. Figure 4 shows the L2 error conver-
gence of the 2D-RST displacement compared to Eq. (79). Classical shell the-
ory (CST) [3, 4], while accurate for stress resultants and bending moments,
exhibits large displacement errors. The 1D-CST solution (Eqs. (75)-(77)) for
mean displacements is
1 π 
u = − pR2 (1 − ν) sin(x2 /R) − 2 ,
2 2 (82)
1 2 π π 
u2 = pR (1 − ν) − cos(x2 /R) − x2 /R ,
2 2 2
which results in a 100% error compared to Eqs. (79) and (81). This discrep-
ancy highlights the limitations of classical shell theory in accurately predict-
ing displacements, as explained in [5].
For Case 2 with clamped bottom edges, an analytical solution of 1D-RST
can be found in terms of exponential functions involving the roots of a cubic

20
Figure 5: Schematic diagram of the half-ring of thickness 1 and radius R̄ under internal
pressure p̄.

polynomial equation. However, this solution is quite cumbersome. Therefore,

we opt for a numerical approach. We first rewrite Eq. (71) as:
1
u2,2 = γ − u,
R
ψ2,2 = −ρ, (83)
1
u,2 = φ + u2 − ψ2 ,
R
Next, we eliminate ρ from Eqs. (74)1 and (74)2 to express γ in terms of n
and m. Differentiating and using Eqs. (73)1 and (73)2 , we obtain:
σ
γ,2 = 1 φ. (84)
2(1 + σ)R[1 − 12 (1 + 65 σ)2 R12 ]

Similarly, eliminating γ from Eqs. (74)1 and (74)2 , expressing ρ in terms of

n and m, differentiating, and using Eqs. (73)1 and (73)2 , we get:
1
5(1 − 12 (1 + 56 σ) R12 )
ρ,2 = − 1 φ. (85)
(1 + σ)[1 − 12 (1 + 65 σ)2 R12 ]

Finally, from Eqs. (73)3 and (74)3 , we have:

6 1 1 6 1 1−σ
φ,2 = [2(1 + σ) γ + (1 + σ)(1 + σ) 2 ρ − p(1 − )]. (86)
5 R 6 5 R 2R

21
 10−5
60
2

L2 Error
40 1D-CST
1D-RST 10 −6 1
u
2D-RST
20 3D-ET
2D-RST
10−7
0
0 5 10 15 20 25 30 102 103 104
x̄2 Degrees of freedom

Figure 6: Normal displacement ū of a semi-cylindrical shell under external pressure with

freely sliding side edges and clamped bottom edges (second case) as function of the rescaled
circumferential coordinate x̄2 = π R̄θ (R̄ = 10).

4
1D-CST
2 1D-RST
ψ 0 2D-RST
3D-ET
−2
−4
0 5 10 15 20 25 30
x̄2

Figure 7: Rescaled rotation angle ψ̄2̄ of a semi-cylindrical shell under external pressure
with freely sliding side edges and clamped bottom edges (second case) as function of the
rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 10).

Introducing v = u,2 and ϑ = ρ,2 , the 1D-CST system (Eqs. (76) and (77))
becomes:
u,2 = v,
1
u2,2 = γ − u,
R
1
γ,2 = ϑ,
12R (87)
1 1
v,2 = ρ + (γ − u),
R R
ρ,2 = ϑ,
6 12
ϑ,2 = p − γ.
1+σ R
Equations (83) and (84)-(86) form a system of six first-order ODEs for

22
u2 , ψ2 , u, γ, ρ, φ, subject to the boundary conditions:

u2 = ψ2 = u = 0 at x2 = 0, W . (88)

For 1D-CST, we have six first-order ODEs (Eq. (87)) for u, u2 , γ, v, ρ, ϑ, sub-
ject to:
u = u2 = v = 0 at x2 = 0, W . (89)
These two-point boundary-value problems are solved using Matlab’s bvp4c
function.
To demonstrate the asymptotic accuracy of our FE-implementation of 2D-
RST, we also compare its numerical solution with that of 3D elasticity theory
(3D-ET). In this case, under plane strain conditions, the rescaled problem
reduces to the two-dimensional analysis of a half-ring having thickness 1 (see
Fig. 5), which can be solved using 2D isogeometrical analysis (IGA) with
solid elements. After solving this 2D problem, the mean radial displacement
and rotation angle over the shell’s thickness are evaluated in accordance with
Eq. (39).
Both the shell analysis and the elastic solid analysis employ cubic-order
NURBS discretization with sufficiently fine meshes. Using overbars to denote
rescaled quantities, we present the results of numerical simulation below. In
all numerical simulations, we set ν = 0.3 and p̄ = 1, so that the solution is
normalized. Due to linearity, the solution for any real p̄ can be obtained by
multiplying the normalized solution by p̄.
Fig. 6 shows, on the right, the convergence rate of the 2D-RST solution
to that of the 1D-RST for the thin shell with R̄ = 10, demonstrating that our
FE-implementation is free from membrane and shear locking. It also plots, on
the left, the normal mean displacement ū versus x̄2 computed using all four
theories (1D-CST, 1D-RST, 2D-RST, 3D-ET). The 2D-RST solution (red
line) matches that of 1D-RST (black dotted-dashed line) perfectly. Since
the shell is thin, this 2D-RST solution agrees quite well with that of 3D-
ET (black plus points). The 1D-CST solution deviates noticeably from the
others. Fig. 7 shows the plot of the rescaled rotation angle ψ̄2̄ versus x̄2
computed using all four theories. For classical shell theory, the rotation
angle ψ̄2̄ is due to bending only and is computed as ψ̄2̄ = −ū,2̄2̄ + ū2̄,2̄ /R̄. We
observe a perfect match between the 2D-RST and 1D-RST solutions, and an
almost perfect match between them and the 3D-ET solution. Here too, the
1D-CST solution deviates noticeably from the others.

23
 6

1D-CST
4 1D-RST
2D-RST
u
3D-ET
2

0
0 1 2 3 4 5 6 7 8 9
x̄2

1 1D-CST
1D-RST
ψ 0 2D-RST
3D-ET

−1

0 1 2 3 4 5 6 7 8 9
x̄2

Figure 8: Normal displacement ū and rescaled rotation angle ψ̄2̄ of a semi-cylindrical shell
under external pressure with freely sliding side edges and clamped bottom edges (second
case) as function of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 3).

We now check whether 2D-RST is applicable to moderately thick shells.

For this purpose, we set R̄ = 3, while keeping all other parameters unchanged.
The plots of the mean normal displacement ū and rescaled rotation angle ψ̄2̄
(−ū,2̄2̄ + ū2̄,2̄ /R̄ for CST) as functions of x̄2 are shown in Fig. 8 for ν = 0.3,
R̄ = 3, and p̄ = 1. The curves for 1D-RST and 2D-RST, which match
each other perfectly, agree well with that of 3D-ET, showing that 2D-RST is
applicable also to moderately thick shells. The classical shell theory shows a
large deviation from 3D-ET, indicating that CST can no longer be used for
moderately thick shells.
For Case 3 with simply supported bottom edges, the equations remain
the same, but the boundary conditions become:
σ
u2 = u = 0, m=− p at x2 = 0, W . (90)
10
σ
The conditions m = − 10 p at x2 = 0, W arise from the vanishing first variation
of functional (72) and the arbitrariness of ψ2 at x2 = 0, W . For 1D-CST, the

24
 60

10−5 2
40

L2 Error
1D-CST
1D-RST 1
u
2D-RST −6
20 10
3D-ET
2D-RST
0
0 5 10 15 20 25 30 102 103 104
x̄2 Degrees of freedom

Figure 9: Normal displacement ū of a semi-cylindrical shell under external pressure with

freely sliding side edges and simply-supported bottom edges (third case) as function of the
rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 10).

5
1D-CST
1D-RST
ψ 0 2D-RST
3D-ET

−5

0 5 10 15 20 25 30
x̄2

Figure 10: Rescaled rotation angle ψ̄2̄ of a semi-cylindrical shell under external pressure
with freely sliding side edges and simply-supported bottom edges (third case) as function
of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 10).

boundary conditions are:

u = u2 = ϑ = 0 at x2 = 0, W . (91)

These two-point boundary-value problems are solved using Matlab’s bvp4c

function. Concerning the boundary conditions within 3D-ET: Since exact
boundary conditions within 3D-ET for the simply supported edge are ab-
sent, we propose boundary conditions that best mimic those of the 2D shell
theory. As such, we impose the following constraints: under the plane strain
conditions, the mean displacements must vanish:

⟨w2 ⟩ = 0, ⟨w⟩ = 0. (92)

Another condition, namely that the rotation ψ2 is not constrained, cannot be

σ
directly imposed within 3D-ET. However, its consequence, m = − 10 p, can

25
 6

1D-CST
4 1D-RST
2D-RST
u
3D-ET
2

0
0 1 2 3 4 5 6 7 8 9
x̄2
2

1 1D-CST
1D-RST
ψ 0 2D-RST
3D-ET
−1

−2
0 1 2 3 4 5 6 7 8 9
x̄2

Figure 11: Normal displacement ū and rescaled rotation angle ψ̄2̄ of a semi-cylindrical shell
under external pressure with free side edges and simply-supported bottom edges (third
case) as function of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 3).

be imposed by specifying at the bottom edges the traction that gives such a
bending moment, namely:
6σ
σ22 (r̄, θ) = − p̄(r̄ − R̄) at θ = 0, π. (93)
5
After solving the plane strain problem for the half-ring, the mean radial
displacement and rotation angle over the shell’s thickness are evaluated in
accordance with Eq. (39), as in the previous case.
The mean normal displacement ū and rescaled rotation angle ψ̄2̄ (−ū,2̄2̄ +
ū2̄,2̄ /R̄ for CST) as functions of x̄2 are shown in Fig. 9 (left) and Fig. 10,
respectively, for ν = 0.3, R̄ = 10, and p = 1. In Fig. 9 (right) we show
also the convergence rate of the 2D-RST solution to that of the 1D-RST.
We observe again the perfect match between 2D-RST and 1D-RST, and the
almost perfect match between them and 3D-ET. The CST differs from them
noticeably.
If we set R̄ = 3 (moderately thick shell) while keeping all other parame-

26
 1D-CST
60
1D-RST
2D-RST(L=10R)
40 2D-RST(L=100R)
3D-ET
20

0
0 5 10 15 20 25 30

4 1D-CST
1D-RST
2 2D-RST (L=10R)
0 2D-RST (L=100R)
3D-ET
−2
−4
0 5 10 15 20 25 30

Figure 12: Normal displacement ū and rescaled rotation angle ψ̄2̄ of a semi-cylindrical
shell under external pressure with free side edges and clamped bottom edges (second case)
as function of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 10).

ters unchanged, we see that 2D-RST and 1D-RST still agree well with 3D-ET
(see Fig. 11, which plots ū and ψ̄2̄ computed according to all four theories).
However, CST deviates considerably from the others. This comparison con-
firms the applicability of 2D-RST for moderately thick shells, but also shows
that 2D-CST can no longer be used.

5.2. Semi-cylindrical shell under internal pressure with free side edges
This subsection investigates the behavior of a semi-cylindrical shell sub-
jected to internal pressure, focusing on the impact of free side edges and
varying bottom edge constraints. We analyze the normal displacement ū and
rotation angle ψ̄2̄ along the circumferential coordinate, leveraging reduced 1-
D models and comparisons with 2D-RST solution in the mid-cross-section
and plane strain ET solution.
For long cylindrical shells with free side edges, where L ≫ R ≫ h, the
1-D models provide accurate benchmark solutions at the shell’s mid-cross-
section. This is due to the negligible influence of the free side edge boundary
effects in the central region, resulting in almost translational invariance along
the x1 -direction. Consequently, a plane strain state can be assumed there,

27
 6
1D-CST
1D-RST
4 2D-RST (L=10R)
2D-RST (L=100R)
3D-ET
2

0
0 1 2 3 4 5 6 7 8 9

1D-CST
1 1D-RST
2D-RST (L=10R)
0 2D-RST (L=100R)
3D-ET

−1

0 1 2 3 4 5 6 7 8 9

Figure 13: Normal displacement ū and rescaled rotation angle ψ̄2̄ of a semi-cylindrical
shell under external pressure with free side edges and clamped bottom edges (second case)
as function of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 3).

enabling the application of the arc-like model (1D-RST). Similarly, the 3D-
ET simplifies to a 2D problem for the half-ring in the shell’s central portion.
Figure 12 illustrates the mean normal displacement ū and rescaled ro-
tation angle ψ̄2̄ as functions of the rescaled circumferential coordinate x̄2 =
π R̄θ, calculated using the 2D-RST at the mid-section (x̄1 = L̄/2) for clamped
bottom edges. Comparisons are made with 1D-CST, 1D-RST, and plane
strain ET solutions. By varying L̄ (setting it to 10R̄ and 100R̄) while main-
taining ν = 0.3, R̄ = 10, and p = 1, we observe that the 2D-RST solution
converges to the 1D-RST solution as L̄ increases, aligning closely with the
3D-ET results.
For moderately thick shells (R̄ = 3) with clamped bottom edges, Figure
13 shows similar trends. The 2D-RST solutions at varying L̄ values confirm
convergence to the 1D-RST solution. Comparisons with CST and ET demon-
strate the continued applicability of RST, while highlighting the inadequacy
of CST in this thickness regime.
We also examine semi-cylindrical shells with free side edges and simply
supported bottom edges. As before, the free side edge effects are negligible in
the central region for L̄ ≫ R̄. Numerical results presented in Figures 14 and

28
 60
1D-CST
1D-RST
40 2D-RST (L=10R)
2D-RST (L=100R)
3D-ET
20

0
0 5 10 15 20 25 30

1D-CST
5 1D-RST
2D-RST (L=10R)
0 2D-RST (L=100R)
3D-ET
−5

0 5 10 15 20 25 30

Figure 14: Normal displacement ū and rescaled rotation angle ψ̄2̄ of a semi-cylindrical shell
under external pressure with free side edges and simply-supported bottom edges (third
case) as function of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 10).

15 validate this observation. These figures further illustrate the 2D-RST’s ap-
plicability to both thin and moderately thick shells, while reinforcing CST’s
limitations for moderately thick shells. Importantly, the FE implementa-
tion of the 2D-RST demonstrates both locking-free behavior and asymptotic
accuracy across all numerical simulations.

5.3. A structure with two plates and a spherical cap under external pressure
To showcase the applicability of developed shell formulation in practical
application, the deformation of a structure consisting of two orthogonal plates
and a spherical cap subjected to external pressure, is investigated. The
structure is defined within a Cartesian coordinate system (z1 , z2 , z3 ), with
the (z1 , z2 )-plane coinciding with the base. The mid-surface of the first plate,
P1 , is described by:
q
P1 = {(z1 , 0, z3 ) | |z1 | ≤ D, 0 ≤ z3 ≤ c − R2 − z12 }, (94)
√
where c = H + R2 − D2 (see Figure 16). Note that the full width of P1
extends to |z1 | ≤ D + h/2, but due to h ≪ D, we focus on the portion

29
 6
1D-CST
1D-RST
4 2D-RST (L=10R)
2D-RST (L=100R)
3D-ET
2

0
0 1 2 3 4 5 6 7 8 9

2
1D-CST
1 1D-RST
2D-RST (L=10R)
0 2D-RST (L=100R)
3D-ET
−1

−2
0 1 2 3 4 5 6 7 8 9

Figure 15: Normal displacement u and rotation ψ of a semi-cylindrical shell under external
pressure with free side edges and simply-supported bottom edges (third case) as function
of the rescaled circumferential coordinate x̄2 = π R̄θ (R̄ = 3).

intersecting the mid-surface. Similarly, the mid-surface of the second plate,

P2 , is: q
P2 = {(0, z2 , z3 ) | |z2 | ≤ D, 0 ≤ z3 ≤ c − R2 − z22 }. (95)
The spherical cap’s mid-surface, P3 , is:
q
P3 = {(z1 , z2 , z3 ) | z12 + z22 ≤ D2 , z3 = c − R2 − z12 − z22 }. (96)

We define Rϕ (azimuthal) and Rθ (polar) as the curvilinear coordinates on

P3 . All components are made of the same material, with a uniform thickness
h, where h ≪ D. The overall mid-surface of the structure is S = P1 ∪P2 ∪P3 .
A constant pressure p is applied to the cap’s top surface, while the plates’
bottom edges are clamped. The remaining boundaries are traction-free.
The energy functional remains consistent with Eq. (11), with Φ(γαβ , ραβ , φα )
defined by Eqs. (12) and (13). The external work is given by:
Z h
σh β 1 i
Aext = −p(1 + h/R)u + pγβ + σh2 pρββ dω . (97)
P3 2 10

30
Figure 16: The spherical cap structure and its (z1 , z3 )-cross section. The unit of dimension
is in cm.

The mid-surfaces P1 , P2 , and P3 intersect along three lines, where continuity

of displacement and rotation vectors is enforced. Rotations about surface
normals are constrained, considering only tangential components. At the
triple point (0, 0, c − R), where all surfaces meet, the rotation vector is zero.

Figure 17: Mesh (left) and z3 -displacement of the spherical cap structure (right). The
deformed structure is visualized with scale factor = 3.

31
 The structure is discretized using 1082 Q9 elements with material prop-
erties: E = 2 × 107 N/cm2 , ν = 0.3, and h = 1 cm. The applied pressure
is p = 103 N/cm2 . Figure 17 depicts the deformed structure and its z3 -
displacement. Figure 18 illustrates the rotation vector of the spherical cap
projected onto tϕ (azimuthal) and tθ (polar), showing that the azimuthal
rotation vanishes along the intersections of the shell and plates.

Figure 18: Rotation angle projected on tϕ (left) and tθ (right) of the spherical cap struc-
ture.

6. Conclusion
This work presents a novel rescaled formulation of refined shell theory that
is inherently free from membrane and shear locking. Combined with an isoge-
ometric finite element implementation, this formulation achieves asymptotic
accuracy in shell structure analysis. This combined achievement - asymp-
totic accuracy and freedom from locking - is key contributions of our paper.
Unlike prior methods that often suggest numerical techniques to gain sta-
bility or accuracy, hence complicates the formulation further, our approach
is mathematical based, concise and rigorous. This enables a straightforward
implementation in the numerical code and yields an efficient pathway for
analysis of shell structures, especially where accurate stress/displacement
prediction in thin/moderately thick shells is crucial. Last but not least, our
formulation still relies on the linear kinematics assumption. Future work will
develop accurate and locking-free FE implementations for nonlinear refined
shell theory, including applications to buckling analysis and shell theories for
complex materials.

32
References
[1] A.L. Goldenveizer, Theory of elastic thin shells: solid and structural
mechanics, vol. 2, Elsevier, 2014.

[2] V.L. Berdichevsky, Variational Principles of Continuum Mechanics,

Springer, 2009.

[3] J.L. Sanders, An improved fisrt-approximation theory for thin shells,

NASA Tech. Rep., 1 (1959).

[4] W.T. Koiter, A consistent first approximation in the general theory of

thin elastic shells, in: Proceedings of the IUTAM Symposium on the
Theory of Thin Elastic Shells, North-Holland Amsterdam, 1960, pp.
12–33.

[5] V.L. Berdichevsky, V. Misyura, Effect of accuracy loss in classical shell

theory, Journal of Applied Mechanics, 59 (2S) (1992) S217–S223.

[6] V.L. Berdichevsky, Variational-asymptotic method of constructing a

theory of shells, Journal of Applied Mathematics and Mechanics, 43
(4) (1979) 711–736.

[7] V.L. Berdichevsky, Variational-asymptotic method of constructing the

nonlinear shell theory, in: W.T. Koiter, G.K. Mikhailov (Eds.), Pro-
ceedings of IUTAM Symposium on Shell Theory, North-Holland, Ams-
terdam, 1979, pp. 137–161.

[8] E. Reissner, The effect of transverse shear deformation on the bending

of elastic plates, J. Appl. Mech., 12 (2) (1945) 69–77.

[9] V.G. Sutyrin, Derivation of plate theory accounting asymptotically cor-

rect shear deformation, J. Appl. Mech., 64 (1997) 905–915.

[10] W. Yu, Mathematical construction of a Reissner–Mindlin plate theory

for composite laminates, International Journal of Solids and Structures,
42 (26) (2005) 6680–6699.

[11] Z. Yifeng, C. Lei, W. Yu, Variational asymptotic modeling of the ther-

momechanical behavior of composite cylindrical shells, Composite Struc-
tures, 94 (3) (2012) 1023–1031.

33
[12] Z. Yifeng, C. Lei, W. Yu, Variational asymptotic modeling of the mul-
tilayer functionally graded cylindrical shells, Composite Structures, 94
(3) (2012) 966–976.
[13] R.G. Burela, D. Harursampath, Vam applied to dimensional reduction
of non-linear hyperelastic plates, International Journal of Engineering
Science, 59 (2012) 90–102.
[14] R.G. Burela, D. Harursampath, Asymptotically-accurate nonlinear hy-
perelastic shell constitutive model using variational asymptotic method,
in: H. Altenbach et al. (Eds.), Recent Developments in the Theory of
Shells, Springer, 2019, pp. 135–156.
[15] S.K. Bhadoria, R.G. Burela, Analytical and computational study of
compressible neo-hookean model using vam for two types of global warp-
ing constraints, International Journal of Non-Linear Mechanics, 160
(2024) 104652.
[16] S.K. Bhadoria, R.G. Burela, Asymptotically correct 3d displacement
of the Mooney–Rivlin model using vam, Thin-Walled Structures, 195
(2024) 111358.
[17] K.C. Le, An asymptotically exact first-order shear deformation theory
for functionally graded plates, International Journal of Engineering Sci-
ence, 190 (2023) 103875.
[18] K.C. Le, B.D. Nguyen, Polygonization: Theory and comparison with
experiments, International Journal of Engineering Science, 59 (2012)
211–218.
[19] K.C. Le, B.D. Nguyen, On bending of single crystal beam with contin-
uously distributed dislocations, International Journal of Plasticity, 48
(2013) 152–167.
[20] K.C. Le, L.T.K. Nguyen, Energy methods in dynamics, Springer, 2014.
[21] K.C. Le, J.-H. Yi, An asymptotically exact theory of smart sandwich
shells, International Journal of Engineering Science, 106 (2016) 179–198.
[22] K.C. Le, An asymptotically exact theory of functionally graded piezo-
electric shells, International Journal of Engineering Science, 112 (2017)
42–62.

34
[23] K.C. Le, Introduction to Micromechanics, Nova Science, 2 edition, 2020.

[24] K.C. Le, T.M. Tran, Asymptotically exact theory of fiber-reinforced

composite beams, Composite Structures, 244 (2020) 112279.

[25] K.C. Le, T.M. Tran, Asymptotically exact theory of functionally graded
elastic beams, International Journal of Engineering Science, 209 (2025)
104214.

[26] T. Belytschko, H. Stolarski, W.K. Liu, N. Carpenter, J.S.J. Ong, Stress

projection for membrane and shear locking in shell finite elements, Com-
puter Methods in Applied Mechanics and Engineering, 51 (1-3) (1985)
221–258.

[27] M.L. Bucalem, K.-J. Bathe, Finite element analysis of shell structures,
Archives of Computational Methods in Engineering, 4 (1997) 3–61.

[28] K.-U. Bletzinger, M. Bischoff, E. Ramm, A unified approach for shear-

locking-free triangular and rectangular shell finite elements, Computers
& Structures, 75 (3) (2000) 321–334.

[29] K.C. Le, Vibrations of Shells and Rods, Springer, Berlin, 1999.

[30] O.C. Zienkiewicz, R.L. Taylor, J.M. Too, Reduced integration technique
in general analysis of plates and shells, International Journal for Numer-
ical Methods in Engineering, 3 (2) (1971) 275–290.

[31] E. Onate, E. Hinton, N. Glover, Techniques for improving the perfor-

mance of Ahmad shell elements, in: Proc. Int. Conf. on Applied Numer-
ical Modelling, Pentech Press, Madrid, 1978, pp. 30–53.

[32] T.J.R. Hughes, M. Cohen, M. Haroun, Reduced and selective integration

techniques in the finite element analysis of plates, Nuclear Engineering
and Design, 46 (1) (1978) 203–222.

[33] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic
Finite Element Analysis, Courier Corporation, 2012.

[34] T.J.R. Hughes, T. Tezduyar, Finite elements based upon Mindlin plate
theory with particular reference to the four-node bilinear isoparametric
element, J. Appl. Mech., 48 (1981) 587–596.

35
[35] A.F. Saleeb, T.Y. Chang, W. Graf, S. Yingyeunyong, A hybrid/mixed
model for non-linear shell analysis and its applications to large-rotation
problems, International Journal for Numerical Methods in Engineering,
29 (2) (1990) 407–446.

[36] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and
the method of incompatible modes, International Journal for Numerical
Methods in Engineering, 29 (8) (1990) 1595–1638.

[37] N. Büchter, E. Ramm, D. Roehl, Three-dimensional extension of non-

linear shell formulation based on the enhanced assumed strain concept,
International Journal for Numerical Methods in Engineering, 37 (15)
(1994) 2551–2568.

[38] F. Koschnick, M. Bischoff, N. Camprubı́, K.-U. Bletzinger, The discrete

strain gap method and membrane locking, Computer Methods in Ap-
plied Mechanics and Engineering, 194 (21-24) (2005) 2444–2463.

[39] T. Nguyen-Thoi, P. Phung-Van, C. Thai-Hoang, H. Nguyen-Xuan, A

cell-based smoothed discrete shear gap method (cs-dsg3) using triangu-
lar elements for static and free vibration analyses of shell structures,
International Journal of Mechanical Sciences, 74 (2013) 32–45.

[40] J.F. Caseiro, R.A.F. Valente, A. Reali, J. Kiendl, F. Auricchio, R.J.

Alves de Sousa, On the assumed natural strain method to alleviate
locking in solid-shell NURBS-based finite elements, Computational Me-
chanics, 53 (2014) 1341–1353.

[41] B. Oesterle, E. Ramm, M. Bischoff, A shear deformable, rotation-free

isogeometric shell formulation, Computer Methods in Applied Mechan-
ics and Engineering, 307 (2016) 235–255.

[42] M. Lavrenčič, B. Brank, Hybrid-mixed low-order finite elements for ge-

ometrically exact shell models: Overview and comparison, Archives of
Computational Methods in Engineering, 28 (2021) 3917–3951.

[43] K.C. Le, H.-G. Bui, Asymptotically accurate and locking-free finite ele-
ment implementation of first order shear deformation theory for plates,
Computers & Structures, 298 (2024) 107387.

36
[44] J. Kiendl, Y. Bazilevs, M.C. Hsu, R. Wüchner, K.U. Bletzinger, The
bending strip method for isogeometric analysis of Kirchhoff–Love shell
structures comprised of multiple patches, Computer Methods in Applied
Mechanics and Engineering, 199 (37) (2010) 2403–2416.

[45] M.J. Borden, M.A. Scott, J.A. Evans, T.J.R. Hughes, Isogeometric finite
element data structures based on Bézier extraction of NURBS, Interna-
tional Journal for Numerical Methods in Engineering, 87 (1-5) (2011)
15–47.

[46] P. Dadvand, R. Rossi, E. Oñate, An Object-oriented Environment for

Developing Finite Element Codes for Multi-disciplinary Applications,
Archives of Computational Methods in Engineering, 17 (3) (2010) 253–
297.

[47] A.E.H. Love, A treatise on the mathematical theory of elasticity, Cam-

bridge University Press, 2013.

37