# Nonlinear static isogeometric analysis of arbitrarily curved Kirchhoff-Love
shells

G. Radenković  Faculty of Civil Engineering, University of Belgrade, Bulevar
kralja Aleksandra 73, 11000 Belgrade, Serbia  University of Banja Luka,
Faculty of Architecture, Civil Engineering and Geodesy, Department of
Mechanics and Theory of Structures, 78000 Banja Luka, Bosnia and Herzegovina,
<EMAIL_ADDRESS><EMAIL_ADDRESS>A. Borković
University of Banja Luka, Faculty of Architecture, Civil Engineering and
Geodesy, Department of Mechanics and Theory of Structures, 78000 Banja Luka,
Bosnia and Herzegovina<EMAIL_ADDRESS><EMAIL_ADDRESS>Institute of Applied Mechanics, Graz University of
Technology, Technikerstraße 4/II, 8010 Graz, Austria  B. Marussig  Institute
of Applied Mechanics, Graz University of Technology, Technikerstraße 4/II,
8010 Graz, Austria

## Abstract

The geometrically rigorous nonlinear analysis of elastic shells is considered
in the context of finite, but small, strain theory. The research is focused on
the introduction of the full shell metric and examination of its influence on
the nonlinear structural response. The exact relation between the reference
and equidistant strains is employed and the complete analytic elastic
constitutive relation between energetically conjugated forces and strains is
derived. Utilizing these strict relations, the geometric stiffness matrix is
derived explicitly by the variation of the unknown metric. Moreover, a compact
form of this matrix is presented. Despite the linear displacement distribution
due to the Kirchhoff-Love hypothesis, a nonlinear strain distribution arises
along the shell thickness. This fact is sometimes disregarded for the
nonlinear analysis of thin shells based on the initial geometry, thereby
ignoring the strong curviness of a shell at some subsequent configuration. We
show that the curviness of a shell at each configuration determines the
appropriate shell formulation. For shells that become strongly curved at some
configurations during deformation, the nonlinear distribution of strain
throughout the thickness must be considered in order to obtain accurate
results. We investigate four computational models: one based on the full
analytical constitutive relation, and three simplified ones. Robustness,
accuracy and relative efficiency of the presented formulation are examined via
selected numerical experiments. Our main finding is that the employment of the
full metric is often required when the complete response of the shells is
sought, even for the initially thin shells. Finally, the simplified model that
provided the best balance between efficiency and accuracy is suggested for the
nonlinear analysis of strongly curved shells.

Keywords: strongly curved Kirchhoff-Love shell; nonlinear analysis; reciprocal
shift tensor; analytical constitutive relation

##  1 Introduction

Thin shells are essential components of various modern-day engineering
structures and their analysis is a classical topic of applied mechanics, [1].
Due to permanent advances in computational capabilities and increasing
requirements from industry, the development of more accurate mathematical and
mechanical models is an everlasting topic. The nonlinear analysis of thin
shells is often performed assuming the linear elastic stress-strain relation
which is appropriate when these structures exhibit finite displacements and
rotations but small strains. The analysis of arbitrarily curved elastic shells
is often simplified by neglecting higher-order terms in the metric and the
strain-displacement relation.

The seminal work by Naghdi, [1], is recommended as a primary source for shell
theory, while a comprehensive overview of contemporary achievements in shell
theories and finite element applications can be found in [2]. Besides the
classical finite element method, [3, 4], several other efficient techniques
exist, such as the meshfree approaches, [5, 6], the generalized quadrature
method, [7], and the finite strip method, [8]. Also, there is ongoing work on
the enhancements of shell theories, see e.g. [9]. In the following, we will
focus on the contributions within the scope of isogeometric analysis (IGA)
[10], since it is the discretization technique employed herein. IGA utilizes
CAD basis functions (e.g. NURBS) for finite element analysis. They are
employed for the description of the geometry as well as the physical fields.
Many useful computational properties stem from this approach, the most
prominent one is increased interelement continuity, which provided new impetus
for the development of thin shell formulations.

The first IGA Kirchhoff-Love (KL) shell element was developed in [11] where
the required continuity between KL elements is easily handled within the scope
of IGA. Since the second-order derivatives of displacements are present in the
variational equations of the KL theory, the approximative functions must have
at least C1 interelement continuity. Using NURBS basis functions enables
C$p$-1 interelement continuity where $p$ is the order of the basis. This
alleviates locking issues which often occur due to low order discretizations.
Furthermore, higher continuity allows the implementation of the KL shell
models without rotational DOFs. A large deformation rotation-free formulation
was presented in [12] where different choices for defining the shell normal
vector are discussed. The formulation for compressible and incompressible
hyperelastic thin shells that can use general 3D constitutive models is
developed in [13]. Numerical tests, including structural dynamics simulations
of a bioprosthetic heart valve, emphasize the good performance and
applicability of the present method. An extended IGA for the analysis of
through-the-thickness cracks is developed in [14] while the application of the
subdivision algorithm is discussed in [15]. In all these papers, a linear
strain distribution is assumed.

Nonlinear buckling of thin shells is considered in [16] where a reduced-order
model based on Koiter perturbation is utilized as a predictor, leading to an
improved efficiency. An isogeometric KL shell formulation is extended to the
large strain elasto-plastic analysis using a stress-based approach in [17].
Furthermore, a reduced model with a simplified strain measure is considered in
[18] and no significant loss of accuracy is detected, even for very large
displacements and composite structures. The authors in [19] analyze shell
models extracted from 3D continua and those directly formulated in a 2D
manifold form in which the 2D Koiter model gives excellent performance.
Notably, a new type of instability named snap-backward is detected in [20]
while examining the postbuckling behavior of imperfect shells.

Although single-patch NURBS surfaces are of primary interest here, let us note
that the construction of more complex geometries is an on-going topic in IGA.
For example, the assembly of structures consisting of multiple shell patches
requires special attention. The bending strip method is introduced in [21]
while the coupling without penalty or stabilization parameters is suggested in
[22]. Furthermore, other CAD basis functions, besides NURBS, may also be used.
For instance, the application of subdivision surfaces [23] or T-splines [24]
was discussed. The latter shows that various trimmed geometries can be
converted to T-splines and used in the analysis. Trimmed geometries are
omnipresent in most engineering CAD models, but they also impose several
computational challenges [25]. In the context of isogeometric KL analysis, the
adequate integration of trimmed surfaces and the application of boundary
conditions using penalty approaches [26], [27] or Nitsche’s method [28, 29,
30] were of particular interest.

While this research deals exclusively with KL shells, we draw attention to the
fact that the first IGA-based thick shell formulation within the scope of
Mindlin-Reissner theory is presented in [31] while implicit dynamics are
considered in [32]. Since the rotation field is independent of the
displacement field, Lagrange-based finite elements with C0 interelement
continuity are well-suited for these shells, [33]. Nevertheless, IGA provides
important enhancements for the modeling of Mindlin-Reissner shells, such as
the exact initial geometry and increased smoothness of both rotational and
displacement fields, [34, 35]. Interesting work is given in [36], where two
hierarchic finite element formulations for geometrically nonlinear shell
analysis including the effects of transverse shear are presented. Both methods
combine a fully nonlinear KL shell model with hierarchically added linearized
transverse shear components. The underlying assumption is that in most
practical applications the transverse shear angles are small, which is
confirmed by various numerical experiments. The hierarchic construction
results in an additive strain decomposition into two parts: one resulting from
the membrane and bending deformation and the other from transverse shear. The
authors conclude that large rotations of Mindlin-Reissner shells can be
modeled without parametrizing the corresponding rotation tensor.

It is often reckoned that KL theory is appropriate for thin shells with
$h/R\leq 0.05$, where $R$ is the minimum shell radius of curvature and $h$ is
the thickness, [37]. This limit is not strict, however, and it is a subject to
be further investigated, [1]. Unlike strongly curved beams, [38, 39, 40, 41,
42], strongly curved shells are rarely considered in the literature. The
influences of the approximation of shift tensor, nonlinear strain distribution
and integration through the thickness of thick shells are discussed in [2].
The classification of thin shell models with emphasis on strongly curved
shells is considered in [43] using the asymptotic expansion technique. The
same method is utilized in [44] for the study of thick KL shells. This
relative lack of literature is attributed to the fact that the KL theory is
appropriate for the thin shells. The ratio $h/R$ is already mentioned as the
one type of a measure of shell slenderness. We suggest a more comprehensive
measure of slenderness, the product $Kh$, where $K$ is the absolute value of
the trace of the curvature tensor. We will refer to this quantity as the
curviness of a shell, analogously to beams, [41]. Therefore, the construction
arbitrarily curved refers to the curviness of a shell and not merely to its
curvature. Since we are dealing with nonlinear analysis, this parameter will
be considered during the whole deformation process, [2]. As strict limits do
not exist, we will refer to the shells with curviness $Kh\leq 0.05$ as small-
curvature shells, while all the others belong to the class of strongly curved,
i.e. large-curvature shells. It is emphasized that this measure has a strictly
local character.

The present research is based on [45] but with significant theoretical
improvements and thorough numerical investigation. In order to study
arbitrarily curved KL shells, the complete metric of the shell is considered
and exact relations between strain at the equidistant surface and reference
strains are employed. Here, the equidistant surface refers to the surface
which consists of points that are equally distanced from the shell midsurface,
measured along the midsurface normal. Analytical integration of the virtual
power along the thickness is complex due to the presence of higher order
rational functions in the constitutive relation, as noted in [19, 2], although
the mathematics is well-posed, [1]. This integration is the main theoretical
contribution of the present research. We propose to utilize the analytical
expression for the reciprocal shift tensor. The full metric can be considered
numerically, as well. Analytical integration, however, allows us to define
reduced constitutive models and compare them with the complete one. The
simplification is not applied in the strain-displacement relation for these
reduced models, but in the constitution itself. We show that the simplest
constitutive plate model that decouples bending and membrane actions gives
reasonably accurate results for small-curvature shells, whereas more involved
coupling is required for strongly curved shells. By utilizing this novel
integration and equations of motion derived in [45], the geometric stiffness
matrix is obtained explicitly by the variation of the current, unknown metric.
This includes the variation of Christoffel symbols, the determinant of the
metric tensor and the normal. Considering the complexity of these expressions,
elegant and compact form of the geometric stiffness matrix is obtained.
Finally, the validity of the present formulation is confirmed via several
numerical examples. By considering the efficiency analysis and accuracy of
results, the most involved reduced model is recommended as the favorable one.

The paper is structured as follows. The next section deals with the metric of
a shell with the emphasis on the reciprocal shift tensor. Subsequently, the
strict definition of the shell’s kinematics is introduced and the rigorous
relation between equidistant and reference strains is highlighted. The virtual
power of an arbitrarily curved shell is postulated and discretized in the
fourth section. It is shown that through-the-thickness integration can be
performed analytically. Thorough numerical experiments are given in the fifth
section and these are followed by the conclusions.

##  2 Metric of the shell continuum

A detailed and rigorous definition of the shell metric is presented in this
section. The classical KL assumption states that a straight line element is
rigid and remains perpendicular to the shell’s midsurface in the deformed
configuration. This assumption leads to the degeneration of a 3D continuum
model into an arbitrarily shaped surface. Note that another approach also
exists in which the surface is primarily observed and the so-called directors
are introduced at each point on the surface, [1, 2]. The material law that
follows from this approach, however, neglects the influence of curvature of
the shell, restricting its applicability solely to small-curvature shells,
[2]. The present analysis is performed with respect to the convective frame of
reference while the complete shell kinematics is defined by the translation of
shell midsurface.

Remark. It is interesting to note that the convective coordinate frame is not
an arbitrary choice among the spatial coordinate frames. We must find
coordinates and base vectors in order to fully describe the metric of the
spatial configuration. This is a significant issue since we are essentially
dealing with one equation with two unknowns. The solution is to keep one
quantity constant, either the coordinates or the base vectors. It is a
standard practice in computational mechanics to keep the coordinates constant,
which gives rise to the convective description, where every particle at each
configuration has the same coordinates.

Regarding the notation, boldface lowercase letters are used for vectors, while
uppercase is used for tensors and matrices. An overbar designates the quantity
at the equidistant surface of the shell while the asterisk sign is used to
designate deformed configuration. Finally, the hat symbol designates the local
component of a vector, with respect to the curvilinear coordinates. The direct
and index notations are applied simultaneously, depending on the context, and
the standard summation convention is adopted. The Greek index letters take
values of 1 and 2 while the Latin indices take values of 1, 2, and 3. Partial
and covariant derivatives with respect to the mth coordinate are designated
with $()_{,m}$ and $()_{|m}$, respectively.

The elaboration on the NURBS-based IGA surface modeling is excluded for
brevity since it is readily available in the literature. For a detailed
insight into the IGA, references [10] and [11] are recommended.

###  2.1 Metric of the midsurface

The midsurface of a shell with respect to the Cartesian coordinates is
determined by its position vector $\textbf{r}=\\{x=x^{1},y=x^{2},z=x^{3}\\}$
which will be described as:

$\textbf{r}=x^{m}\textbf{i}_{m}=\sum\limits_{I=1}^{N}\sum\limits_{J=1}^{M}R_{IJ}(\xi,\eta)\textbf{r}_{IJ},\quad
x^{m}=\sum\limits_{I=1}^{N}\sum\limits_{J=1}^{M}R_{IJ}(\theta^{\alpha})x^{m}_{IJ},$
(1)

where $R_{IJ}$ are bivariate NURBS basis functions,
$\textbf{r}_{IJ}=\\{x_{IJ},y_{IJ},z_{IJ}\\}$ are the position vectors of the
control points $IJ$, while $N$ and $M$ are the total numbers of control points
along the $\xi$ and $\eta$ directions, respectively, [10]. Furthermore,
$\textbf{i}^{m}=\textbf{i}_{m}$ are the base vectors of the Cartesian
coordinate system while $\theta^{\alpha}$ are coordinates $\xi$ and $\eta$
such that $\theta^{1}=\xi$ and $\theta^{2}=\eta$, see Fig. 1.

Figure 1: The mapping between (center) the parameter domain $\mathscr{P}$,
(left) reference and (right) deformed configurations of a Kirchhoff-Love
shell. The shell boundaries are colored with solid green lines, the
equidistant surface is indicated by dashed red lines, and the midsurface is
colored in gray with black boundaries. Quantities at the deformed
configuration are denoted with an asterisk while those at the equidistant
surfaces have an overbar.

Tangent vectors of the midsurface are given by

$\textbf{g}_{\alpha}=\frac{\partial\textbf{r}}{\partial\theta^{\alpha}}=x^{m}_{,\alpha}\textbf{i}_{m},\quad
x^{m}_{,\alpha}=\sum\limits_{I=1}^{N}\sum\limits_{J=1}^{M}R_{IJ,\alpha}x^{m}_{IJ},$
(2)

where $x^{m}_{,\alpha}$ are the partial derivatives of the Cartesian
components of the midsurface with respect to the $\theta^{\alpha}$
coordinates. In order to describe the complete geometry of a shell, the third
coordinate $\theta^{3}=\zeta$ is required. Due to the KL hypothesis, this
coordinate line is rigid, straight and perpendicular to the $\theta^{\alpha}$
coordinate lines. It follows that the base vector $\textbf{g}_{3}$ and its
reciprocal counterpart $\textbf{g}^{3}$ are the same. They have the unit
length and represent the normal at each point of the shell midsurface:

$\textbf{g}^{3}=\textbf{g}_{3}=\frac{1}{\sqrt{g}}(\textbf{g}_{1}\times\textbf{g}_{2})=x^{m}_{,3}\textbf{i}_{m},\quad
x^{m}_{,3}=x_{m,3}=\frac{1}{\sqrt{g}}x^{k}_{,1}x^{l}_{,2}e_{klm},$ (3)

where $e_{klm}$ is permutation symbol while $g$ is the determinant of the
metric tensor of shell midsurface:

$g_{ij}=\begin{bmatrix}x^{m}_{,1}x_{m,1}&x^{m}_{,1}x_{m,2}&0\\\
x^{m}_{,2}x_{m,1}&x^{m}_{,2}x_{m,2}&0\\\
0&0&1\end{bmatrix}=\begin{bmatrix}g_{11}&g_{12}&0\\\ g_{21}&g_{22}&0\\\
0&0&1\end{bmatrix},\;g_{12}=g_{21},\;g=\det g_{ij}=g_{11}g_{22}-{g_{12}}^{2}.$
(4)

The reciprocal metric tensor of midsurface is:

$g^{ij}={g_{ij}}^{-1}=\frac{1}{g}\begin{bmatrix}g_{22}&-g_{12}&0\\\
-g_{12}&g_{11}&0\\\ 0&0&g\end{bmatrix}=\begin{bmatrix}g^{11}&g^{12}&0\\\
g^{21}&g^{22}&0\\\ 0&0&1\end{bmatrix},\;\det g^{ij}=1/g,$ (5)

while the reciprocal tangent vectors are:

$\displaystyle\textbf{g}^{1}$
$\displaystyle=\frac{1}{\sqrt{g}}(\textbf{g}_{2}\times\textbf{g}_{3})=x^{,1}_{n}\textbf{i}^{n},\quad
x^{,1}_{n}=\frac{1}{\sqrt{g}}x^{m}_{,2}x^{k}_{,3}e_{mkn},$ (6)
$\displaystyle\textbf{g}^{2}$
$\displaystyle=\frac{1}{\sqrt{g}}(\textbf{g}_{3}\times\textbf{g}_{1})=x^{,2}_{n}\textbf{i}^{n},\quad
x^{,2}_{n}=\frac{1}{\sqrt{g}}x^{m}_{,3}x^{k}_{,1}e_{mkn}.$

In order to completely define the metric of the midsurface, the Christoffel
symbols are required. They stem from the differentiation of the base vectors
with respect to the curvilinear coordinates, [1]:

$\textbf{g}_{\alpha,\beta}=x^{m}_{,\alpha\beta}\textbf{i}_{m}=\Gamma_{\alpha\beta
k}\textbf{g}^{k}=\Gamma^{n}_{\alpha\beta}\textbf{g}_{n}\Rightarrow\begin{cases}\begin{aligned}
\Gamma_{\alpha\beta k}&=x^{m}_{,\alpha\beta}x_{m,k}\\\
\Gamma^{n}_{\alpha\beta}&=x^{m}_{,\alpha\beta}x^{,n}_{m}\end{aligned}\end{cases}\Rightarrow\Gamma^{n}_{\alpha\beta}=g^{nk}\Gamma_{\alpha\beta
k}.$ (7)

Here, $\Gamma_{\alpha\beta k}$ and $\Gamma^{n}_{\alpha\beta}$ are the
Christoffel symbols of the first and the second kind, respectively. The
components of the differentiated base vectors with respect to the Cartesian
coordinates are:

$x^{m}_{,\alpha\beta}=\Gamma^{n}_{\alpha\beta}x^{m}_{,n}=\Gamma_{\alpha\beta
n}x^{m,n}.$ (8)

The Christoffel symbols of the second kind related to the third coordinate
line $\Gamma^{3}_{\alpha\beta}$ are of special importance since they represent
curvatures of the surface with respect to the coordinate system
$\theta^{\alpha}$. They determine the curvature tensor and it is a common
practice to designate a particular symbol for them, i.e.,
$\Gamma^{3}_{\alpha\beta}=b_{\alpha\beta}$. The components of the curvature
tensor of a surface can be determined by differentiating the normal with
respect to the curvilinear coordinates:

$\textbf{g}_{3,\alpha}=-\Gamma^{\mu}_{3\alpha}\textbf{g}_{\mu}=-b^{\mu}_{\alpha}\textbf{g}_{\mu},\;(\Gamma^{3}_{3\alpha}=0),$
(9)

where $b^{\mu}_{\alpha}$ are mixed components of the curvature tensor which
relate to the covariant ones by $b^{\mu}_{\alpha}=g^{\mu\nu}b_{\nu\alpha}$,
[1].

###  2.2 Metric of an equidistant surface

The position vector of an equidistant surface and its appropriate tangent
vectors, Fig. 1, are:

$\displaystyle\bar{\textbf{r}}$
$\displaystyle=\textbf{r}+\zeta\textbf{g}_{3},$ (10)
$\displaystyle\bar{\textbf{g}}_{\alpha}$
$\displaystyle=\frac{\partial\bar{\textbf{r}}}{\partial\theta^{\alpha}}=\textbf{g}_{\alpha}-\zeta
b^{\nu}_{\alpha}\textbf{g}_{\nu}=(\delta^{\nu}_{\alpha}-\zeta
b^{\nu}_{\alpha})\textbf{g}_{\nu}=\bar{C}^{\nu}_{\alpha}\textbf{g}_{\nu},\;\bar{C}^{\nu}_{\alpha}=\delta^{\nu}_{\alpha}-\zeta
b^{\nu}_{\alpha},$

where $\bar{C}^{\nu}_{\alpha}$ are components of the shift tensor, [46], while
$\delta^{\nu}_{\alpha}$ are components of the Kronecker tensor. This
expression reveals the fact that the normal
$\textbf{g}_{3}=\bar{\textbf{g}}_{3}$ and the tangent vectors of an
equidistant surface are mutually orthogonal:

$\textbf{g}_{3}\cdotp\bar{\textbf{g}}_{\alpha}=\bar{C}^{\nu}_{\alpha}g_{3\nu}=0.$
(11)

Furthermore, the metric tensor and its determinant at an equidistant surface
are:

$\displaystyle\bar{g}_{\alpha\beta}$
$\displaystyle=\bar{\textbf{g}}_{\alpha}\cdotp\bar{\textbf{g}}_{\beta}=\bar{C}^{\mu}_{\alpha}\bar{C}^{\nu}_{\beta}g_{\mu\nu},$
(12) $\displaystyle\bar{g}$
$\displaystyle=\det\bar{\textbf{g}}_{\alpha\beta}=\left(\det\bar{C}^{\mu}_{\alpha}\right)\left(\det\bar{C}^{\nu}_{\beta}\right)g.$

Using the relations

$\det\bar{C}^{\mu}_{\alpha}=\det(\delta^{\mu}_{\alpha}-\zeta
b^{\mu}_{\alpha})=1-\zeta b^{\nu}_{\nu}+\zeta^{2}b,\quad b=\det
b^{\mu}_{\alpha},\quad b^{\nu}_{\nu}=\Tr b^{\mu}_{\alpha},$ (13)

the determinant of the metric tensor at an equidistant surface reduces to:

$\bar{g}=\left(1-\zeta b^{\nu}_{\nu}+\zeta^{2}b\right)^{2}g={g_{0}}^{2}g,\quad
g_{0}=\det\bar{C}^{\mu}_{\alpha},$ (14)

where $g_{0}$ is the determinant of the shift tensor, also known as the
shifter, [19], while $b^{\nu}_{\nu}$ and $b$ are the trace and determinant of
the curvature tensor, respectively. Note that the trace of the curvature
tensor is equal to the sum of the principal curvatures,
$b^{\nu}_{\nu}=K_{1}+K_{2}=2H$, where $H$ is the mean curvature, [46].

For the explicit derivation of the constitutive relation which follows, it is
crucial to establish the analytical relation between the tangent vectors of
the equidistant and middle surfaces. From Eq. (10) it follows that the
reciprocal shift tensor
$C^{\alpha}_{\nu}=\left(\bar{C}^{\nu}_{\alpha}\right)^{-1}$ could be
introduced as well:

$\bar{\textbf{g}}^{\alpha}=C^{\alpha}_{\nu}\textbf{g}^{\nu},$ (15)

and it can be easily numerically calculated. However, we are interested in its
analytical form. Let us utilize the fact that the base vector $\textbf{g}_{3}$
is independent of $\zeta$ coordinate and observe the vector product:

$\bar{\textbf{g}}_{3}\times\bar{\textbf{g}}_{\mu}=\textbf{g}_{3}\times\bar{C}^{\gamma}_{\mu}\textbf{g}_{\gamma}=\sqrt{g}\bar{C}^{\gamma}_{\mu}e_{3\gamma\nu}\textbf{g}^{\nu}.$
(16)

If we multiply this expression with $e^{3\mu\alpha}$:

$\sqrt{\bar{g}}\delta^{\alpha}_{\beta}\bar{\textbf{g}}^{\beta}=\sqrt{g}\bar{C}^{\gamma}_{\mu}\delta^{\mu\alpha}_{\gamma\nu}\textbf{g}^{\nu},$
(17)

the analytical form of the reciprocal shift tensor follows as:

$C^{\alpha}_{\nu}=\frac{1}{g_{0}}\delta^{\alpha\mu}_{\nu\gamma}\bar{C}^{\gamma}_{\mu}=\frac{1}{g_{0}}\left[\left(1-\zeta
b^{\gamma}_{\gamma}\right)\delta^{\alpha}_{\nu}+\zeta
b^{\alpha}_{\nu}\right].$ (18)

The same result can be found in [1] and it allows us to represent the
reciprocal metric tensor at an arbitrary point as an analytic function of the
metric and curvature of the midsurface:

$\bar{g}^{\alpha\beta}=\bar{\textbf{g}}^{\alpha}\bar{\textbf{g}}^{\beta}=C^{\alpha}_{\mu}C^{\beta}_{\nu}g^{\mu\nu}.$
(19)

The tensor $C^{\alpha}_{\mu}C^{\beta}_{\nu}$ is a rational function of the
$\zeta$ coordinate where the numerator and denominator are polynomials of the
$2^{nd}$ and $4^{th}$ degree, respectively. The relation between the metrics
of the equidistant and reference surfaces directly influences the equivalent
strain relations, [2]. In order to simplify the computational model, many
types of approximations are readily introduced by researchers, mostly based on
first- or second-order Taylor approximations of the shift tensor, [47]. The
exact analytical form of the reciprocal shift tensor, on the other hand, is
often disregarded. An interesting discussion on the influence of the exact
shift tensor is given in [2].

Finally, let us define the relations between a differential area at the
equidistant and middle surfaces which follow from (14). If we notice that
these areas are:

$\displaystyle\differential{\bar{\textbf{a}}_{\alpha}}$
$\displaystyle=\differential{\bar{a}_{\alpha}\bar{\textbf{g}}^{\alpha}}=\differential{\zeta}\textbf{g}_{3}\times\bar{\textbf{g}}_{\beta}=\differential{\zeta}\sqrt{\bar{g}}e_{3\beta\alpha}\bar{\textbf{g}}^{\alpha},\quad\differential{\bar{a}_{\alpha}}=\sqrt{\bar{g}}\differential{\zeta}e_{3\beta\alpha},$
(20) $\displaystyle\differential{\textbf{a}_{\alpha}}$
$\displaystyle=\differential{a_{\alpha}\textbf{g}^{\alpha}}=\differential{\zeta}\textbf{g}_{3}\times\textbf{g}_{\beta}=\differential{\zeta}\sqrt{g}e_{3\beta\alpha}\textbf{g}^{\alpha},\quad\differential{a_{\alpha}}=\sqrt{g}\differential{\zeta}e_{3\beta\alpha},$

the required relation follows as:

$\differential{\bar{a}_{\alpha}}=\sqrt{\frac{\bar{g}}{g}}\differential{a_{\alpha}}=g_{0}\differential{a}_{\alpha}.$
(21)

Additionally, the differential volume of a shell is:

$\differential{v}=\left(\differential{\xi}\textbf{g}_{1}\times\differential{\eta}\textbf{g}_{2}\right)\cdotp\differential{\zeta}\textbf{g}_{3}=\sqrt{\bar{g}}\differential{\xi}\differential{\eta}\differential{\zeta}=g_{0}\differential{a}\differential{\zeta},$
(22)

where $\differential{a}=\sqrt{g}\differential{\xi}\differential{\eta}$ is the
differential area of the midsurface.

##  3 Kirchhoff-Love shell theory

After the metric of the shell continuum is defined, the next step is to
introduce some strain measure. For the convective coordinate frame, the
Lagrange strain equals the difference between the current and reference
metrics:

${\epsilon}_{\alpha\beta}=\frac{1}{2}\left({g}^{*}_{\alpha\beta}-{g}_{\alpha\beta}\right).$
(23)

Due to the KL hypothesis, the components of the shear strain rate along the
normal direction vanish, i.e. $d_{\alpha 3}=0$. Consequently, the position of
the line element is completely determined by the translation of the
midsurface, which is the only generalized coordinate of the KL shell. This
fact gives rise to the so-called rotation-free shell theories, [11]. These
formulations are especially efficient for the finite rotation analysis since
the parametrization of the finite rotation variables is not required.

###  3.1 Kinematics of the shell

As indicated in Fig. 1, the position vector of an equidistant surface in the
deformed configuration is:

$\bar{\textbf{r}}^{*}={\textbf{r}}^{*}(\zeta)={\textbf{r}}^{*}+\zeta{\textbf{g}}^{*}_{3},$
(24)

where ${\textbf{r}}^{*}$ is the position vector of the deformed midsurface
while the base vectors ${\textbf{g}}^{*}_{m}$ are calculated analogously to
the Eqs. (2) and (3). The position vector of the deformed midsurface is:

${\textbf{r}}^{*}=\textbf{r}+\textbf{u},$ (25)

where the displacement u is discretized with NURBS, in the same way as the
geometry, Eq. (1):

$\textbf{u}=\sum\limits_{I=1}^{N}\sum\limits_{J=1}^{M}R_{IJ}(\theta^{\alpha})\textbf{u}_{IJ}.$
(26)

$\textbf{u}_{IJ}$ is the vector that consists of displacement components of a
control point $IJ$ with respect to the Cartesian system. For the sake of the
seamless transition to the discrete equation of motion which follows in
Section 4, let us introduce the matrix of basis functions N such that Eq. (26)
can be written as:

$\textbf{u}=\textbf{N}\textbf{q},$ (27)

where:

$\displaystyle\textbf{q}^{\mathsf{T}}$
$\displaystyle=\begin{bmatrix}\textbf{u}^{\mathsf{T}}_{11}&\textbf{u}^{\mathsf{T}}_{12}&...&\textbf{u}^{\mathsf{T}}_{IJ}&...&\textbf{u}^{\mathsf{T}}_{NM}\end{bmatrix},\quad\textbf{u}_{IJ}=\begin{bmatrix}u^{1}_{IJ}&u^{2}_{IJ}&u^{3}_{IJ}\end{bmatrix},$
(28) $\displaystyle\textbf{N}^{\mathsf{T}}$
$\displaystyle=\begin{bmatrix}\textbf{R}^{\mathsf{T}}&\textbf{R}^{\mathsf{T}}&\textbf{R}^{\mathsf{T}}\end{bmatrix},\quad\textbf{R}=\begin{bmatrix}R_{11}&R_{12}&...&R_{IJ}&...&R_{NM}\end{bmatrix}.$

The velocity field is obtained as the material derivative of the displacement
field:

$\textbf{v}=\dot{\textbf{r}}^{*}=\dot{\textbf{u}}=\dot{u}^{n}\textbf{i}_{n}=v^{n}\textbf{i}_{n}=\dot{x}^{*n}\textbf{i}_{n},$
(29)

and it is represented analogously:

$\textbf{v}=\dot{\textbf{u}}=\textbf{N}\dot{\textbf{q}},\;\;\textbf{$\dot{\textbf{q}}$}^{\mathsf{T}}=\begin{bmatrix}\textbf{$\dot{\textbf{u}}_{11}$}^{\mathsf{T}}&\textbf{$\dot{\textbf{u}}_{12}$}^{\mathsf{T}}&...&\textbf{$\dot{\textbf{u}}_{IJ}$}^{\mathsf{T}}&...&\textbf{$\dot{\textbf{u}}_{NM}$}^{\mathsf{T}}\end{bmatrix},\;\;\dot{\textbf{u}}_{IJ}=\begin{bmatrix}v^{1}_{IJ}&v^{2}_{IJ}&v^{3}_{IJ}\end{bmatrix}.$
(30)

The strain rate tensor is the symmetric part of the velocity gradient and it
equals the material derivative of (23):

${d}_{\alpha\beta}=\dot{\epsilon}_{\alpha\beta}=\frac{1}{2}\dot{g}^{*}_{\alpha\beta}=\frac{1}{2}\left(\hat{v}_{\alpha|\beta}+\hat{v}_{\beta|\alpha}\right)=\frac{1}{2}\left(\hat{v}_{\alpha,\beta}+\hat{v}_{\beta,\alpha}-2\hat{v}_{m}{\Gamma}^{*m}_{\alpha\beta}\right),$
(31)

where $\hat{v}_{m}$ are the components of the velocity vector with respect to
the local curvilinear coordinates, which are related to the global ones as:

$\textbf{v}=\hat{v}_{m}\textbf{g}^{*m}=\hat{v}^{m}\textbf{g}^{*}_{m}={v}_{k}\textbf{i}^{k}={v}^{k}\textbf{i}_{k}\Rightarrow\begin{cases}\hat{v}_{m}&={v}_{k}\textbf{i}^{k}\cdotp\textbf{g}^{*}_{m}={x}^{*k}_{,m}{v}_{k}\\\
\hat{v}^{m}&={v}^{k}\textbf{i}_{k}\cdotp\textbf{g}^{*m}={x}^{*,m}_{k}{v}^{k}\end{cases}.$
(32)

The covariant components of the strain rate with respect to the Cartesian
coordinates can be written as:

$\displaystyle{d}_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\left({x}^{*k}_{,\alpha}{v}_{k,\beta}+{x}^{*k}_{,\beta}{v}_{k,\alpha}+2{x}^{*k}_{,\alpha\beta}{v}_{k}-2{\Gamma}^{*m}_{\alpha\beta}{x}^{*n}_{,m}{v}_{n}\right)$
(33)
$\displaystyle=\frac{1}{2}\left({x}^{*k}_{,\alpha}{v}_{k,\beta}+{x}^{*k}_{,\beta}{v}_{k,\alpha}\right)=\frac{1}{2}\left(\textbf{g}^{*}_{\alpha}\cdotp\textbf{v}_{\beta}+\textbf{g}^{*}_{\beta}\cdotp\textbf{v}_{\alpha}\right),$

while the trace of the strain rate tensor is:

${d}^{\alpha}_{\alpha}={g}^{*\alpha\beta}{d}_{\alpha\beta}=\frac{1}{2}\left({x}^{*,\beta}_{n}{v}^{n}_{,\beta}+{x}^{*,\alpha}_{n}{v}^{n}_{,\alpha}\right)={x}^{*,\alpha}_{n}{v}^{n}_{,\alpha}=\textbf{g}^{*\alpha}\cdotp\textbf{v}_{\alpha}.$
(34)

The material derivative of (24) gives the velocity vector of an equidistant
surface:

$\bar{\textbf{v}}=\textbf{v}(\zeta)=\textbf{v}+\zeta\textbf{v}_{3},$ (35)

where the following notation is introduced:

$\bar{\textbf{v}}=\dot{\bar{\textbf{r}}}^{*}=\bar{v}^{n}\textbf{i}_{n}=\dot{\bar{x}}^{*n}\textbf{i}_{n},\;\textbf{v}_{m}=\dot{\textbf{g}}^{*}_{m}={v}^{n}_{,m}\textbf{i}_{n}=\dot{x}^{*n}_{,m}\textbf{i}_{n}.$
(36)

Thus, the vector $\textbf{v}_{3}$ is the material derivative of the deformed
normal:

$\displaystyle\textbf{v}_{3}$
$\displaystyle=\dot{\textbf{g}}^{*}_{3}={v}_{k,3}\textbf{i}^{k}=\frac{1}{\sqrt{{g}^{*}}}\left(\textbf{v}_{1}\times\textbf{g}^{*}_{2}+\textbf{g}^{*}_{1}\times\textbf{v}_{2}-\frac{\dot{g}^{*}}{2\sqrt{{g}^{*}}}\textbf{g}^{*}_{3}\right),$
(37) $\displaystyle{v}_{k,3}$
$\displaystyle=\frac{1}{\sqrt{{g}^{*}}}\left[\left({v}^{m}_{,1}{x}^{*n}_{,2}+{x}^{*m}_{,1}{v}^{n}_{,2}\right){e}_{mnk}-\frac{\dot{g}^{*}}{2\sqrt{{g}^{*}}}{x}^{*}_{k,3}\right],$

and its derivation requires special attention. Firstly, the material
derivative of the determinant of the metric tensor, utilizing Eqs. (4), (31),
and (33), is computed:

$\displaystyle\dot{g}^{*}$
$\displaystyle=2\left({d}_{11}{g}^{*}_{22}+{d}_{22}{g}^{*}_{11}-2{d}_{12}{g}^{*}_{12}\right)$
(38)
$\displaystyle=2\left[\left(\textbf{g}^{*}_{1}\cdotp\textbf{v}_{1}\right){g}^{*}_{22}+\left(\textbf{g}^{*}_{2}\cdotp\textbf{v}_{2}\right){g}^{*}_{11}-2\left(\textbf{g}^{*}_{1}\cdotp\textbf{v}_{2}+\textbf{g}^{*}_{2}\cdotp\textbf{v}_{1}\right){g}^{*}_{12}\right],$

which leads to:

$\frac{\dot{g}^{*}}{2{g}^{*}}={d}_{11}{g}^{*11}+{d}_{12}{g}^{*12}+{d}_{21}{g}^{*21}+{d}_{22}{g}^{*22}={d}^{1}_{1}+{d}^{2}_{2}={d}^{\alpha}_{\alpha}.$
(39)

If we represent the tangent vectors of the deformed configuration as:

$\displaystyle\textbf{g}^{*}_{1}$
$\displaystyle={x}^{*n}_{,1}\textbf{i}_{n}=\sqrt{{g}^{*}}\left(\textbf{g}^{*2}\times\textbf{g}^{*}_{3}\right),\quad{x}^{*n}_{,1}=\sqrt{{g}^{*}}{x}^{*,2}_{k}{x}^{*}_{l,3}{e}^{kln},$
(40) $\displaystyle\textbf{g}^{*}_{2}$
$\displaystyle={x}^{*n}_{,2}\textbf{i}_{n}=\sqrt{{g}^{*}}\left(\textbf{g}^{*}_{3}\times\textbf{g}^{*1}\right),\quad{x}^{*n}_{,2}=\sqrt{{g}^{*}}{x}^{*}_{k,3}{x}^{*,1}_{l}{e}^{kln},$

the first part of Eq. (37) can be rewritten as:

$\displaystyle\frac{1}{\sqrt{{g}^{*}}}\left(\textbf{v}_{1}\times\textbf{g}^{*}_{2}+\textbf{g}^{*}_{1}\times\textbf{v}_{2}\right)$
$\displaystyle=\textbf{v}_{1}\times\left(\textbf{g}^{*}_{3}\times\textbf{g}^{*1}\right)-\textbf{v}_{2}\times\left(\textbf{g}^{*2}\times\textbf{g}^{*}_{3}\right)$
(41)
$\displaystyle=\left(\textbf{g}^{*\alpha}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*}_{3}-\left(\textbf{g}^{*}_{3}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*\alpha}.$

By the insertion of Eqs. (39), (34), and (41) into (37), we can finally obtain
the velocity of the normal:

$\displaystyle\textbf{v}_{3}$
$\displaystyle=\left(\textbf{g}^{*\alpha}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*}_{3}-\left(\textbf{g}^{*}_{3}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*\alpha}-\left(\textbf{g}^{*\alpha}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*}_{3}=-\left(\textbf{g}^{*}_{3}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*\alpha},$
(42) $\displaystyle{v}_{k,3}$
$\displaystyle={x}^{*,\alpha}_{m}{v}^{m}_{,\alpha}{x}^{*}_{k,3}-{x}^{*,\alpha}_{k}{v}^{m}_{,\alpha}{x}^{*}_{m,3}-{x}^{*,\alpha}_{m}{v}^{m}_{,\alpha}{x}^{*}_{k,3}=-{x}^{*,\alpha}_{k}{v}^{m}_{,\alpha}{x}^{*}_{m,3},$

which is the same result as calculated in [1]. Finally, the components of
velocity at an equidistant surface (35) can be represented as:

$\bar{\textbf{v}}=\bar{v}_{n}\textbf{i}^{n}=\textbf{v}-\zeta\left(\textbf{g}^{*}_{3}\cdotp\textbf{v}_{\alpha}\right)\textbf{g}^{*\alpha},\;\bar{v}_{n}={v}_{n}-\zeta{x}^{*,\alpha}_{n}{v}^{m}_{,\alpha}{x}^{*}_{m,3}.$
(43)

Remark. Proving that the line element remains perpendicular to the midsurface
at each configuration is now straightforward. That is, multiplying Eq. (42)
with ${x}^{*k}_{,\alpha}$ gives:

${x}^{*k}_{,\alpha}{v}_{k,3}=-{x}^{*k}_{,\alpha}{x}^{*,\beta}_{k}{v}^{m}_{,\beta}{x}^{*}_{m,3}=-{\delta}^{\beta}_{\alpha}{v}^{m}_{,\beta}{x}^{*}_{m,3}=-{v}^{m}_{,\alpha}{x}^{*}_{m,3},$
(44)

and the shear strain rates in the planes perpendicular to the tangent planes
of the midsurface vanish:

$2{d}_{\alpha 3}={x}^{*n}_{,\alpha}{v}_{n,3}+{x}^{*n}_{,3}{v}_{n,\alpha}=0.$
(45)

Additionally, the strain rate along the normal also vanishes at each
configuration:

${d}_{33}={x}^{*n}_{,3}{v}_{n,3}=-{x}^{*n}_{,3}{x}^{*,\alpha}_{n}{v}^{k}_{,\alpha}{x}^{*}_{k,3}=-{\delta}^{\alpha}_{3}{v}^{k}_{,\alpha}{x}^{*}_{k,3}=0.$
(46)

###  3.2 Strain rate at an equidistant surface

The crucial step for the derivation of structural theories is the definition
of the strain at each point of a structure via chosen reference quantities. It
allows a dimensional reduction from 3D to 2D. The strain rate at an
equidistant surface is, analogously to Eq. (33):

$\bar{d}_{\alpha\beta}={d}_{\alpha\beta}(\zeta)=\frac{1}{2}\left(\bar{\textbf{g}}^{*}_{\alpha}\cdotp\bar{\textbf{v}}_{\beta}+\bar{\textbf{g}}^{*}_{\beta}\cdotp\bar{\textbf{v}}_{\alpha}\right)=\frac{1}{2}\left(\bar{x}^{*n}_{,\alpha}\bar{v}_{n,\beta}+\bar{x}^{*n}_{,\beta}\bar{v}_{n,\alpha}\right),$
(47)

where the components of the velocity gradient at an arbitrary point,
$\bar{v}_{n,\beta}$, follow from Eq. (43):

$\bar{v}_{n,\alpha}={v}_{n,\alpha}+\zeta{v}_{n,3\alpha},$ (48)

while the derivative of the velocity of normal stems from Eq. (42):

$\displaystyle{v}_{n,3\alpha}$
$\displaystyle=-\left({x}^{*,\mu}_{n,\alpha}{v}_{k,\mu}{x}^{*k}_{,3}+{x}^{*,\mu}_{n}{v}_{k,\mu\alpha}{x}^{*k}_{,3}+{x}^{*,\mu}_{n}{v}_{k,\mu}{x}^{*k}_{,3\alpha}\right)$
(49)
$\displaystyle={\Gamma}^{*\mu}_{\alpha\nu}{x}^{*,\nu}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{b}^{*\mu}_{\alpha}{x}^{*,3}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{x}^{*,\mu}_{n}{v}_{k,\mu\alpha}{x}^{*k}_{,3}+{x}^{*,\mu}_{n}{v}_{k,\mu}{b}^{*\nu}_{\alpha}{x}^{*k}_{,\nu}.$

Remark. Utilizing the last two expressions, we can show that the shear strain
rate ${d}_{\alpha 3}$ is zero at an equidistant surface, as well as at the
midsurface:

$\displaystyle 2\bar{d}_{\alpha 3}$
$\displaystyle=\left(\bar{\textbf{g}}^{*}_{\alpha}\cdotp\textbf{v}_{3}+\textbf{g}^{*}_{3}\cdotp\textbf{v}_{\alpha}\right)=\bar{x}^{*n}_{,\alpha}{v}_{n,3}+{x}^{*n}_{,3}\bar{v}_{n,\alpha}$
(50)
$\displaystyle=\left({\delta}^{\mu}_{\alpha}-\zeta{b}^{*\mu}_{\alpha}\right){x}^{*n}_{,\mu}{v}_{n,3}+{x}^{*n}_{,3}\left({v}_{n,\alpha}+\zeta{v}_{n,3\alpha}\right)$
$\displaystyle={x}^{*n}_{,\alpha}{v}_{n,3}-\zeta{b}^{*\mu}_{\alpha}{x}^{*n}_{,\mu}{v}_{n,3}+{x}^{*n}_{,3}{v}_{n,\alpha}$
$\displaystyle\;\;\;\;+\zeta{x}^{*n}_{,3}\left({\Gamma}^{*\mu}_{\alpha\nu}{x}^{*,\nu}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{b}^{*\mu}_{\alpha}{x}^{*,3}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{x}^{*,\mu}_{n}{v}_{k,\mu\alpha}{x}^{*k}_{,3}-{x}^{*,\mu}_{n}{v}_{k,\mu}{b}^{*\nu}_{\alpha}{x}^{*k}_{,\nu}\right)$
$\displaystyle={x}^{*n}_{,\alpha}{v}_{n,3}-\zeta{b}^{*\mu}_{\alpha}{x}^{*n}_{,\mu}{v}_{n,3}+{x}^{*n}_{,3}{v}_{n,\alpha}-\zeta{b}^{*\mu}_{\alpha}{v}_{k,\mu}{x}^{*k}_{,3}$
$\displaystyle={x}^{*n}_{,\alpha}{v}_{n,3}+{x}^{*n}_{,3}{v}_{n,\alpha}-\zeta{b}^{*\mu}_{\alpha}\left({x}^{*n}_{,\mu}{v}_{n,3}+{x}^{*n}_{,3}{v}_{n,\mu}\right)=0,$

which proves that the line segment is perpendicular to the equidistant
surfaces as well as to the midsurface.

Now we have all the ingredients required for the exact relation between
equidistant and reference strain rates. Inserting Eqs. (48) and (49) into Eq.
(47), returns the components of the equidistant strain rate:

$\bar{d}_{\alpha\beta}={A}^{\mu\nu}_{\alpha\beta}{d}_{\mu\nu}-\zeta{B}^{\mu\nu}_{\alpha\beta}\dot{\kappa}_{\mu\nu},$
(51)

where:

$\displaystyle{A}^{\mu\nu}_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\left[\bar{C}^{*\mu}_{\alpha}\left({\delta}^{\nu}_{\beta}+\zeta{b}^{*\nu}_{\beta}\right)+\bar{C}^{*\nu}_{\beta}\left({\delta}^{\mu}_{\alpha}+\zeta{b}^{*\mu}_{\alpha}\right)\right]={\delta}^{\mu}_{\alpha}{\delta}^{\nu}_{\beta}-\zeta^{2}{b}^{*\mu}_{\alpha}{b}^{*\nu}_{\beta},$
(52) $\displaystyle{B}^{\mu\nu}_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\left({\delta}^{\nu}_{\alpha}\bar{C}^{*\mu}_{\beta}+{\delta}^{\mu}_{\beta}\bar{C}^{*\nu}_{\alpha}\right),$

while the tensor of the rate of curvature change of the midsurface,
$\dot{\kappa}_{\alpha\beta}$, is introduced as:

$\displaystyle\dot{\kappa}_{\alpha\beta}$
$\displaystyle=\dot{b}^{*}_{\alpha\beta}=\dot{\textbf{g}}^{*}_{\alpha,\beta}\cdotp\textbf{g}^{*}_{3}+\textbf{g}^{*}_{\alpha,\beta}\cdotp\dot{\textbf{g}}^{*}_{3}=\textbf{v}_{\alpha,\beta}\cdotp\textbf{g}^{*}_{3}+\textbf{g}^{*}_{\alpha,\beta}\cdotp\textbf{v}_{3}={v}_{k,\alpha\beta}{x}^{*k}_{,3}+{\Gamma}^{*n}_{\alpha\beta}{x}^{*k}_{,n}{v}_{k,3}$
(53)
$\displaystyle={x}^{*k}_{,3}{v}_{k,\alpha\beta}-{\Gamma}^{*\nu}_{\alpha\beta}{x}^{*k}_{,\nu}{x}^{*,\mu}_{k}{v}^{m}_{,\mu}{x}^{*}_{m,3}-{b}^{*}_{\alpha\beta}{x}^{*k}_{,3}{x}^{*,\mu}_{k}{v}^{m}_{,\mu}{x}^{*}_{m,3}={x}^{*k}_{,3}\left({v}_{k,\alpha\beta}-{\Gamma}^{*\mu}_{\alpha\beta}{v}_{k,\mu}\right).$

The detailed derivation of Eqs. (51) and (52) is given in Appendix A.
Equations (51) and (52) define the exact relation of the equidistant and
reference strain rates of a KL shell within the scope of finite, but small,
strain theory. These expressions can be found in [1] where they are derived as
a special case of the general shell theory.

As already noted, many existing IGA shell formulations use a simplified form
of the expression (51). This is justified by the fact that the influence of
the second-order terms with respect to the thickness coordinate can be
neglected for thin shells, [13, 14, 24]. This assumption corresponds to the
following approximation:

${A}^{\mu\nu}_{\alpha\beta}={B}^{\mu\nu}_{\alpha\beta}\approx{\delta}^{\mu\nu}_{\alpha\beta},$
(54)

which reduces Eq. (51) to the linear form:

$\bar{d}_{\alpha\beta}={d}_{\alpha\beta}-\zeta\dot{\kappa}_{\alpha\beta}.$
(55)

A notable exception is the reference [19] where one type of simplified, but
nonlinear, strain distribution is utilized. However, the comparison of
different models is not considered.

It is interesting to note that the physical components of strain are
distributed nonlinearly due to the metric, although the previous equation
imposes a linear strain distribution:

$\bar{d}_{(\alpha\beta)}=\frac{\bar{d}_{\alpha\beta}}{\sqrt{\bar{g}_{\alpha\alpha}\bar{g}_{\beta\beta}}}.$
(56)

The linear analysis performed in [2] shows that for thick shells, the
utilization of a nonlinear strain distribution and an exact shifter is
recommended, along with exact integration. It is also noted that for $Kh<0.1$,
the error introduced by the assumption of a linear strain distribution is not
significant.

Equations (33), (45), (46), (51), and (52) are the basis for the formulation
of the isogeometric element for the finite strain and finite rotation analysis
of an arbitrarily curved KL shell.

##  4 Finite element formulation

In line with the previous derivation, we will formulate isogeometric finite
elements of the KL shell using the principle of virtual power. The present
formulation follows the approach given in [45]. The generalized coordinates
are the components of the velocities of the control points.

Let us start from the generalized Hooke law for the linear elastic material,
also known as the Saint Venant-Kirchhoff material model. This is well-suited
for the small strain and large rotation analysis, [2]. The components of the
stress and strain rates are related as:

$\dot{\sigma}^{i}_{j}=2\mu{d}^{i}_{j}+\lambda{\delta}^{i}_{j}{d}^{k}_{k},$
(57)

where $\mu$ and $\lambda$ are the Lamé material parameters. The plane stress
assumption states that the normal component of the stress rate in the
direction of the shell thickness is negligible compared to the other
components:

$\dot{\sigma}^{3}_{3}=2\mu{d}^{3}_{3}+\lambda\left({d}^{\alpha}_{\alpha}+{d}^{3}_{3}\right)=0,$
(58)

and the corresponding strain rate becomes:

${d}^{3}_{3}=-\frac{\nu}{1-\nu}{d}^{\alpha}_{\alpha}.$ (59)

The existence of the ${d}^{3}_{3}$ component violates the KL hypothesis and it
is neglected here, together with the $\dot{\sigma}^{3}_{3}$. This is a known
inconsistency of the KL theory caused by its approximate nature. Nevertheless,
a great quantity of literature accounts for this component of strain rate
since it increases the accuracy of the mechanical model and is easily
introduced via the static condensation of the material tensor, [13].

Here, the plane stress and the plane strain conditions are enforced and the
stress-strain relation takes its final form:

$\dot{\sigma}^{\alpha\beta}=2\mu\left({g}^{*\alpha\nu}{g}^{*\beta\gamma}+\frac{\nu}{1-\nu}{g}^{*\alpha\beta}{g}^{*\nu\gamma}\right){d}_{\nu\gamma}={D}^{\alpha\beta\nu\gamma}{d}_{\nu\gamma},$
(60)

where ${D}^{\alpha\beta\nu\gamma}$ are the components of the constitutive,
material, tensor.

###  4.1 Principle of virtual power

The principle of virtual power states that at any time, the total power of the
external, internal and inertial forces is zero for any admissible virtual
state of motion. If the inertial effects are neglected and only the surface
loads are considered, it can be written as:

$\delta
P=\int_{v}{\sigma}^{ij}\delta{d}_{ij}\differential{v}-\int_{a}\hat{p}^{i}\delta\hat{v}_{i}\differential{a}=\int_{v}\bm{\sigma}:\delta\textbf{d}\differential{v}-\int_{a}\textbf{p}\delta\textbf{v}\differential{a}=0,$
(61)

where $\bm{\sigma}$ is the Cauchy stress tensor, d is the strain rate tensor,
and p is the vector of external surface loads. All these quantities are
measured with respect to the current configuration. This configuration is
unknown and its determination is the main subject of the nonlinear analysis.
For problems with deformation-independent loads, it is only required to
linearize the Cauchy stress. At the current $(n+1)$ configuration this stress
is approximated by:

$\displaystyle{}^{(n+1)}{\bm{\sigma}}$
$\displaystyle=^{(n)}_{(n+1)}{\bm{\sigma}}+\frac{{\operatorname{d}}^{t}\>{}^{(n)}{\bm{\sigma}}}{\operatorname{d}t}\Delta{t}_{(n+1)},$
(62) $\displaystyle{}^{(n)}_{(n+1)}{\bar{\sigma}}_{\alpha\beta}$
$\displaystyle=\sqrt{\frac{{}^{(n)}{\bar{g}}}{{}^{(n+1)}{\bar{g}}}}\>^{(n)}{\bar{\sigma}}_{\alpha\beta}=\frac{{}^{(n)}{g}_{0}}{{}^{(n+1)}{g}_{0}}\sqrt{\frac{{}^{(n)}{g}}{{}^{(n+1)}{g}}}\>^{(n)}{\bar{\sigma}}_{\alpha\beta},$

where ${}^{(n)}_{(n+1)}{\bar{\sigma}}_{\alpha\beta}$ are the components of the
equidistant stress from the previous $(n)$ configuration expressed with
respect to the metric of the current $(n+1)$ configuration. Similarly,
${}^{(n)}{\bar{\sigma}}_{\alpha\beta}$ is the same quantity but expressed with
respect to the previous $(n)$ configuration. The additive decomposition of
stress, as in Eq. (62), is only valid if all of the terms are expressed with
respect to the same metric. However, the rigorous integration of this
expression is far from trivial due to the ratio of shifters in two different
configurations. If the incremental change of curvature is small, this ratio
can be approximated as ${}^{(n)}{g}_{0}/^{(n+1)}{g}_{0}\approx 1$. An
appropriate time derivative is designated with
$\operatorname{d}^{t}()/\operatorname{d}t$ while
$\Delta{t}_{(n+1)}={t}_{(n+1)}-{t}_{(n)}$ is the time increment.

The rate form of stress-strain relations requires an objective time
derivative. Note that the material derivative is not objective, while the
corotational (Jaumann) and convective time derivatives fulfill this property,
[48, 49]. The relation (62) is rewritten as:

$^{(n+1)}{\bm{\sigma}}=^{(n)}{\bm{\sigma}}+\left({L}_{v}\>^{(n)}{\bm{\sigma}}\right)\Delta{t}_{(n+1)}$
(63)

where ${L}_{v}$ designates Lie time derivative which matches the convective
derivative, for convective coordinates. With respect to the covariant metric
base it is:

${L}_{v}\bm{\sigma}=\dot{\bm{\sigma}}-\textbf{L}\bm{\sigma}-\bm{\sigma}\textbf{L}^{\mathsf{T}},\quad\textbf{L}=\hat{v}_{i|j}\textbf{g}^{i}\otimes\textbf{g}^{j},$
(64)

where L is the spatial velocity gradient. On the other hand, since the
material time derivative of the stress tensor reads:

$\dot{\bm{\sigma}}=\dot{\sigma}^{ij}\textbf{g}_{i}\otimes\textbf{g}_{j}+\textbf{L}\bm{\sigma}+\bm{\sigma}\textbf{L}^{\mathsf{T}},$

relation (64) returns:

${L}_{v}\bm{\sigma}=\dot{\sigma}^{ij}\textbf{g}_{i}\otimes\textbf{g}_{j}\Rightarrow\left({L}_{v}\bm{\sigma}\right)^{ij}=\dot{\sigma}^{ij}.$
(65)

Therefore, the components of the stress rate tensor equal the material
derivative of the components of the stress tensor. This is an important fact
which will be utilized further.

After the insertion of Eqs. (63) and (65) into Eq. (61), the linearized form
of the principle of virtual power at the current configuration is obtained:

$\displaystyle\int_{v}\dot{\bm{\sigma}}:\delta\textbf{d}\differential{v}\Delta
t+\int_{v}\bm{\sigma}:\delta\textbf{d}\differential{v}=\int_{a}\textbf{p}\delta\textbf{v}\differential{a},$
(66)
$\displaystyle\int_{v}\dot{\bar{\sigma}}^{\alpha\beta}\delta\bar{d}_{\alpha\beta}\differential{v}\Delta
t+\int_{v}\bar{\sigma}^{\alpha\beta}\delta\bar{d}_{\alpha\beta}\differential{v}=\int_{a}\hat{p}^{i}\delta\hat{v}_{i}\differential{a}.$

From this equation and further in the text, the notation is slightly
simplified for the sake of readability by removing the time indices and
asterisks. This can be done without introducing any notational ambiguity since
(i) the stress and strain rates are instantaneous quantities, while the known
stress $\bm{\sigma}$ is calculated at the previous configuration and
transformed to the metric of the current configuration, see Eqs. (62) and
(66), and (ii) all integrations are performed with respect to the metric of
the current configuration, in accordance with the updated Lagrangian
procedure, [50].

In order to reduce the dimension from 3D to 2D, it is necessary to integrate
the left-hand side of Eq. (66) along the thickness. Thus, integrals over the
area of midsurface are obtained:

$\displaystyle\int_{v}\dot{\bar{\sigma}}^{\alpha\beta}\delta\bar{d}_{\alpha\beta}\differential{v}\Delta
t+\int_{v}\bar{\sigma}^{\alpha\beta}\delta\bar{d}_{\alpha\beta}\differential{v}=$
$\displaystyle\int_{a}\left(\dot{N}^{\mu\nu}\delta{d}_{\mu\nu}+\dot{M}^{\mu\nu}\delta\dot{\kappa}_{\mu\nu}\right)\differential{a}\Delta
t$ (67)
$\displaystyle+\int_{a}\left({N}^{\mu\nu}\delta{d}_{\mu\nu}+{M}^{\mu\nu}\delta\dot{\kappa}_{\mu\nu}\right)\differential{a},$

where ${N}^{\mu\nu}$ and ${M}^{\mu\nu}$ are stress resultants and stress
couples, which are energetically conjugated with the reference strain rates of
the midsurface, ${d}_{\mu\nu}$ and $\dot{\kappa}_{\mu\nu}$, while
$\dot{N}^{\mu\nu}$ and $\dot{M}^{\mu\nu}$ are their respective rates:

$\displaystyle\dot{N}^{\mu\nu}$
$\displaystyle=\int_{-h/2}^{h/2}\dot{\bar{\sigma}}^{\alpha\beta}{A}^{\mu\nu}_{\alpha\beta}{g}_{0}\differential{\zeta}=\int_{-h/2}^{h/2}\dot{\bar{\sigma}}^{\alpha\beta}\left({\delta}^{\mu}_{\alpha}{\delta}^{\nu}_{\beta}-\zeta^{2}{b}^{\mu}_{\alpha}{b}^{\nu}_{\beta}\right){g}_{0}\differential{\zeta},$
(68) $\displaystyle\dot{M}^{\mu\nu}$
$\displaystyle=-\frac{1}{2}\int_{-h/2}^{h/2}\zeta\dot{\bar{\sigma}}^{\alpha\beta}{B}^{\mu\nu}_{\alpha\beta}{g}_{0}\differential{\zeta}=-\frac{1}{2}\int_{-h/2}^{h/2}\zeta\dot{\bar{\sigma}}^{\alpha\beta}\left({\delta}^{\nu}_{\alpha}\bar{C}^{\mu}_{\beta}+{\delta}^{\mu}_{\beta}\bar{C}^{\nu}_{\alpha}\right){g}_{0}\differential{\zeta}.$

These expressions are the same as those obtained in [51] but they are rarely
applied for the nonlinear analysis due to the complexity. If we introduce the
vectors of generalized section forces and strain rates of the shell
midsurface:

$\displaystyle\textbf{f}^{\mathsf{T}}$
$\displaystyle=\begin{bmatrix}{N}^{11}&{N}^{22}&{N}^{12}&{M}^{11}&{M}^{22}&{M}^{12}\end{bmatrix},$
(69) $\displaystyle\textbf{e}^{\mathsf{T}}$
$\displaystyle=\begin{bmatrix}{d}_{11}&{d}_{22}&2{d}_{12}&\dot{\kappa}_{11}&\dot{\kappa}_{22}&2\dot{\kappa}_{12}\end{bmatrix},$

Eq. (66) can be written in compact matrix form as:

$\int_{a}\dot{\textbf{f}}^{\mathsf{T}}\delta\textbf{e}\differential{a}\Delta
t+\int_{a}\textbf{f}^{\mathsf{T}}\delta\textbf{e}\differential{a}=\int_{a}\textbf{p}\delta\textbf{v}\differential{a}.$
(70)

The geometrically exact relations (67) and (68) are crucial for the accurate
formulation of structural shell theories. In particular, energetically
conjugated pairs are defined rigorously and the appropriate constitutive
matrices are guaranteed to be symmetric, as will be shown in the next section.

###  4.2 Analytical relation between energetically conjugated pairs

The integration represented with Eqs. (67) and (68) is an essential step for
the development of the structural theory of shells, since it allows the
calculation of virtual power by the reference quantities of the midsurface.
This requires the constitutive relation between the stress and strain rates at
an arbitrary point of the shell continuum. Analogously to Eq. (60), this
relation is:

$\dot{\bar{\sigma}}^{\alpha\beta}=2\mu\left(\bar{g}^{\alpha\nu}\bar{g}^{\beta\gamma}+\frac{\nu}{1-\nu}\bar{g}^{\alpha\beta}\bar{g}^{\nu\gamma}\right)\bar{d}_{\nu\gamma}=\bar{D}^{\alpha\beta\nu\gamma}\bar{d}_{\nu\gamma}.$
(71)

Using Eq. (19), it is evident that the following transformation is valid:

$\displaystyle\bar{D}^{\mu\nu\epsilon\varphi}$
$\displaystyle=2\mu\left({C}^{\mu}_{\rho}{C}^{\epsilon}_{\delta}{g}^{\rho\delta}{C}^{\nu}_{\alpha}{C}^{\varphi}_{\beta}{g}^{\alpha\beta}+\frac{\nu}{1-\nu}{C}^{\mu}_{\rho}{C}^{\nu}_{\alpha}{g}^{\rho\alpha}{C}^{\epsilon}_{\delta}{C}^{\varphi}_{\beta}{g}^{\delta\beta}\right)$
(72)
$\displaystyle={C}^{\mu}_{\rho}{C}^{\nu}_{\alpha}{C}^{\epsilon}_{\delta}{C}^{\varphi}_{\beta}2\mu\left({g}^{\rho\delta}{g}^{\alpha\beta}+\frac{\nu}{1-\nu}{g}^{\rho\alpha}{g}^{\delta\beta}\right)={C}^{\mu}_{\rho}{C}^{\nu}_{\alpha}{C}^{\epsilon}_{\delta}{C}^{\varphi}_{\beta}{D}^{\rho\alpha\delta\beta}.$

This relation explicitly relates the components of constitutive tensor at an
arbitrary point with the ones at the midsurface via the reciprocal shift
tensor. In this way, the constitutive relations between the energetically
conjugated rates of generalized section forces and reference strain rates of
the midsurface are:

$\displaystyle\dot{N}^{\alpha\beta}$
$\displaystyle={D}^{\alpha\beta\gamma\lambda}_{M}{d}_{\gamma\lambda}+{D}^{\alpha\beta\gamma\lambda}_{MB}\dot{\kappa}_{\gamma\lambda},$
(73) $\displaystyle\dot{M}^{\alpha\beta}$
$\displaystyle={D}^{\alpha\beta\gamma\lambda}_{BM}{d}_{\gamma\lambda}+{D}^{\alpha\beta\gamma\lambda}_{B}\dot{\kappa}_{\gamma\lambda},$

where:

$\displaystyle{D}^{\alpha\beta\gamma\lambda}_{M}$
$\displaystyle=\int_{-h/2}^{h/2}{A}^{\alpha\beta}_{\mu\nu}\bar{D}^{\mu\nu\delta\varphi}{A}^{\gamma\lambda}_{\delta\varphi}{g}_{0}\differential{\zeta},$
(74) $\displaystyle{D}^{\alpha\beta\gamma\lambda}_{MB}$
$\displaystyle=-\int_{-h/2}^{h/2}\zeta{A}^{\alpha\beta}_{\mu\nu}\bar{D}^{\mu\nu\delta\varphi}{B}^{\gamma\lambda}_{\delta\varphi}{g}_{0}\differential{\zeta},$
$\displaystyle{D}^{\alpha\beta\gamma\lambda}_{BM}$
$\displaystyle=-\int_{-h/2}^{h/2}\zeta{B}^{\alpha\beta}_{\mu\nu}\bar{D}^{\mu\nu\delta\varphi}{A}^{\gamma\lambda}_{\delta\varphi}{g}_{0}\differential{\zeta}={D}^{\gamma\lambda\alpha\beta}_{MB},$
$\displaystyle{D}^{\alpha\beta\gamma\lambda}_{B}$
$\displaystyle=\int_{-h/2}^{h/2}\zeta^{2}{B}^{\alpha\beta}_{\mu\nu}\bar{D}^{\mu\nu\delta\varphi}{B}^{\gamma\lambda}_{\delta\varphi}{g}_{0}\differential{\zeta}.$

Now, we can represent Eq. (73) in a compact matrix form as:

$\dot{\textbf{f}}=\textbf{D}\textbf{e},$ (75)

where D is the symmetric constitutive tensor at the current configuration:

$\textbf{D}=\begin{bmatrix}\textbf{D}_{M}&\textbf{D}_{MB}\\\
\textbf{$\;\textbf{D}_{MB}$}^{\mathsf{T}}&\textbf{D}_{B}\end{bmatrix}=\begin{bmatrix}{D}^{\alpha\beta\gamma\lambda}_{M}&{D}^{\alpha\beta\gamma\lambda}_{MB}\\\
{D}^{\gamma\lambda\alpha\beta}_{MB}&{D}^{\alpha\beta\gamma\lambda}_{B}\end{bmatrix}\textbf{g}_{\alpha}\otimes\textbf{g}_{\beta}\otimes\textbf{g}_{\gamma}\otimes\textbf{g}_{\lambda}.$
(76)

Complete expressions for the components of the material tensor D are given in
Appendix B along with the closed-form solutions of the appropriate integrals.
To the best of the authors’ knowledge, this is the first complete analytical
representation of the constitutive relation between the energetically
conjugated pairs for the elastic KL shell.

###  4.3 Variation of strains

Since the strain rate is a function of the generalized coordinates as well as
the metric, we must vary it with respect to both arguments. By noting that the
variation of the tangent vectors can be expressed as:

$\delta\textbf{g}_{\alpha}=\delta\textbf{v}_{\alpha}\Delta t,$ (77)

the variations of the reference strains are given by:

$\displaystyle\delta{d}_{\alpha\beta}$
$\displaystyle=\delta\frac{1}{2}\left(\textbf{g}_{\alpha}\cdotp\textbf{v}_{\beta}+\textbf{g}_{\beta}\cdotp\textbf{v}_{\alpha}\right)=\frac{1}{2}\left(\delta\textbf{v}_{\alpha}\cdotp\textbf{v}_{\beta}+\delta\textbf{v}_{\beta}\cdotp\textbf{v}_{\alpha}\right)\Delta
t+\frac{1}{2}\left(\textbf{g}_{\alpha}\cdotp\delta\textbf{v}_{\beta}+\textbf{g}_{\beta}\cdotp\delta\textbf{v}_{\alpha}\right),$
(78) $\displaystyle\delta\dot{\kappa}_{\alpha\beta}$
$\displaystyle=\delta\left[\textbf{g}_{3}\cdotp\left(\textbf{v}_{\alpha,\beta}-{\Gamma}^{\mu}_{\alpha\beta}\textbf{v}_{\mu}\right)\right]$
$\displaystyle=\left[\delta\textbf{v}_{3}\cdotp\left(\textbf{v}_{\alpha,\beta}-{\Gamma}^{\lambda}_{\alpha\beta}\textbf{v}_{\lambda}\right)-\delta{\Gamma}^{\lambda}_{\alpha\beta}\left(\textbf{g}_{3}\cdotp\textbf{v}_{\lambda}\right)\right]\Delta
t+\textbf{g}_{3}\cdotp\left(\delta\textbf{v}_{\alpha.\beta}-{\Gamma}^{\mu}_{\alpha\beta}\delta\textbf{v}_{\mu}\right).$

Most of the terms in the previous expression can be easily computed since the
variation is explicitly performed with respect to the unknown variables.
However, the first addend of $\dot{\kappa}_{\alpha\beta}$ is an exception
since it requires the variation of the normal and the Christoffel symbols.
These variations must be represented via the variations of generalized
coordinates.

The variation of the normal can be performed analogously as the material
differentiation given in Eq. (42):

$\delta\textbf{v}_{3}=-\left(\textbf{g}_{3}\cdotp\delta\textbf{v}_{\alpha}\right)\textbf{g}^{\alpha}.$
(79)

For the variation of the Christoffel symbols we need the variations of
reciprocal base vectors (6):

$\displaystyle\textbf{g}^{1}$
$\displaystyle=\left(\textbf{g}_{2}\times\textbf{g}_{3}\right)/\sqrt{g}\Rightarrow\delta\textbf{g}^{1}=-\left[\left(\delta\textbf{v}_{2}\cdotp\textbf{g}^{1}\right)\textbf{g}^{2}+\left(\delta\textbf{v}_{3}\cdotp\textbf{g}^{1}\right)\textbf{g}^{3}+\left(\delta\textbf{v}_{1}\cdotp\textbf{g}^{1}\right)\textbf{g}^{1}\right]\Delta
t,$ (80) $\displaystyle\textbf{g}^{2}$
$\displaystyle=\left(\textbf{g}_{3}\times\textbf{g}_{1}\right)/\sqrt{g}\Rightarrow\delta\textbf{g}^{2}=-\left[\left(\delta\textbf{v}_{3}\cdotp\textbf{g}^{2}\right)\textbf{g}^{3}+\left(\delta\textbf{v}_{1}\cdotp\textbf{g}^{2}\right)\textbf{g}^{1}+\left(\delta\textbf{v}_{2}\cdotp\textbf{g}^{2}\right)\textbf{g}^{2}\right]\Delta
t,$

which allows us to write:

$\delta{\Gamma}^{\lambda}_{\alpha\beta}=\delta\left(\textbf{g}_{\alpha,\beta}\cdotp\textbf{g}^{\lambda}\right)=\delta\textbf{g}_{\alpha,\beta}\cdotp\textbf{g}^{\lambda}+\textbf{g}_{\alpha,\beta}\cdotp\delta\textbf{g}^{\lambda}=\left[\textbf{g}^{\lambda}\cdotp\delta\textbf{v}_{\alpha,\beta}-{\Gamma}^{\nu}_{\alpha\beta}\left(\textbf{g}^{\lambda}\cdot\delta\textbf{v}_{\nu}\right)\right]\Delta
t.$ (81)

By inserting Eqs. (79) and (81) into Eq. (78), we obtain the final expression
for the variation of the rate of the change of curvature:

$\displaystyle\delta\dot{\kappa}_{\alpha\beta}$
$\displaystyle=\textbf{g}_{3}\cdotp\left(\delta\textbf{v}_{\alpha,\beta}-{\Gamma}^{\mu}_{\alpha\beta}\delta\textbf{v}_{\mu}\right)$
(82)
$\displaystyle\;\;\;\;-\left\\{\left(\textbf{g}_{3}\cdotp\delta\textbf{v}_{\mu}\right)\textbf{g}^{\mu}\cdotp\left(\textbf{v}_{\alpha,\beta}-{\Gamma}^{\nu}_{\alpha\beta}\textbf{v}_{\nu}\right)+\left[\textbf{g}^{\nu}\cdotp\delta\textbf{v}_{\alpha,\beta}-{\Gamma}^{\mu}_{\alpha\beta}\left(\textbf{g}^{\nu}\cdotp\delta\textbf{v}_{\mu}\right)\right]\left(\textbf{g}_{3}\cdotp\textbf{v}_{\nu}\right)\right\\}\Delta
t.$

###  4.4 Discrete equation of motion

Let us find the matrix $\textbf{B}_{L}$ which relates the reference strain
rates of the midsurface with the velocities of control points:

$\textbf{e}=\textbf{B}_{L}\dot{\textbf{q}}.$ (83)

The vector of reference strain rates, using Eqs. (33) and (53), can be
represented as:

$\begin{gathered}\textbf{e}=\textbf{H}\textbf{w}=\sum_{m=1}^{3}\textbf{H}_{m}\textbf{w}^{m},\quad\textbf{H}=\begin{bmatrix}\textbf{H}_{1}&\textbf{H}_{2}&\textbf{H}_{3}\end{bmatrix},\\\
\textbf{H}_{m}=\begin{bmatrix}{x}_{m,1}&0&0&0&0\\\ 0&{x}_{m,2}&0&0&0\\\
{x}_{m,2}&{x}_{m,1}&0&0&0\\\
-{\Gamma}^{1}_{11}{x}_{m,3}&-{\Gamma}^{2}_{11}{x}_{m,3}&{x}_{m,3}&0&0\\\
-{\Gamma}^{1}_{22}{x}_{m,3}&-{\Gamma}^{2}_{22}{x}_{m,3}&0&{x}_{m,3}&0\\\
-2{\Gamma}^{1}_{12}{x}_{m,3}&-2{\Gamma}^{2}_{12}{x}_{m,3}&0&0&2{x}_{m,3}\end{bmatrix},\quad\textbf{w}^{m}=\begin{bmatrix}{v}^{m}_{,1}\\\
{v}^{m}_{,2}\\\ {v}^{m}_{,11}\\\ {v}^{m}_{,22}\\\ {v}^{m}_{,12}\\\
\end{bmatrix},\end{gathered}$ (84)

where the vector w is:

$\textbf{w}=\textbf{B}\dot{\textbf{q}},\;\textbf{w}=\begin{bmatrix}\textbf{w}^{1}\\\
\textbf{w}^{2}\\\
\textbf{w}^{3}\end{bmatrix},\quad\textbf{B}=\begin{bmatrix}\textbf{B}^{1}_{11}&...&\textbf{B}^{1}_{IJ}&...&\textbf{B}^{1}_{NM}\\\
\textbf{B}^{2}_{11}&...&\textbf{B}^{2}_{IJ}&...&\textbf{B}^{2}_{NM}\\\
\textbf{B}^{3}_{11}&...&\textbf{B}^{3}_{IJ}&...&\textbf{B}^{3}_{NM}\end{bmatrix}.$
(85)

The submatrices $\textbf{B}^{m}_{IJ}$ for an arbitrary control point $IJ$
consists of the derivatives of basis functions:

$\textbf{B}^{1}_{IJ}=\begin{bmatrix}{R}_{IJ,1}&0&0\\\ {R}_{IJ,2}&0&0\\\
{R}_{IJ,11}&0&0\\\ {R}_{IJ,22}&0&0\\\ {R}_{IJ,12}&0&0\\\
\end{bmatrix},\quad\textbf{B}^{2}_{IJ}=\begin{bmatrix}0&{R}_{IJ,1}&0\\\
0&{R}_{IJ,2}&0\\\ 0&{R}_{IJ,11}&0\\\ 0&{R}_{IJ,22}&0\\\ 0&{R}_{IJ,12}&0\\\
\end{bmatrix},\quad\textbf{B}^{3}_{IJ}=\begin{bmatrix}0&0&{R}_{IJ,1}\\\
0&0&{R}_{IJ,2}\\\ 0&0&{R}_{IJ,11}\\\ 0&0&{R}_{IJ,22}\\\ 0&0&{R}_{IJ,12}\\\
\end{bmatrix}.$ (86)

The matrix $\textbf{B}_{L}$ now follows as:

$\textbf{e}=\textbf{H}\textbf{w}=\textbf{H}\textbf{B}\dot{\textbf{q}}=\textbf{B}_{L}\dot{\textbf{q}},\quad\textbf{B}_{L}=\textbf{H}\textbf{B}.$
(87)

Furthermore the requirement, in order to rigorously define the geometric
stiffness term, is to find the explicit matrix form of the virtual power
generated by the known stress and the variation of strain rate with respect to
the metric, Eqs. (70) and (78):

$\displaystyle\int_{a}\textbf{f}^{\mathsf{T}}\delta\textbf{B}_{L}\dot{\textbf{q}}\differential{a}\Delta
t$
$\displaystyle=\int_{a}\sum_{m=1}^{3}\sum_{n=1}^{3}\textbf{$(\textbf{w}^{m})$}^{\mathsf{T}}\textbf{G}^{n}_{m}\delta\textbf{w}_{m}\differential{a}\Delta
t$ (88)
$\displaystyle=\frac{1}{2}\int_{a}{N}^{\alpha\beta}\left({v}^{m}_{,\beta}\delta{v}_{m,\alpha}+{v}^{m}_{,\alpha}\delta{v}_{m,\beta}\right)\differential{a}\Delta
t$
$\displaystyle\;\;\;\;+\int_{a}{M}^{\alpha\beta}\left[{v}^{k}_{,\nu}\left({\Gamma}^{\nu}_{\alpha\beta}{x}^{n}_{,3}{x}^{,\mu}_{k}+{\Gamma}^{\mu}_{\alpha\beta}{x}^{n,\nu}{x}_{k,3}\right)\delta{v}_{n,\mu}\right.$
$\displaystyle\left.\;\;\;\;-{v}^{k}_{,\alpha\beta}{x}^{n}_{,3}{x}^{,\mu}_{k}\delta{v}_{n,\mu}-{v}^{k}_{,\nu}{x}^{n,\nu}{x}_{k,3}\delta{v}_{n,\alpha\beta}\right]\differential{a}\Delta
t,$

where the matrix $\textbf{G}^{n}_{m}$ is:

$\textbf{G}^{n}_{m}=\begin{bmatrix}\left({G}^{n}_{m}\right)_{11}&\left({G}^{n}_{m}\right)_{12}&-{M}^{11}{x}^{n,1}{x}_{m,3}&-{M}^{22}{x}^{n,1}{x}_{m,3}&-2{M}^{12}{x}^{n,1}{x}_{m,3}\\\
\left({G}^{n}_{m}\right)_{21}&\left({G}^{n}_{m}\right)_{22}&-{M}^{11}{x}^{n,2}{x}_{m,3}&-{M}^{22}{x}^{n,2}{x}_{m,3}&-2{M}^{12}{x}^{n,2}{x}_{m,3}\\\
-{M}^{11}{x}^{n}_{,3}{x}^{,1}_{m}&-{M}^{11}{x}^{n}_{,3}{x}^{,2}_{m}&0&0&0\\\
-{M}^{22}{x}^{n}_{,3}{x}^{,1}_{m}&-{M}^{22}{x}^{n}_{,3}{x}^{,2}_{m}&0&0&0\\\
-2{M}^{12}{x}^{n}_{,3}{x}^{,1}_{m}&-2{M}^{12}{x}^{n}_{,3}{x}^{,2}_{m}&0&0&0\end{bmatrix},$
(89)

with:

$\displaystyle\left({G}^{n}_{m}\right)_{11}$
$\displaystyle={\Gamma}^{1}_{M}\left({x}^{n}_{,3}{x}^{,1}_{m}+{x}^{n,1}{x}_{m,3}\right)+{\delta}^{n}_{m}{N}^{11},$
(90) $\displaystyle\left({G}^{n}_{m}\right)_{12}$
$\displaystyle={\Gamma}^{1}_{M}{x}^{n}_{,3}{x}^{,2}_{m}+{\Gamma}^{2}_{M}{x}^{n,1}{x}_{m,3}+{\delta}^{n}_{m}{N}^{12},$
$\displaystyle\left({G}^{n}_{m}\right)_{21}$
$\displaystyle={\Gamma}^{2}_{M}{x}^{n}_{,3}{x}^{,1}_{m}+{\Gamma}^{1}_{M}{x}^{n,2}{x}_{m,3}+{\delta}^{n}_{m}{N}^{21},$
$\displaystyle\left({G}^{n}_{m}\right)_{22}$
$\displaystyle={\Gamma}^{2}_{M}\left({x}^{n}_{,3}{x}^{,2}_{m}+{x}^{n,2}{x}_{m,3}\right)+{\delta}^{n}_{m}{N}^{22},$
$\displaystyle{\Gamma}^{\nu}_{M}$
$\displaystyle={\Gamma}^{\nu}_{11}{M}^{11}+2{\Gamma}^{\nu}_{12}{M}^{12}+{\Gamma}^{\nu}_{22}{M}^{22}.$

By introducing the total matrix of the generalized section forces:

$\textbf{G}=\begin{bmatrix}\textbf{G}^{1}_{1}&\textbf{G}^{2}_{1}&\textbf{G}^{3}_{1}\\\
\textbf{G}^{1}_{2}&\textbf{G}^{2}_{2}&\textbf{G}^{3}_{2}\\\
\textbf{G}^{1}_{3}&\textbf{G}^{2}_{3}&\textbf{G}^{3}_{3}\end{bmatrix},$ (91)

expression (88) can be written as:

$\int_{a}\textbf{f}^{\mathsf{T}}\delta\textbf{B}_{L}\dot{\textbf{q}}\differential{a}\Delta
t=\int_{a}\textbf{w}^{\mathsf{T}}\textbf{G}\delta\textbf{w}\differential{a}\Delta
t.$ (92)

A careful inspection of Eqs. (89) and (90) reveals the fact that the matrix of
generalized section forces, G, is symmetric. This form was obtained in [45]
first. Now, the terms in the equation of the virtual power, (70), reduce to:

$\displaystyle\int_{a}\dot{\textbf{f}}^{\mathsf{T}}\delta\textbf{e}\differential{a}$
$\displaystyle\approx\int_{a}\dot{\textbf{q}}^{\mathsf{T}}\textbf{$\textbf{B}_{L}$}^{\mathsf{T}}\textbf{D}\textbf{B}_{L}\delta\dot{\textbf{q}}\differential{a},\quad\left(\delta\textbf{e}\approx\frac{1}{2}\left(\textbf{g}_{\alpha}\cdotp\delta\textbf{v}_{\beta}+\textbf{g}_{\beta}\cdotp\delta\textbf{v}_{\alpha}\right)=\textbf{B}_{L}\delta\dot{\textbf{q}}\right)$
(93)
$\displaystyle\int_{a}\textbf{f}^{\mathsf{T}}\delta\textbf{e}\differential{a}$
$\displaystyle=\int_{a}\textbf{f}^{\mathsf{T}}\left(\delta\textbf{B}_{L}\dot{\textbf{q}}+\textbf{B}_{L}\delta\dot{\textbf{q}}\right)\differential{a}=\int_{a}\left(\dot{\textbf{q}}^{\mathsf{T}}\textbf{B}^{\mathsf{T}}\textbf{G}\textbf{B}\delta\dot{\textbf{q}}\Delta
t+\textbf{f}^{\mathsf{T}}\textbf{B}_{L}\delta\dot{\textbf{q}}\right)\differential{a},$

and the equation of equilibrium becomes:

$\dot{\textbf{q}}^{\mathsf{T}}\int_{a}\left(\textbf{$\textbf{B}_{L}$}^{\mathsf{T}}\textbf{D}\textbf{B}_{L}+\textbf{B}^{\mathsf{T}}\textbf{G}\textbf{B}\right)\differential{a}\delta\dot{\textbf{q}}\Delta
t=\int_{a}\textbf{p}^{\mathsf{T}}\textbf{N}\differential{a}\delta\dot{\textbf{q}}-\int_{a}\textbf{f}^{\mathsf{T}}\textbf{B}_{L}\differential{a}\delta\dot{\textbf{q}}.$
(94)

Eq. (94) reduces to the standard form:

$\textbf{K}_{T}\Delta\textbf{q}=\textbf{Q}-\textbf{F},\quad\left(\Delta\textbf{q}=\dot{\textbf{q}}\Delta
t\right),$ (95)

where:

$\textbf{K}_{T}=\int_{a}\textbf{$\textbf{B}_{L}$}^{\mathsf{T}}\textbf{D}\textbf{B}_{L}\differential{a}+\int_{a}\textbf{B}^{\mathsf{T}}\textbf{G}\textbf{B}\differential{a},$
(96)

is the tangent stiffness matrix and:

$\textbf{Q}=\int_{a}\textbf{N}^{\mathsf{T}}\textbf{p}\differential{a},\quad\textbf{F}=\int_{a}\textbf{$\textbf{B}_{L}$}^{\mathsf{T}}\textbf{f}\differential{a},$
(97)

are the vectors of the external and internal forces, respectively. The vector
$\Delta\textbf{q}$ in Eq. (95) contains increments of displacements of control
points with respect to the previous configuration. Due to the approximations
introduced with Eq. (93), the solution of Eq. (95) does not satisfy the
principle of virtual power and some additional numerical procedure is
required, see Section 4.4.1.

Although the derivation of the geometric stiffness matrix represents a
standard procedure for nonlinear shell formulations, [13], the present
derivation differs. In particular, the full KL shell metric is incorporated in
Eq. (88) and the symmetric nature of the geometric stiffness term is stressed
by the elegant and compact form of Eqs. (91) and (96). This confirms that the
formulation is correct and that the adopted force and strain quantities are
energetically conjugated.

####  4.4.1 Numerical solution procedure

The incremental form of Eq. (95) emphasizes that the equilibrium path must be
considered as a set of discrete equilibrium points
$(\textbf{q}^{i},\textbf{Q}^{i})$ which correspond to a set of fictitious
moments in time. In order to find these points, some iterative procedure must
be employed. Here, the arc-length method is primarily used for all examples.
It is worth noting that the Newton-Raphson method is applicable for the
examples that do not exhibit snap-through behavior. We compared the
performance of the linearized arc-length method, the cylindrical arc-length
method, and the modified Riks method, [52, 53]. For the numerical tests which
follow in the next section, all of these variants return similar results. The
exception is the last example where only the modified Riks method, as
implemented in Abaqus, found the complete structural response, [54]. Hence,
this procedure is briefly presented.

The main idea of the arc-length method is to introduce a load proportionality
factor $\lambda^{i}$ which scales the value of total external load,
$\textbf{Q}^{f}$, and allows it to increase or decrease during the deformation
process:

$\textbf{Q}^{i}={\lambda}^{i}\textbf{Q}^{f}.$ (98)

This approach results in an indeterminate system of $3NM$ equations with
$3NM+1$ unknowns where the additional unknown is the load proportionality
factor itself, [52]. We start from some known, converged, configuration
defined in a point $(\textbf{q}^{i},{\lambda}^{i})$ and look for the next
equilibrium point $(\textbf{q}^{i+1},{\lambda}^{i+1})$. The first step is to
find a predictor solution. Therefore, the predictor tangential displacement
vector, ${}^{0}{\textbf{q}}_{T}$, is calculated from the known tangent
stiffness and the total external load:

${\textbf{K}}^{i}_{T}\>{}^{0}{\textbf{q}}^{i+1}_{T}=\textbf{Q}^{f}.$ (99)

The predictor load increment $\Delta\>^{0}{\lambda}^{i+1}$ follows from the
appropriate constraint equation, [52, 54]:

$\Delta\>^{0}{\lambda}^{i+1}=\pm\frac{\Delta\>{l}^{i+1}}{\sqrt{{}^{0}{\tilde{\textbf{q}}}^{i+1}_{T}\>\cdotp\>^{0}{\tilde{\textbf{q}}}^{i+1}_{T}+1}}.$
(100)

Here, $\Delta{l}^{i+1}$ is the arc-length value which is estimated for the
first increment, and calculated for all the other increments as a function of
desired number of iterations. For convenience, the displacement vectors are
scaled by the maximum displacement component of the corresponding linear
solution and designated with tilde, [54]. The sign of the predictor solution
is calculated from the condition that the projection of the predictor
tangential displacement onto the previously converged displacement increment
must be positive, ensuring continuation along the equilibrium path. This
condition can be written as:

$\Delta\>^{0}{\lambda}^{i+1}\left({}^{0}{\tilde{\textbf{q}}}^{i+1}_{T}\Delta\>{\tilde{\textbf{q}}}^{i}+\Delta\>{\lambda}^{i}\right)>0.$
(101)

In general, the predictor solution does not satisfy the equilibrium and the
arc-length iterations act as correctors. For the first iteration, the
increments of displacement and load are initialized from the predictor
solutions:

$\Delta\>^{1}{\lambda}^{i+1}=\Delta\>^{0}{\lambda}^{i+1},\;\Delta\>^{1}{\textbf{q}}^{i+1}=\Delta\>^{0}{\lambda}^{i+1}\>{}^{0}{\textbf{q}}^{i+1}_{T}.$
(102)

Now, the current tangent stiffness, ${}^{j}{\textbf{K}}^{i+1}_{T}$, and
internal forces, ${}^{j}{\textbf{F}}^{i+1}$, are calculated and convergence is
checked for each iteration $j=1,2,3...$. The convergence criteria can be
defined with respect to the current unbalanced forces,
${}^{j}{\bm{\Psi}}^{i+1}=\textbf{Q}^{i}-\>^{j}{\textbf{F}}^{i+1}$, or/and with
respect to the previously calculated iterative displacement. If the defined
criteria is not met, new tangential, ${}^{j}{\textbf{q}}^{i+1}_{T}$, and
residual, $\delta\>^{j}{\textbf{q}}^{i+1}$, displacements are calculated from:

$^{j}{\textbf{K}}^{i+1}_{T}\>{}^{j}{\textbf{q}}^{i+1}_{T}={\textbf{Q}}^{f},\;\;^{j}{\textbf{K}}^{i+1}_{T}\>\delta\>^{j}{\textbf{q}}^{i+1}=\>^{j}{\bm{\Psi}}^{i+1}.$
(103)

The iterative load factor, $\delta\lambda$, follows from the condition that a
new potential equilibrium point is orthogonal to the tangential displacement
in the arc-length solution space:

$\delta\lambda=-\frac{\delta\>^{j}{\tilde{\textbf{q}}}^{i+1}\>\cdotp\>^{0}{\tilde{\textbf{q}}}^{i+1}_{T}}{{}^{j}{\tilde{\textbf{q}}}^{i+1}_{T}\>\cdotp\>^{0}{\tilde{\textbf{q}}}^{i+1}_{T}+1}.$
(104)

Finally, the new solution is:

$\displaystyle{}^{j+1}{\textbf{q}}^{i+1}$
$\displaystyle={\textbf{q}}^{i}+\Delta\>^{j}{\textbf{q}}^{i+1}+\delta\>^{j}{\textbf{q}}^{i+1}+\delta\lambda\>^{j}{\textbf{q}}^{i+1}_{T},$
(105) $\displaystyle{}^{j+1}{\lambda}^{i+1}$
$\displaystyle={\lambda}^{i}+\Delta\>^{j}{\lambda}^{i+1}+\delta\lambda,$

and the procedure is repeated until the convergence criteria is fulfilled,
[54].

##  5 Numerical examples

There are four main objectives of the present numerical analysis: i)
verification and validation of the developed nonlinear formulation and its
implementation, ii) examination of the influence of full constitutive relation
on the structural response, iii) comparison of the linear and nonlinear
relation between reference and equidistant strains, iv) assessment of the
influence of the curviness.

To examine the influence of the different constitutive relations, we introduce
four models $D$:

  * •

${D}^{a}$ is the full analytical constitutive model. It is given in an
expanded form in Appendix B.

  * •

${D}^{0}$ decouples the bending and membrane actions (e.g., [11]).

  * •

${D}^{1}$ employs a membrane-bending coupling based on the first-order Taylor
approximation of the determinant of the reciprocal shift tensor $1/{g}_{0}$ as
suggested in [47], (C3).

  * •

${D}^{2}$ is like ${D}^{1}$ but introduces an additional effect due to
curvature (similar to the Flugge-Lure-Byrne-Truesdell equations, [1]), (C3).

In order to examine the influence of the constitutive relation on the
structural response, reference strains are observed and compared for different
models. For ${D}^{a}$ model, the influence of the relation between the
reference and equidistant strains is also examined by evaluating the strains
on the outer surface of shell using the exact (51) and the linearized relation
(55). Furthermore, the stresses are calculated by Eq. (71). Since the
equilibrium paths of stresses and strains for nonlinear shells are rarely
found in the literature, these quantities are here compared with the results
from Abaqus simulations, which utilize the standard assumption of the linear
strain distribution, [54]. Since the time is a fictitious quantity in the
present static analysis, strain and stress rates are equal to strains and
stresses, respectively.

We have analyzed each nonlinear example with various different meshes. For the
sake of brevity, only a few results are reported here. Moreover, the solutions
have singularities due to point forces and no optimal higher-order convergence
can be achieved.

Due to the ${C}^{1}$ interelement continuity, reference membrane strains and
${\kappa}_{12}$ are continuous at the interelement boundaries, Eqs. (33) and
(53). However, the curvature changes ${\kappa}_{\alpha\alpha}$ require at
least ${C}^{2}$ continuous meshes and they are linearly averaged here for
${C}^{1}$ meshes. Generally, the post-processing of results in IGA is an
ongoing topic, [55]. The models in Abaqus are meshed with S4R shell elements
and element nodal results are again linearly averaged. This general-purpose
four-node quadrilateral element with reduced integration is readily utilized
for the benchmarks of shell formulations, [56]. First, a classical linear
shell example is investigated. The other examples address nonlinear responses
of the different shell constitutive models.

The boundary conditions are imposed in a well-known manner where rotations are
treated with special care since they are not utilized as DOFs, [11]. The
symmetry boundary conditions are enforced by constraining the displacements of
control points that influences rotation. Numerical integration is performed on
an element level using the Gauss quadrature with $p\times p$ integration
points. There are more efficient quadrature schemes for IGA [57], but they are
not considered in this paper. Furthermore, locking effects are not studied
here.

All the results are presented with respect to the load proportionality factor
(LPF), rather than the load intensity itself.

###  5.1 Linear analysis of pinched cylinder

The first example deals with the linear analysis of a pinched cylinder with
rigid diaphragms at its ends as shown in Fig. 2a. The deflection at point $A$
where the force is applied commonly serves as a benchmark test for shell
elements. Since the structure has three planes of symmetry, only one eighth of
the shell is considered. The analytical solution based on the Fourier series
approximation exists, but it involves some numerical issues, as discussed in
[12]. The reference value often used in the literature is $1.82488\times
10^{-5}$ and it stems from the approximation of the solution with $80\times
80$ Fourier terms. Here, the reference value of $1.827158\times 10^{-5}$
obtained with $8192\times 8192$ Fourier terms is adopted, [19]. The analysis
results for different IGA meshes are shown in Fig. 2b where the expected
convergence behavior can be observed.

Furthermore, the influence of the different constitutive models on the
deflection at point $A$ is examined. The midsurface is represented by
$36\times 36$ cubic $C^{1}$ elements and the curviness $Kh$ is varied. Since
the geometry of midsurface is kept constant, the curviness becomes a function
of thickness only, $Kh=h/300$. In Fig. 3, the obtained values of displacement
component $w_{A}$ for the three reduced constitutive models are related to the
one computed by the full analytical $D^{a}$ model. As expected, these relative
differences increase with the curviness and the results of the $D^{2}$ model
are the closest to the $D^{a}$ model. Note that for strongly curved shells,
the $D^{0}$ model is more accurate than the $D^{1}$ model, which implies that
the simplified coupling of membrane and bending actions introduced in the
$D^{1}$ model does not improve the accuracy in this example.

Figure 2: Linear analysis of the pinched cylinder. a) Geometry, material
characteristics, and applied load. b) Convergence of the radial displacement
at point A for different NURBS orders and interelement continuities. Figure 3:
The relative difference of displacement for different constitutive material
models with respect to the full analytical one.

###  5.2 Snap behavior of a shallow shell

The second problem deals with the shallow shell illustrated in Fig. 4a. This
is a standard example which exhibits both snap-through and snap-back behavior.
It serves as a rigorous test for algorithms that solve complex nonlinear
behavior of shells. Several variants of this example exist, [52], and we are
referring to the ones discussed in [56]. By enforcing the symmetry boundary
conditions, only one fourth of the structure is analyzed. A mesh with $4\times
4$ quartic elements with $C^{3}$ interelement continuity is employed and the
results are compared with those from [56] in Fig. 4b. Almost full agreement is
observed for the two different values of thickness. The results for all four
constitutive models are practically the same due to the small curviness of
this shell at all configurations. Overall, these results confirm that the
presented formulation can deal with complex equilibrium paths, but it should
be noted that these shallow shells do not undergo large rotations.

The following efficiency comparison of the different constitutive models is
performed on the shallow shell with thickness $t=12.7$. Three series of
calculations are done for each constitutive model using the Intel Core i7-8700
CPU. The results are linearly averaged and shown in Table 1. In each test, the
prescribed accuracy of a solution is reached after 21 load increments and a
total of 52 arc-length iterations. Evidently, the full analytical constitutive
model $D^{a}$ requires significantly more computational time than the reduced
models. This is expected since the expressions for the $D^{a}$ model are
considerably more involved, see Appendix B, Eqs. (B1), (B2) and (B3). On the
other hand, the computational time of all three reduced models is quite
similar.

Figure 4: Shallow shell: a) Geometry, material characteristics, and applied load. b) Comparison of results obtained with the present formulation and those from [56]. Table 1: Computational times for different constitutive models of the shallow shell example. For all models, the solution procedure requires 21 load increments and 52 arc-length iterations. Model | $D^{0}$ | $D^{1}$ | $D^{2}$ | $D^{a}$
---|---|---|---|---
Time [s] | 28.70 | 29.54 | 30.08 | 98.36

###  5.3 Pinched semi-cylinder

The pinched semi-cylinder as in Fig. 5a is the well-documented nonlinear
benchmark for thin shells which exhibit both large displacements and
rotations. A mesh with $20\times 20$ quadratic $C^{1}$ elements is utilized
for the discretization of the one half of the symmetric structure and the
deformed configuration is displayed in Fig. 5b. A comparison of the deflection
at point $A$ is shown in Fig. 6a. All results are in good agreement, but small
differences can be detected for LPF=1. Models $D^{a}$ and $D^{2}$ give
virtually indistinguishable results, similarly as models $D^{0}$ and $D^{1}$.
The largest relative difference between these models is close to 0.3 %.
Comparison with the reference solution given in [56] reveals that the largest
relative difference is less than 0.9 %.

In order to examine the influence of constitutive relation on this
characteristic displacement component, the shell thickness and the load
intensity are increased and three cases are considered: (i) $h/P=6/14000$,
(ii) $h/P=12/48000$, and (iii) $h/P=24/190000$. The differences between the
results obtained with the $D^{a}$ constitutive model and the reduced ones are
studied. The corresponding results are given in Fig. 6b, Fig. 6c, and Fig. 6d.
Evidently, all of the constitutive model types return virtually
indistinguishable results for the first part of the equilibrium path,
approximately $\textrm{LPF}<0.3$, regardless of the thickness. The
discrepancies become visible for the other part of the equilibrium path,
especially when $w_{A}>R$. This is due to the cumulative nonlinear effect and
also because the curviness value $Kh$ becomes significant at some parts of the
shell. It is emphasized that the curviness is measured at an arbitrary
configuration. Note that the $D^{2}$ model starts to visibly deviate from the
analytical $D^{a}$ model only for the shell with reference curviness
$Kh=24/101.6\approx 0.236$ which corresponds to initially strongly curved
shell. The distribution of curviness for two values of thickness and
$\textrm{LPF}=1$ is shown in Fig. 7 where the threshold value of $Kh=0.25$ is
utilized in order to depict this characteristic more clearly. To be precise,
all areas with $Kh\geq 0.25$ are colored with the same color. It is evident
how some parts of the shell become strongly curved.

Figure 5: Pinched semi-cylinder. a) Geometry, material characteristics, and
applied load. b) Deformed configuration for $\textrm{LPF}=1$. The colors show
the total value of the displacement vector. Figure 6: Pinched semi-cylinder.
a) Comparison of present formulation with [56] for thickness $h$ and load $P$;
Comparison of deflection for different constitutive models for: b) $2h$ and
$7P$; c) $4h$ and $24P$; d) $8h$ and $95P$. Insets show zoomed parts for
$\textrm{LPF}>0.8$. Figure 7: Pinched semi-cylinder. Distribution of curviness
for $\textrm{LPF}=1$ from three points of view for a shell with two different
values of thickness, i.e., $h$ and $2h$. For a clearer representation, a
threshold of $Kh=0.25$ is applied.

###  5.4 Pullout of a cylinder

The geometric nonlinear analysis of a cylinder with free ends, as in Fig. 8a,
is a well-established test for shell formulations. The shell is pulled out by
the two opposing forces and three displacement components are compared, [56].
Due to the symmetry of the problem, only one eighth of the shell is
discretized using a mesh with $60\times 60$ cubic $C^{2}$ elements. The
comparison with the reference solutions is displayed in Fig. 8b where the good
agreement is apparent. Since all constitutive models return virtually
indistinguishable results, their comparison is omitted. The deformed
configuration from three different points of view is displayed in Fig. 9.

The initial curviness of the observed shell is $Kh\approx 0.019$, but for
$\textrm{LPF}=1$ the maximum curviness of $Kh=1.42$ develops at point $A$. The
distribution of the curviness at the final configuration is given in Fig. 10a
using the threshold value of $Kh=0.25$. It is evident that some parts of this
shell become strongly curved locally. The change of the curviness with respect
to the load proportionality factor is shown in Fig. 10b for two characteristic
points.

In order to examine the influence of the curviness on the structural response,
strains at point $D$ are considered in detail. This point exhibits strong
local curviness, i.e., $Kh\approx 0.154$, and it is far enough away from the
applied force to avoid singularity issues. The four reference strains are
calculated using the four different constitutive models. Furthermore, Abaqus
results using a mesh with 15700 S4R elements are used for reference. All the
results are compared in Fig. 11 and interesting complex equilibrium paths are
detected.

Most importantly, the dominant reference strain is ${\kappa}_{1}$ and its
equilibrium path appears to be practically identical for all the considered
models, Fig. 11c. However, a close inspection of the zoomed part for
$\textrm{LPF}>0.8$ reveals discrepancies, Fig. 11c. The models $D^{a}$ and
$D^{2}$ return practically identical results, while the same observation holds
for the models $D^{0}$ and $D^{1}$. The relative difference of results for all
reduced models and Abaqus, with respect to $D^{a}$ model, is less than 1 % for
LPF=1. The other change of curvature, ${\kappa}_{2}$, is smaller for an order
of magnitude and the results obtained with different models start to visibly
deviate for $\textrm{LPF}>0.3$, Fig. 11d. Regarding the membrane reference
strains, all models return different equilibrium paths for ${\epsilon}_{11}$.
For ${\epsilon}_{22}$, the reduced models $D^{0}$ and $D^{1}$ differ from the
full analytical $D^{a}$ solution. This is expected since the deformed
curviness at point $D$ of this shell is dominated by the curvature component
along the $\xi$ direction $\left({b}^{1}_{1}\gg{b}^{2}_{2}\right)$, which
reflects on these membrane strain components, see Eq. (33). The differences in
comparison with Abaqus are evident and they are mainly attributed to the
application of the more rigorous constitutive models within the present
formulation. Additional influences are mesh density, continuity, and strain
averaging at the interelement boundaries, which are standard issues in FEA and
which are not of primary interest here.

Moreover, the $D^{a}$ constitutive model is employed for the calculation of
the strains and stresses at the outer surface at point $D$. The strains are
calculated using the approximative linear (55) and full nonlinear (51)
relations between reference and equidistant strains. The stresses are
calculated from these strains using (71). Comparison with strains and stresses
from Abaqus is given in Fig. 12. The dominant strain component is
${\epsilon}^{out}_{11}$ and the results obtained with Abaqus and Eq. (55)
differ from the rigorous ones, obtained with Eq. (51). The observed difference
with respect to Abaqus is mainly due to the difference in curvature change
${\kappa}_{1}$, Fig. 11c. The other component of strain is smaller for an
order of magnitude and all models return similar results, which is again a
consequence of the fact that ${b}^{1}_{1}\gg{b}^{2}_{2}$ at this point.
Regrading the stresses, discrepancies are observed for both components. The
noted differences with respect to Abaqus exist mainly due to the differences
in reference strains and assumed strain distributions over thickness. Another
influence comes from the calculation of the physical components, which are
here determined rigorously by the Eq. (56). Remarkably, these differences
practically vanish in the areas of shell that are not strongly curved. As an
example, let us observe the equidistant strains and stresses on the outer
surface at point $B$, Fig. 13. In this area, the curviness for
$\textrm{LPF}=1$ is small, $Kh\approx 0.018$, and the results for strains and
stresses agree well. This underscores that the observed discrepancies at the
point $C$ are mainly due to the large curviness of shell at that particular
point.

Figure 8: Pullout of a cylinder. a) Geometry, material characteristics, and
applied load. b) Comparison of displacements between two present models and
[56] for three displacement components.  Figure 9: Pullout of a cylinder.
Deformed configuration for $\textrm{LPF}=1$ from three different perspectives.
The color-coding represents the total value of the displacement vector. Figure
10: Pullout of a cylinder. a) The distribution of the curviness for
$\textrm{LPF}=1$ from three points of view. The threshold value of $Kh=0.25$
is utilized. b) Curviness at points $B$ and $D$ vs. LPF. Figure 11: Pullout of
a cylinder. Comparison of reference strains at point $D$ obtained with the
four different constitutive models and Abaqus: a) ${\epsilon}_{11}$; b)
${\epsilon}_{22}$; c) ${\kappa}_{1}$; d) ${\kappa}_{2}$. Figure 12: Pullout of
a cylinder. Comparison of strains and stresses at the outer surface at point
$D$ obtained with simplified linear and complete nonlinear relations between
the reference and equidistant strains, and Abaqus: a) ${\epsilon}^{out}_{11}$;
b) ${\epsilon}^{out}_{22}$; c) ${\sigma}^{11}_{out}$; d)
${\sigma}^{22}_{out}$. Figure 13: Pullout of a cylinder. Comparison of strains
and stresses at the outer surface at point $B$ obtained with simplified linear
and complete nonlinear relations between the reference and equidistant
strains, and Abaqus: a) ${\epsilon}^{out}_{11}$; b) ${\epsilon}^{out}_{22}$;
c) ${\sigma}^{11}_{out}$; d) ${\sigma}^{22}_{out}$.

###  5.5 Pinched cylinder

As the final example, the pinched cylinder with rigid diaphragms is considered
again, Fig. 14a. However, the emphasis is now on the geometric nonlinear
response, [56]. The reference curviness is $Kh=0.01$, as in the linear
example. A mesh with $50\times 50$ quadratic $C^{1}$ elements is employed for
the discretization of one eighth of the shell, similar to that in [19]. Two
characteristic displacement components are compared in Fig. 14b. For the
complete equilibrium path, there are no practical differences between the
results obtained with the different constitutive models, so they are not
displayed separately. Furthermore, the IGA results are in agreement with the
reference solution given in [56].

It is already noted that this example is very challenging for the solution
using the arc-length method. The tangent stiffness matrix has one negative
eigenvalue for $\textrm{LPF}>0.089$ which results with the erroneous negative
sign of the predictor solution. As noted, several types of arc-length
procedure are tested here, and the one implemented in the Abaqus proved to be
the most robust, [54].

Although the initial curviness of the observed shell is $Kh=0.01$, the finite
deformation resulted in the final configuration with a maximum curviness at
point $A$ of $Kh=0.346$. The distribution of the curviness for
$\textrm{LPF}=1$ is given in Fig. 16a. It is evident that some parts of this
shell become strongly curved locally. The change of curviness with respect to
the load proportionality factor is shown in Fig. 16b for two characteristic
points.

Next, the strains at point $C$ are considered in detail. There, the shell is
strongly curved, i.e., $Kh=0.101$ for $\textrm{LPF}=1$. The four reference
strains are calculated using the four different constitutive models and
compared in Fig. 17 with the results from Abaqus, which employs a mesh with
62800 S4R elements. In comparison to the previous example, even more complex
equilibrium paths are observed. The dominant reference strain is
${\kappa}_{2}$ and it is almost identical for all studied models, Fig. 17d. An
exception is the Abaqus model which gives different results for
$\textrm{LPF}>0.4$, with a maximum relative difference of nearly $3\%$ for
LPF=1. Considering the other component of the change of curvature, all models
return similar results, Fig. 17c. Regarding the reference membrane strains,
they are in agreement for $\textrm{LPF}<0.1$ but clear differences exist for
the greater part of equilibrium path which is especially pronounced for
${\epsilon}_{22}$, Fig. 17b. The constitutive models $D^{a}$ and $D^{2}$
return practically the same results, while there are small differences of
membrane strains for $\textrm{LPF}>0.7$. Additionally, the models $D^{0}$ and
$D^{1}$ are in close agreement.

Finally, the $D^{a}$ constitutive model is employed and the strains at the
outer surface are calculated using approximative (55) and exact (51) relations
between reference and equidistant strains, Fig. 18. Good agreement of all
models for $\textrm{LPF}<0.2$ is observed for the dominant strain,
${\epsilon}^{out}_{22}$. However, as the load and curviness increase, the
differences become noticeable. The results obtained with approximative
relation and Abaqus are close since the dominant reference strain is
${\kappa}_{2}$, Fig. 17c. The results of the other equidistant strain
component are indifferent with respect to the Eqs. (51) and (55) due to the
fact that ${b}^{1}_{1}\approx 0$ at point $C$.

Figure 14: Pinched cylinder. a) Geometry, material characteristics, and
applied load. b) Comparison of characteristic displacements components between
present formulation and [56].  Figure 15: Pinched cylinder. Deformed
configuration for $\textrm{LPF}=1$ from three perspectives. Color shows the
total value of displacement vector.  Figure 16: Pinched cylinder. a)
Distribution of curviness for $\textrm{LPF}=1$ from two points of view. The
threshold value of $Kh=0.25$ is utilized. b) Curviness at points $B$ and $C$
vs. LPF. Figure 17: Pinched cylinder. Comparison of reference strains at point
$C$ obtained with four constitutive models and Abaqus: a) ${\epsilon}_{11}$;
b) ${\epsilon}_{22}$; c) ${\kappa}_{1}$; d) ${\kappa}_{2}$. Figure 18: Pinched
cylinder. Comparison of strains at outer surface at point $C$ obtained with
linear and complete nonlinear relations between the reference and equidistant
strains, and Abaqus: a) ${\epsilon}^{out}_{11}$; b) ${\epsilon}^{out}_{22}$.

##  6 Conclusions

A rigorous geometric nonlinear IGA shell formulation is presented and
validated by thorough numerical experiments. A key finding is that the initial
configuration is not enough to determine the slenderness of the structure; all
subsequent deformed configurations must be considered. As a result, some
standard approximations of applied shell theories lead to erroneous results
for strongly curved shells.

We strictly derive the nonlinear equation of equilibrium by considering the
geometrically exact relations between the equidistant and reference strains of
a KL shell. The isogeometric approach is applied for the spatial
discretization and the arc-length method is employed as a solution procedure.
The geometric stiffness matrix is derived by the variation of strains with
respect to the metric. The symmetry of this matrix is maintained which proves
that the selected force and strain quantities are energetically conjugated.

The main contribution of the paper is the full analytical integration of the
virtual power along the thickness direction which is performed using the
reciprocal shift tensor. The resulting $D^{a}$ model is complex and thus,
computationally expensive. However, this constitutive model is exact, and
allows us to derive different reduced constitutive models $D^{0}$, $D^{1}$,
and $D^{2}$. Based on the present numerical experiments, it follows that the
approximative constitutive model $D^{2}$ is well-suited for the nonlinear
modeling of strongly curved shells. It is as efficient as the simple $D^{0}$
model and returns results which are virtually indistinguishable from the exact
ones. On the other hand, the displacement field of thin, small-curvature
shells, is practically indifferent to the more rigorous constitutive
relations, no matter how large the displacements and rotations are and the
simple decoupled equations of the $D^{0}$ model can return reasonably accurate
results. Regardless of the constitutive model used, nonlinear distribution of
strain through the thickness should be considered during the post-processing
phase in order to obtain accurate strain field of strongly curved shells.

A detailed inspection of the observed differences between IGA and Abaqus was
beyond the scope of the present research and it should be a subject of further
research. The development of an adaptive refinement scheme that applies higher
regularity on smooth parts of the deformed structure, and lower regularity at
the kinks and similar irregularities would be an interesting improvement.
Additionally, the different constitutive relations could be applied, depending
on the curviness of the observed area, like the blended shell approach
suggested in [58]. More accurate shell models should also include the change
of thickness as well, since this directly influences the curviness.

## Acknowledgments

The first author expresses his gratitude to Professor Sava Vukelić (1938-1992)
for the inspiring discussions and infectious enthusiasm regarding the
mechanics.

During this work, our beloved colleague and friend, Professor Gligor
Radenković (1956-2019), passed away. The second author acknowledges that his
unprecedented enthusiasm and love for mechanics were crucial for much of his
previous, present, and future research. His students and colleagues will
remember him by his incredibly infectious energy for both life and science.

The presented research is partially supported by the Ministry for Scientific-
Technological Development, Higher Education and Information Society of
Republic of Srpska through the project Dynamic stability of thin-walled
structures using the Isogeometric finite strip method. This support is
gratefully acknowledged.

We also acknowledge the support of the Austrian Science Fund (FWF): M 2806-N.

## Appendix A. Strain at an equidistant surface

If we insert Eqs. (48) and (49) into Eq. (47), the components of
$\bar{d}_{\alpha\beta}$ can be written as:

$\displaystyle 2\bar{d}_{\alpha\beta}$
$\displaystyle=\left({x}^{*n}_{,\alpha}-\zeta{b}^{*\nu}_{\alpha}{x}^{*n}_{,\nu}\right)\left({v}_{n,\beta}+\zeta{v}_{n,3\beta}\right)+\left({x}^{*n}_{,\beta}-\zeta{b}^{*\nu}_{\beta}{x}^{*n}_{,\nu}\right)\left({v}_{n,\alpha}+\zeta{v}_{n,3\alpha}\right)$
(A1)
$\displaystyle={x}^{*n}_{,\alpha}{v}_{n,\beta}+{x}^{*n}_{,\beta}{v}_{n,\alpha}-\zeta\left({b}^{*\nu}_{\alpha}{x}^{*n}_{,\nu}{v}_{n,\beta}+{b}^{*\nu}_{\beta}{x}^{*n}_{,\nu}{v}_{n,\alpha}-{x}^{*n}_{,\alpha}{v}_{n,3\beta}-{x}^{*n}_{,\beta}{v}_{n,3\alpha}\right)$
$\displaystyle\;\;\;\;-\zeta^{2}{x}^{*n}_{,\nu}\left({b}^{*\nu}_{\alpha}{v}_{n,3\beta}+{b}^{*\nu}_{\beta}{v}_{n,3\alpha}\right)$
$\displaystyle=2{d}_{\alpha\beta}-\zeta\left({b}^{*\nu}_{\alpha}{x}^{*n}_{,\nu}{v}_{n,\beta}+{b}^{*\nu}_{\beta}{x}^{*n}_{,\nu}{v}_{n,\alpha}\right.$
$\displaystyle\left.\;\;\;\;-{x}^{*n}_{,\alpha}{v}_{n,3\beta}-{x}^{*n}_{,\beta}{v}_{n,3\alpha}\right)-\zeta^{2}{x}^{*n}_{,\nu}\left({b}^{*\nu}_{\alpha}{v}_{n,3\beta}+{b}^{*\nu}_{\beta}{v}_{n,3\alpha}\right).$

By utilizing Eq. (49), the terms that are multiplied by $\zeta$ can be
transformed to:

$\displaystyle{x}^{*n}_{,\nu}\left({b}^{*\nu}_{\alpha}{v}_{n,\beta}+{b}^{*\nu}_{\beta}{v}_{n,\alpha}\right)-{x}^{*n}_{,\alpha}\left({\Gamma}^{*\mu}_{\beta\gamma}{x}^{*,\gamma}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{b}^{*\mu}_{\beta}{x}^{*,3}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{x}^{*,\mu}_{n}{v}_{k,\mu\beta}{x}^{*k}_{,3}\right.$
(A2)
$\displaystyle\left.+{x}^{*,\mu}_{n}{v}_{k,\mu}{b}^{*\nu}_{\beta}{x}^{*k}_{,\nu}\right)-{x}^{*n}_{,\beta}\left({\Gamma}^{*\mu}_{\alpha\gamma}{x}^{*,\gamma}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{b}^{*\mu}_{\alpha}{x}^{*,3}_{n}{v}_{k,\mu}{x}^{*k}_{,3}-{x}^{*,\mu}_{n}{v}_{k,\mu\alpha}{x}^{*k}_{,3}\right.$
$\displaystyle\left.+{x}^{*,\mu}_{n}{v}_{k,\mu}{b}^{*\nu}_{\alpha}{x}^{*k}_{,\nu}\right)={x}^{*n}_{,\nu}\left({b}^{*\nu}_{\alpha}{v}_{n,\beta}+{b}^{*\nu}_{\beta}{v}_{n,\alpha}\right)-{\Gamma}^{*\mu}_{\beta\alpha}{v}_{k,\mu}{x}^{*k}_{,3}+{v}_{k,\alpha\beta}{x}^{*k}_{,3}-{v}_{k,\alpha}{b}^{*\nu}_{\beta}{x}^{*k}_{,\nu}$
$\displaystyle-{\Gamma}^{*\mu}_{\beta\alpha}{v}_{k,\mu}{x}^{*k}_{,3}+{v}_{k,\alpha\beta}{x}^{*k}_{,3}-{v}_{k,\beta}{b}^{*\nu}_{\alpha}{x}^{*k}_{,\nu}=2{x}^{*k}_{,3}\left({v}_{k,\alpha\beta}-{\Gamma}^{*\mu}_{\alpha\beta}{v}_{k,\mu}\right).$

Analogously, the expression in Eq. (A1) that is multiplied by $\zeta^{2}$ can
be written as:

$\displaystyle{x}^{*n}_{,\nu}\left({b}^{*\nu}_{\alpha}{v}_{n,3\beta}+{b}^{*\nu}_{\beta}{v}_{n,3\alpha}\right)$
(A3)
$\displaystyle={b}^{*\nu}_{\alpha}\left({\Gamma}^{*\mu}_{\beta\nu}{v}_{k,\mu}{x}^{*k}_{,3}-{v}_{k,\nu\beta}{x}^{*k}_{,3}+{v}_{k,\nu}{b}^{*\mu}_{\beta}{x}^{*k}_{,\mu}\right)$
$\displaystyle\;\;\;\;+{b}^{*\nu}_{\beta}\left({\Gamma}^{*\mu}_{\alpha\nu}{v}_{k,\mu}{x}^{*k}_{,3}-{v}_{k,\nu\alpha}{x}^{*k}_{,3}+{v}_{k,\nu}{b}^{*\mu}_{\alpha}{x}^{*k}_{,\mu}\right)$
$\displaystyle=\left({b}^{*\nu}_{\alpha}{b}^{*\mu}_{\beta}+{b}^{*\nu}_{\beta}{b}^{*\mu}_{\alpha}\right){v}_{k,\nu}{x}^{*k}_{,\mu}-{b}^{*\nu}_{\alpha}{x}^{*k}_{,3}\left({v}_{k,\nu\beta}-{\Gamma}^{*\mu}_{\nu\beta}{v}_{k,\mu}\right)-{b}^{*\nu}_{\beta}{x}^{*k}_{,3}\left({v}_{k,\nu\alpha}-{\Gamma}^{*\mu}_{\nu\alpha}{v}_{k,\mu}\right).$

If we introduce designations for the following quantities:

$\dot{\kappa}_{\alpha\beta}={x}^{*k}_{,3}\left({v}_{k,\alpha\beta}-{\Gamma}^{*\mu}_{\alpha\beta}{v}_{k,\mu}\right),$
(A4)
$2\dot{\chi}_{\alpha\beta}=\left({b}^{*\nu}_{\alpha}{b}^{*\mu}_{\beta}+{b}^{*\nu}_{\beta}{b}^{*\mu}_{\alpha}\right){v}_{k,\nu}{x}^{*k}_{,\mu}-{b}^{*\nu}_{\alpha}\dot{\kappa}_{\nu\beta}-{b}^{*\nu}_{\beta}\dot{\kappa}_{\nu\alpha},$
(A5)

the strain rate at an arbitrary point can be written as:

$\bar{d}_{\alpha\beta}={d}_{\alpha\beta}-\zeta\dot{\kappa}_{\alpha\beta}-{\zeta}^{2}\dot{\chi}_{\alpha\beta}.$
(A6)

It is evident that $\dot{\chi}_{\alpha\beta}$ is a complicated quantity which
consists of the products of curvatures, curvature change rates, and gradients
of velocity. If we notice that $\nu$ and $\mu$ are dummy indices, Eq. (A5)
reduces to:

$2\dot{\chi}_{\alpha\beta}={b}^{*\nu}_{\alpha}{b}^{*\mu}_{\beta}\left({v}_{k,\nu}{x}^{*k}_{,\mu}+{v}_{k,\mu}{x}^{*k}_{,\nu}\right)-{b}^{*\nu}_{\alpha}\dot{\kappa}_{\nu\beta}-{b}^{*\nu}_{\beta}\dot{\kappa}_{\nu\alpha}=2{b}^{*\nu}_{\alpha}{b}^{*\mu}_{\beta}{d}_{\nu\mu}-{b}^{*\nu}_{\alpha}\dot{\kappa}_{\nu\beta}-{b}^{*\nu}_{\beta}\dot{\kappa}_{\nu\alpha}.$
(A7)

Using this relation, the strain rate at an arbitrary point can be expressed as
a function of the strain and curvature change rates of the shell midsurface
as:

$\bar{d}_{\alpha\beta}=\langle\bar{C}^{*\mu}_{\alpha}\left[\left({\delta}^{\nu}_{\beta}+\zeta{b}^{*\nu}_{\beta}\right){d}_{\mu\nu}-\zeta\dot{\kappa}_{\mu\beta}\right]\rangle^{sym},\quad\bar{C}^{*\mu}_{\alpha}={\delta}^{\mu}_{\alpha}-\zeta{b}^{*\mu}_{\alpha},$
(A8)

where the symmetry of the expression refers to the symmetry by indices
$\alpha$ and $\beta$, and $\mu$ and $\nu$, that is:

$2\bar{d}_{\alpha\beta}=\bar{C}^{*\mu}_{\alpha}\left[\left({\delta}^{\nu}_{\beta}+\zeta{b}^{*\nu}_{\beta}\right){d}_{\mu\nu}-\zeta\dot{\kappa}_{\mu\beta}\right]+\bar{C}^{*\nu}_{\beta}\left[\left({\delta}^{\mu}_{\alpha}+\zeta{b}^{*\mu}_{\alpha}\right){d}_{\mu\nu}-\zeta\dot{\kappa}_{\nu\alpha}\right].$
(A9)

By introducing:

$\displaystyle{A}^{\mu\nu}_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\left[\bar{C}^{*\mu}_{\alpha}\left({\delta}^{\nu}_{\beta}+\zeta{b}^{*\nu}_{\beta}\right)+\bar{C}^{*\nu}_{\beta}\left({\delta}^{\mu}_{\alpha}+\zeta{b}^{*\mu}_{\alpha}\right)\right]={\delta}^{\mu}_{\alpha}{\delta}^{\nu}_{\beta}-\zeta^{2}{b}^{*\mu}_{\alpha}{b}^{*\nu}_{\beta},$
(A10) $\displaystyle{B}^{\mu\nu}_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\left({\delta}^{\nu}_{\alpha}\bar{C}^{*\mu}_{\beta}+{\delta}^{\mu}_{\beta}\bar{C}^{*\nu}_{\alpha}\right),$

we finally obtain the compact form for the strain rate at an equidistant
surface:

$\bar{d}_{\alpha\beta}={A}^{\mu\nu}_{\alpha\beta}{d}_{\mu\nu}-\zeta{B}^{\mu\nu}_{\alpha\beta}\dot{\kappa}_{\mu\nu}.$
(A11)

## Appendix B. Closed-form of the constitutive matrix

The elements of the constitutive matrix given by Eq. (74) are now presented in
an expanded form.

$\displaystyle{D}^{\gamma\lambda\chi\omega}_{M}$
$\displaystyle=\left({I}_{0}-2T{I}_{1}+{T}^{2}{I}_{2}\right){D}^{\gamma\lambda\chi\omega}$
(B1)
$\displaystyle\;\;\;\;+\frac{1}{2}\left({b}^{\lambda}_{\nu}{K}^{1}_{M}+{b}^{\lambda}_{\mu}{b}^{\mu}_{\nu}{K}^{2}_{M}\right){D}^{\gamma\nu\chi\omega}+\frac{1}{2}\left({b}^{\omega}_{\varphi}{K}^{1}_{M}+{b}^{\omega}_{\mu}{b}^{\mu}_{\varphi}{K}^{2}_{M}\right){D}^{\gamma\lambda\chi\varphi}$
$\displaystyle\;\;\;\;+\frac{1}{2}\left({b}^{\chi}_{\varphi}{K}^{1}_{M}+{b}^{\chi}_{\mu}{b}^{\mu}_{\varphi}{K}^{2}_{M}\right){D}^{\gamma\lambda\varphi\omega}+\frac{1}{2}\left({b}^{\gamma}_{\nu}{K}^{1}_{M}+{b}^{\gamma}_{\mu}{b}^{\mu}_{\nu}{K}^{2}_{M}\right){D}^{\nu\lambda\chi\omega}$
$\displaystyle\;\;\;\;+\frac{1}{4}\left[{b}^{\omega}_{\varphi}{b}^{\lambda}_{\nu}{K}^{3}_{M}+\left({b}^{\omega}_{\mu}{b}^{\mu}_{\varphi}{b}^{\lambda}_{\nu}+{b}^{\omega}_{\varphi}{b}^{\lambda}_{\mu}{b}^{\mu}_{\nu}\right){K}^{4}_{M}+{b}^{\omega}_{\epsilon}{b}^{\epsilon}_{\varphi}{b}^{\lambda}_{\mu}{b}^{\mu}_{\nu}{I}_{4}\right]{D}^{\gamma\nu\chi\varphi}$
$\displaystyle\;\;\;\;+\frac{1}{4}\left[{b}^{\omega}_{\varphi}{b}^{\gamma}_{\nu}{K}^{3}_{M}+\left({b}^{\omega}_{\mu}{b}^{\mu}_{\varphi}{b}^{\gamma}_{\nu}+{b}^{\omega}_{\varphi}{b}^{\gamma}_{\mu}{b}^{\mu}_{\nu}\right){K}^{4}_{M}+{b}^{\omega}_{\epsilon}{b}^{\epsilon}_{\varphi}{b}^{\gamma}_{\mu}{b}^{\mu}_{\nu}{I}_{4}\right]{D}^{\nu\lambda\chi\varphi}$
$\displaystyle\;\;\;\;+\frac{1}{4}\left[{b}^{\chi}_{\varphi}{b}^{\lambda}_{\nu}{K}^{3}_{M}+\left({b}^{\chi}_{\mu}{b}^{\mu}_{\varphi}{b}^{\lambda}_{\nu}+{b}^{\chi}_{\varphi}{b}^{\lambda}_{\mu}{b}^{\mu}_{\nu}\right){K}^{4}_{M}+{b}^{\chi}_{\epsilon}{b}^{\epsilon}_{\varphi}{b}^{\lambda}_{\mu}{b}^{\mu}_{\nu}{I}_{4}\right]{D}^{\gamma\nu\varphi\omega}$
$\displaystyle\;\;\;\;+\frac{1}{4}\left[{b}^{\chi}_{\varphi}{b}^{\gamma}_{\nu}{K}^{3}_{M}+\left({b}^{\chi}_{\mu}{b}^{\mu}_{\varphi}{b}^{\gamma}_{\nu}+{b}^{\chi}_{\varphi}{b}^{\gamma}_{\mu}{b}^{\mu}_{\nu}\right){K}^{4}_{M}+{b}^{\chi}_{\epsilon}{b}^{\epsilon}_{\varphi}{b}^{\gamma}_{\mu}{b}^{\mu}_{\nu}{I}_{4}\right]{D}^{\nu\lambda\varphi\omega},$
$\displaystyle-{D}^{\lambda\gamma\chi\omega}_{MB}$
$\displaystyle=\left({I}_{1}-2T{I}_{2}+{T}^{2}{I}_{3}\right){D}^{\lambda\gamma\chi\omega}+\frac{1}{2}{K}^{2}_{M}\left({b}^{\gamma}_{\nu}{D}^{\lambda\nu\chi\omega}+{b}^{\lambda}_{\nu}{D}^{\nu\gamma\chi\omega}\right)$
(B2)
$\displaystyle\;\;\;\;+\frac{1}{2}\left({b}^{\omega}_{\varphi}{K}^{1}_{MB}+{b}^{\omega}_{\mu}{b}^{\mu}_{\varphi}{K}^{2}_{MB}\right){D}^{\lambda\gamma\chi\varphi}+\frac{1}{2}\left({b}^{\chi}_{\varphi}{K}^{1}_{MB}+{b}^{\chi}_{\mu}{b}^{\mu}_{\varphi}{K}^{2}_{MB}\right){D}^{\lambda\gamma\varphi\omega}$
$\displaystyle\;\;\;\;+\frac{1}{4}{b}^{\gamma}_{\nu}\left({b}^{\omega}_{\varphi}{K}^{3}_{MB}+{b}^{\omega}_{\mu}{b}^{\mu}_{\varphi}{I}_{4}\right){D}^{\lambda\nu\chi\varphi}+\frac{1}{4}{b}^{\lambda}_{\nu}\left({b}^{\omega}_{\varphi}{K}^{3}_{MB}+{b}^{\omega}_{\mu}{b}^{\mu}_{\varphi}{I}_{4}\right){D}^{\nu\gamma\chi\varphi}$
$\displaystyle\;\;\;\;+\frac{1}{4}{b}^{\gamma}_{\nu}\left({b}^{\chi}_{\varphi}{K}^{3}_{MB}+{b}^{\chi}_{\mu}{b}^{\mu}_{\varphi}{I}_{4}\right){D}^{\lambda\nu\varphi\omega}+\frac{1}{4}{b}^{\lambda}_{\nu}\left({b}^{\chi}_{\varphi}{K}^{3}_{MB}+{b}^{\chi}_{\mu}{b}^{\mu}_{\varphi}{I}_{4}\right){D}^{\nu\gamma\varphi\omega},$
$\displaystyle{D}^{\lambda\gamma\chi\omega}_{B}$
$\displaystyle=\left({I}_{2}-2T{I}_{3}+{T}^{2}{I}_{4}\right){D}^{\lambda\gamma\chi\omega}$
(B3)
$\displaystyle\;\;\;\;+\frac{1}{2}\left({I}_{3}-T{I}_{4}\right)\left({b}^{\lambda}_{\nu}{D}^{\nu\gamma\chi\omega}+{b}^{\gamma}_{\nu}{D}^{\lambda\nu\chi\omega}+{b}^{\chi}_{\nu}{D}^{\lambda\gamma\nu\omega}+{b}^{\omega}_{\nu}{D}^{\lambda\gamma\chi\nu}\right)$
$\displaystyle\;\;\;\;+\frac{1}{4}{I}_{4}\left({b}^{\lambda}_{\nu}{b}^{\chi}_{\varphi}{D}^{\nu\gamma\varphi\omega}+{b}^{\lambda}_{\nu}{b}^{\omega}_{\varphi}{D}^{\nu\gamma\chi\varphi}+{b}^{\gamma}_{\nu}{b}^{\chi}_{\varphi}{D}^{\lambda\nu\varphi\omega}+{b}^{\gamma}_{\nu}{b}^{\omega}_{\varphi}{D}^{\lambda\nu\chi\varphi}\right).$

Here, the following designations are introduced:

$\displaystyle{K}^{1}_{M}$
$\displaystyle=2{I}_{1}-3T{I}_{2}+{T}^{2}{I}_{3},\quad{K}^{2}_{M}=\
{I}_{2}-T{I}_{3},\;{K}^{3}_{M}=4{I}_{2}-4T{I}_{3}+{T}^{2}{I}_{4},$ (B4)
$\displaystyle{K}^{4}_{M}$
$\displaystyle=2{I}_{3}-T{I}_{4},\quad{K}^{1}_{MB}=2{I}_{2}-3T{I}_{3}+{T}^{2}{I}_{4},\quad{K}^{2}_{MB}={I}_{3}-T{I}_{4},$
$\displaystyle{K}^{3}_{MB}$ $\displaystyle=2{I}_{3}-T{I}_{4},$

while $T$ denotes the trace of the curvature tensor. Finally, utilizing
symbolic capabilities of Wolfram Mathematica, following closed-form solutions
for the integrals $I_{0}-I_{4}$ are obtained:

$\displaystyle p$ $\displaystyle=\sqrt{4b-{T}^{2}},\quad
m=\frac{bh-T}{p},\quad n=\frac{bh+T}{p},$ (B5) $\displaystyle r$
$\displaystyle=4+b{h}^{2}+2hT,\quad s=4+bh^{2}-2hT,$ $\displaystyle{I}_{0}$
$\displaystyle=\int_{-h/2}^{h/2}\frac{1}{1-\zeta
T+\zeta^{2}b}\differential{\zeta}=\frac{2}{p}\left[\arctan(m)+\arctan(n)\right],$
$\displaystyle{I}_{1}$ $\displaystyle=\int_{-h/2}^{h/2}\frac{\zeta}{1-\zeta
T+\zeta^{2}b}\differential{\zeta}$
$\displaystyle=\frac{1}{2bp}\left\\{2T\left[\arctan(m)+\arctan(n)\right]+p\ln(\frac{8+2bh^{2}}{r}-1)\right\\},$
$\displaystyle{I}_{2}$
$\displaystyle=\int_{-h/2}^{h/2}\frac{\zeta^{2}}{1-\zeta
T+\zeta^{2}b}\differential{\zeta}$
$\displaystyle=\frac{1}{b^{2}p}\left\\{bhp-\left(2b-T^{2}\right)\left[\arctan(m)+\arctan(n)\right]-Tp\arctan(\frac{2hT}{4+bh^{2}})\right\\},$
$\displaystyle{I}_{3}$
$\displaystyle=\int_{-h/2}^{h/2}\frac{\zeta^{3}}{1-\zeta
T+\zeta^{2}b}\differential{\zeta}$
$\displaystyle=\frac{1}{2b^{3}p}\left\\{\left(6bT-2T^{3}\right)\arctan(-n)+\left(2T^{3}-6bT\right)\arctan(m)\right.$
$\displaystyle\left.\;\;\;\;+p\left[2bhT+\left(T^{2}-b\right)\ln(s)+\left(b-T^{2}\right)\ln(r)\right]\right\\},$
$\displaystyle{I}_{4}$
$\displaystyle=\int_{-h/2}^{h/2}\frac{\zeta^{4}}{1-\zeta
T+\zeta^{2}b}\differential{\zeta}$
$\displaystyle=\frac{1}{12b^{4}p}\left\\{-12\left(2b^{2}-4bT^{2}+T^{4}\right)\arctan(-n)+12\left(2b^{2}-4bT^{2}+T^{4}\right)\arctan(m)\right.$
$\displaystyle\left.\;\;\;\;+p\left[bh\left(-12b+b^{2}h^{2}+12T^{2}\right)+6T\left(T^{2}-2b\right)\ln(s)+6\left(2bT-h^{3}\right)\ln(r)\right]\right\\}.$

However, it should be noted that these closed-form integrals can become
unstable for small-curvature shells and an appropriate switch in the algorithm
is required in order to include reduced expressions for which the standard
thin shell assumptions are valid. Therefore, the presented expressions are
ideally suited for strongly curved parts of a shell, but the exact limit
between small and large curvature shells is arbitrary and some criteria must
be introduced, e.g., $Kh>0.05$.

## Appendix C. Different reduced constitutive models

Three simple reduced models are utilized in this paper, $D^{0}$, $D^{1}$, and
$D^{2}$. $D^{0}$ model is analogous to the flat plate model for which:

$\begin{gathered}\frac{1}{1-\zeta T+\zeta^{2}b}\approx
1\Rightarrow{I}_{0}=h,\quad{I}_{1}=0,{I}_{2}=\frac{h^{3}}{12},\\\
{D}^{\lambda\gamma\chi\omega}_{M}=h{D}^{\lambda\gamma\chi\omega},\quad{D}^{\lambda\gamma\chi\omega}_{MB}=0,\quad{D}^{\lambda\gamma\chi\omega}_{B}=\frac{h^{3}}{12}{D}^{\lambda\gamma\chi\omega}.\end{gathered}$
(C1)

The other two models are based on the first-order Taylor approximation:

$\frac{1}{1-\zeta T+\zeta^{2}b}\approx 1+\zeta
T\Rightarrow{I}_{0}=h,\quad{I}_{1}=\frac{h^{3}}{12}T=T{I}_{2},\quad{I}_{2}=\frac{h^{3}}{12},$
(C2)

where the ${D}_{M}$ and ${D}_{B}$ terms are the same as those of the ${D}^{0}$
model, but the coupling term ${D}_{MB}$ differs:

$\displaystyle\textrm{for}\;D^{1}\;\textrm{model:}\quad-{D}^{\lambda\gamma\chi\omega}_{MB}$
$\displaystyle=-\frac{h^{3}}{12}T{D}^{\lambda\gamma\chi\omega},$ (C3)
$\displaystyle\textrm{for}\;D^{2}\;\textrm{model:}\quad-{D}^{\lambda\gamma\chi\omega}_{MB}$
$\displaystyle=-\frac{h^{3}}{12}T{D}^{\lambda\gamma\chi\omega}$
$\displaystyle\;\;\;\;+\frac{h^{3}}{24}\left({b}^{\gamma}_{\nu}{D}^{\lambda\nu\chi\omega}+{b}^{\lambda}_{\nu}{D}^{\nu\gamma\chi\omega}+2{b}^{\omega}_{\varphi}{D}^{\lambda\gamma\chi\varphi}+2{b}^{\chi}_{\varphi}{D}^{\lambda\gamma\varphi\omega}\right).$

## References

  * [1] P. Naghdi, “The Theory of Shells and Plates,” in Truesdell C. (Eds) Linear Theories of Elasticity and Thermoelasticity, pp. 425–640, Berlin, Heidelberg: Springer, 1973.
  * [2] M. Bischoff, K.-U. Bletzinger, W. A. Wall, and E. Ramm, “Models and Finite Elements for Thin-Walled Structures,” in Encyclopedia of Computational Mechanics, American Cancer Society, 2004.
  * [3] H.-M. Jeon, Y. Lee, P.-S. Lee, and K.-J. Bathe, “The MITC3+ shell element in geometric nonlinear analysis,” Computers & Structures, vol. 146, pp. 91–104, Jan. 2015.
  * [4] Y. Ko, P.-S. Lee, and K.-J. Bathe, “The MITC4+ shell element in geometric nonlinear analysis,” Computers & Structures, vol. 185, pp. 1–14, June 2017.
  * [5] S. Sadamoto, M. Ozdemir, S. Tanaka, T. Q. Bui, and S. Okazawa, “Finite rotation meshfree formulation for geometrically nonlinear analysis of flat, curved and folded shells,” International Journal of Non-Linear Mechanics, vol. 119, p. 103300, Mar. 2020.
  * [6] M.-J. Dai, S. Tanaka, S. Sadamoto, T. Yu, and T. Q. Bui, “Advanced reproducing kernel meshfree modeling of cracked curved shells for mixed-mode stress resultant intensity factors,” Engineering Fracture Mechanics, vol. 233, p. 107012, June 2020.
  * [7] H. Kurtaran, “Geometrically nonlinear transient analysis of moderately thick laminated composite shallow shells with generalized differential quadrature method,” Composite Structures, vol. 125, pp. 605–614, July 2015.
  * [8] A. Borković, S. Kovačević, D. D. Milašinović, G. Radenković, O. Mijatović, and V. Golubović-Bugarski, “Geometric nonlinear analysis of prismatic shells using the semi-analytical finite strip method,” Thin-Walled Structures, vol. 117, pp. 63–88, Aug. 2017.
  * [9] B. Wu, A. Pagani, W. Q. Chen, and E. Carrera, “Geometrically nonlinear refined shell theories by Carrera Unified Formulation,” Mechanics of Advanced Materials and Structures, vol. 0, pp. 1–21, Dec. 2019.
  * [10] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement,” Computer Methods in Applied Mechanics and Engineering, vol. 194, pp. 4135–4195, Oct. 2005.
  * [11] J. Kiendl, K. U. Bletzinger, J. Linhard, and R. Wüchner, “Isogeometric shell analysis with Kirchhoff–Love elements,” Computer Methods in Applied Mechanics and Engineering, vol. 198, pp. 3902–3914, Nov. 2009.
  * [12] D. J. Benson, Y. Bazilevs, M. C. Hsu, and T. J. R. Hughes, “A large deformation, rotation-free, isogeometric shell,” Computer Methods in Applied Mechanics and Engineering, vol. 200, pp. 1367–1378, Mar. 2011.
  * [13] J. Kiendl, M.-C. Hsu, M. C. H. Wu, and A. Reali, “Isogeometric Kirchhoff–Love shell formulations for general hyperelastic materials,” Computer Methods in Applied Mechanics and Engineering, vol. 291, pp. 280–303, July 2015.
  * [14] N. Nguyen-Thanh, N. Valizadeh, M. N. Nguyen, H. Nguyen-Xuan, X. Zhuang, P. Areias, G. Zi, Y. Bazilevs, L. De Lorenzis, and T. Rabczuk, “An extended isogeometric thin shell analysis based on Kirchhoff–Love theory,” Computer Methods in Applied Mechanics and Engineering, vol. 284, pp. 265–291, Feb. 2015.
  * [15] A. Riffnaller-Schiefer, U. H. Augsdörfer, and D. W. Fellner, “Isogeometric shell analysis with NURBS compatible subdivision surfaces,” Applied Mathematics and Computation, vol. 272, pp. 139–147, Jan. 2016.
  * [16] K. Luo, C. Liu, Q. Tian, and H. Hu, “An efficient model reduction method for buckling analyses of thin shells based on IGA,” Computer Methods in Applied Mechanics and Engineering, vol. 309, pp. 243–268, Sept. 2016.
  * [17] M. Ambati, J. Kiendl, and L. De Lorenzis, “Isogeometric Kirchhoff–Love shell formulation for elasto-plasticity,” Computer Methods in Applied Mechanics and Engineering, vol. 340, pp. 320–339, Oct. 2018.
  * [18] L. Leonetti, D. Magisano, A. Madeo, G. Garcea, J. Kiendl, and A. Reali, “A simplified Kirchhoff–Love large deformation model for elastic shells and its effective isogeometric formulation,” Computer Methods in Applied Mechanics and Engineering, vol. 354, pp. 369–396, Sept. 2019\.
  * [19] T. X. Duong, F. Roohbakhshan, and R. A. Sauer, “A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries,” Computer Methods in Applied Mechanics and Engineering, vol. 316, pp. 43–83, Apr. 2017.
  * [20] T. N. Nguyen, S. Lee, P.-C. Nguyen, H. Nguyen-Xuan, and J. Lee, “Geometrically nonlinear postbuckling behavior of imperfect FG-CNTRC shells under axial compression using isogeometric analysis,” European Journal of Mechanics - A/Solids, vol. 84, p. 104066, Nov. 2020.
  * [21] J. Kiendl, Y. Bazilevs, M. C. Hsu, R. Wüchner, and 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, vol. 199, pp. 2403–2416, Aug. 2010.
  * [22] L. Coox, F. Maurin, F. Greco, E. Deckers, D. Vandepitte, and W. Desmet, “A flexible approach for coupling NURBS patches in rotationless isogeometric analysis of Kirchhoff–Love shells,” Computer Methods in Applied Mechanics and Engineering, vol. 325, pp. 505–531, Oct. 2017.
  * [23] F. Cirak and Q. Long, “Subdivision shells with exact boundary control and non-manifold geometry,” International Journal for Numerical Methods in Engineering, vol. 88, no. 9, pp. 897–923, 2011.
  * [24] H. Casquero, L. Liu, Y. Zhang, A. Reali, J. Kiendl, and H. Gomez, “Arbitrary-degree T-splines for isogeometric analysis of fully nonlinear Kirchhoff–Love shells,” Computer-Aided Design, vol. 82, pp. 140–153, Jan. 2017.
  * [25] B. Marussig and T. J. R. Hughes, “A Review of Trimming in Isogeometric Analysis: Challenges, Data Exchange and Simulation Aspects,” Archives of Computational Methods in Engineering, vol. 25, pp. 1059–1127, Nov. 2018.
  * [26] R. Schmidt, R. Wüchner, and K.-U. Bletzinger, “Isogeometric analysis of trimmed NURBS geometries,” Computer Methods in Applied Mechanics and Engineering, vol. 241-244, pp. 93–111, Oct. 2012.
  * [27] M. Breitenberger, A. Apostolatos, B. Philipp, R. Wüchner, and K. U. Bletzinger, “Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures,” Computer Methods in Applied Mechanics and Engineering, vol. 284, pp. 401–457, Feb. 2015.
  * [28] Y. Guo and M. Ruess, “Weak Dirichlet boundary conditions for trimmed thin isogeometric shells,” Computers & Mathematics with Applications, vol. 70, pp. 1425–1440, Oct. 2015.
  * [29] Y. Guo, M. Ruess, and D. Schillinger, “A parameter-free variational coupling approach for trimmed isogeometric thin shells,” Computational Mechanics, vol. 59, pp. 693–715, Apr. 2017.
  * [30] Y. Guo, J. Heller, T. J. R. Hughes, M. Ruess, and D. Schillinger, “Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format,” Computer Methods in Applied Mechanics and Engineering, vol. 336, pp. 39–79, July 2018\.
  * [31] D. J. Benson, Y. Bazilevs, M. C. Hsu, and T. J. R. Hughes, “Isogeometric shell analysis: The Reissner–Mindlin shell,” Computer Methods in Applied Mechanics and Engineering, vol. 199, pp. 276–289, Jan. 2010\.
  * [32] P. M. Sobota, W. Dornisch, R. Müller, and S. Klinkel, “Implicit dynamic analysis using an isogeometric Reissner–Mindlin shell formulation,” International Journal for Numerical Methods in Engineering, vol. 110, no. 9, pp. 803–825, 2017.
  * [33] P.-S. Lee and K.-J. Bathe, “Development of MITC isotropic triangular shell finite elements,” Computers & Structures, vol. 82, pp. 945–962, May 2004\.
  * [34] W. Dornisch, S. Klinkel, and B. Simeon, “Isogeometric Reissner–Mindlin shell analysis with exactly calculated director vectors,” Computer Methods in Applied Mechanics and Engineering, vol. 253, pp. 491–504, Jan. 2013.
  * [35] J. Kiendl, E. Marino, and L. De Lorenzis, “Isogeometric collocation for the Reissner–Mindlin shell problem,” Computer Methods in Applied Mechanics and Engineering, vol. 325, pp. 645–665, Oct. 2017.
  * [36] B. Oesterle, R. Sachse, E. Ramm, and M. Bischoff, “Hierarchic isogeometric large rotation shell elements including linearized transverse shear parametrization,” Computer Methods in Applied Mechanics and Engineering, vol. 321, pp. 383–405, July 2017.
  * [37] W. Li, N. Nguyen-Thanh, and K. Zhou, “Geometrically nonlinear analysis of thin-shell structures based on an isogeometric-meshfree coupling approach,” Computer Methods in Applied Mechanics and Engineering, vol. 336, pp. 111–134, July 2018.
  * [38] A. Cazzani, M. Malagù, E. Turco, and F. Stochino, “Constitutive models for strongly curved beams in the frame of isogeometric analysis:,” Mathematics and Mechanics of Solids, Mar. 2015.
  * [39] A. Borković, S. Kovačević, G. Radenković, S. Milovanović, and M. Guzijan-Dilber, “Rotation-free isogeometric analysis of an arbitrarily curved plane Bernoulli–Euler beam,” Computer Methods in Applied Mechanics and Engineering, vol. 334, pp. 238–267, June 2018.
  * [40] A. Borković, S. Kovačević, G. Radenković, S. Milovanović, and D. Majstorović, “Rotation-free isogeometric dynamic analysis of an arbitrarily curved plane Bernoulli-Euler beam,” Engineering Structures, vol. 181, pp. 192–215, Feb. 2019.
  * [41] G. Radenković and A. Borković, “Linear static isogeometric analysis of an arbitrarily curved spatial Bernoulli–Euler beam,” Computer Methods in Applied Mechanics and Engineering, vol. 341, pp. 360–396, Nov. 2018.
  * [42] G. Radenković and A. Borković, “On the analytical approach to the linear analysis of an arbitrarily curved spatial Bernoulli–Euler beam,” Applied Mathematical Modelling, vol. 77, pp. 1603–1624, Jan. 2020.
  * [43] A. Hamdouni and O. Millet, “Classification of thin shell models deduced from the nonlinear three-dimensional elasticity. Part II : The strongly curved shells,” Archives of Mechanics, vol. 55, pp. 177–219, Jan. 2003.
  * [44] D. Gamby and B. Na, “Generalized Kirchhoff-Love assumptions for thick shells: Application to thermal stresses in rectilinear tubes,” Acta Mechanica, vol. 74, pp. 95–106, Oct. 1988.
  * [45] G. Radenković, Finite rotation and finite strain isogeometric structural analysis (in Serbian).  Faculty of Architecture Belgrad, 2017.
  * [46] S. Naomis and P. C. M. Lau, Computational Tensor Analysis of Shell Structures.  Lecture Notes in Engineering, Berlin Heidelberg: Springer-Verlag, 1990.
  * [47] E. H. C. Baker, “Shell analysis manual,” tech. rep., Apr. 1968.
  * [48] P. Wriggers, Nonlinear Finite Element Methods.  Berlin Heidelberg: Springer-Verlag, 2008.
  * [49] G. C. Johnson and D. J. Bammann, “A discussion of stress rates in finite deformation problems,” International Journal of Solids and Structures, vol. 20, pp. 725–737, Jan. 1984.
  * [50] K. J. Bathe, Finite Element Procedures.  Boston, Mass.: Klaus-Jurgen Bathe, Feb. 2007.
  * [51] W. T. Koiter and J. G. Simmonds, “Foundations of shell theory,” in Theoretical and Applied Mechanics (E. Becker and G. K. Mikhailov, eds.), IUTAM Symposia, (Berlin, Heidelberg), pp. 150–176, Springer, 1973.
  * [52] M. A. Crisfield, “A fast incremental/iterative solution procedure that handles “snap-through”,” Computers & Structures, vol. 13, pp. 55–62, June 1981\.
  * [53] M. Ritto-Corrêa and D. Camotim, “On the arc-length and other quadratic control methods: Established, less known and new implementation procedures,” Computers & Structures, vol. 86, pp. 1353–1368, June 2008\.
  * [54] M. Smith, “ABAQUS/Standard User’s Manual, Version 6.9,” 2009.
  * [55] A. Stahl, T. Kvamsdal, and C. Schellewald, “Post-processing and visualization techniques for isogeometric analysis results,” Computer Methods in Applied Mechanics and Engineering, vol. 316, pp. 880–943, Apr. 2017.
  * [56] K. Y. Sze, X. H. Liu, and S. H. Lo, “Popular benchmark problems for geometric nonlinear analysis of shells,” Finite Elements in Analysis and Design, vol. 40, pp. 1551–1569, July 2004.
  * [57] F. Auricchio, F. Calabrò, T. J. R. Hughes, A. Reali, and G. Sangalli, “A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis,” Computer Methods in Applied Mechanics and Engineering, vol. 249-252, pp. 15–27, Dec. 2012.
  * [58] D. J. Benson, S. Hartmann, Y. Bazilevs, M. C. Hsu, and T. J. R. Hughes, “Blended isogeometric shells,” Computer Methods in Applied Mechanics and Engineering, vol. 255, pp. 133–146, Mar. 2013.