# An embedding-aware continuum thin shell formulation

Abhishek Ghosh  James Watt School of Engineering, University of Glasgow,
Glasgow G12 8LT, United Kingdom Andrew McBride  James Watt School of
Engineering, University of Glasgow, Glasgow G12 8LT, United Kingdom Zhaowei
Liu  College of Mechanics and Engineering Science, Hohai University, Nanjing,
211100, China  Luca Heltai  University of Pisa, Via Buonarroti 1/c 56127 Pisa
(PI), Italy  Paul Steinmann  James Watt School of Engineering, University of
Glasgow, Glasgow G12 8LT, United Kingdom Institute of Applied Mechanics,
Friedrich-Alexander Universität Erlangen-Nürnberg, D-91052, Erlangen, Germany
Prashant Saxena Corresponding author email<EMAIL_ADDRESS>James Watt School of Engineering, University of Glasgow, Glasgow G12 8LT,
United Kingdom

###### Abstract

Cutting-edge smart materials are transforming the domains of soft robotics,
actuators, and sensors by harnessing diverse non-mechanical stimuli, such as
electric and magnetic fields. Accurately modelling their physical behaviour
necessitates an understanding of the complex interactions between the
structural deformation and the fields in the surrounding medium. For thin
shell structures, this challenge is addressed by developing a shell model that
effectively incorporates the three-dimensional field it is embedded in by
appropriately accounting for the relevant boundary conditions. This study
presents a model for the nonlinear deformation of thin hyperelastic shells,
incorporating Kirchhoff-Love assumptions and a rigorous variational approach.
The shell theory is derived from 3D nonlinear elasticity by dimension
reduction while preserving the boundary conditions at the top and bottom
surfaces of the shell. Consequently, unlike classical shell theories, this
approach can distinguish between pressure loads applied at the top and bottom
surfaces, and delivers a platform to include multi-physics coupling. Numerical
examples are presented to illustrate the theory and provide a physical
interpretation of the novel mechanical variables of the model.

Keywords: Thin shells; Kirchhoff-Love; Nonlinear elasticity; Dimension
reduction

##  1 Introduction

Advancements in smart materials based on multi-physics coupling are
transforming the technological landscape in the field of soft robotics,
actuators, and sensors. These cutting-edge materials harness non-mechanical
stimuli from electric, magnetic, thermal, or chemical fields (Jolly et al.,
1996; McKay et al., 2010; Kim et al., 2012; Sheng et al., 2012), presenting
considerable challenges in modelling due to their complex physics.
Specifically, the interplay between the surrounding space and the soft,
deformable body is pivotal in electroelasticity, magnetoelasticity, and fluid-
structure interaction scenarios (Dorfmann and Ogden, 2014; Pelteret and
Steinmann, 2020). These problems often involve thin structures, where one
characteristic dimension is negligible compared to the others, ensuring an
appreciable deformation response to external stimuli in the surrounding free
space. Classically, thin structures are modelled as lower-dimensional
manifolds embedded in three-dimensional space, incorporating appropriate
kinematic simplifications. While these simplifications are often adequate to
model such structures in isolation, they pose some challenges when coupling
with the external ambient is required. Consequently, there is a need for a
thin shell formulation that is aware of the physical fields
(electric/magnetic/pressure) in the embedding free-space. Working towards this
aim, we present a novel thin shell theory derived from 3D nonlinear elasticity
via dimension reduction that preserves the knowledge of fields and boundary
conditions at the top and bottom surfaces. This theory can be extended to
modelling interactions between the thin shell and the free-space in
electroelastic and magnetoelastic scenarios as has been preliminarily
demonstrated in (Ghosh et al., 2023).

Modelling slender structures, such as rods, membranes, plates, and shells,
which exhibit both material and geometric nonlinearities is a particularly
relevant due to their extensive use in engineering applications such as tyres,
airbags, air springs, buffers, pneumatic actuators (Galley et al., 2019), and
soft grippers (Hao et al., 2017). When modelling physical phenomena on curved
surfaces, definitions for geometric quantities (normal vectors, curvatures,
etc.) and differential surface operators (gradients, divergence, etc.) are the
key ingredients (Steinmann, 2015). These quantities may be defined based on
either two-dimensional, curvilinear local coordinates living on the manifold
or on global coordinates of the surrounding, three-dimensional space. A
substantial body of literature exists on the modelling of nonlinear membranes.
Using semi-analytical approaches and simplifications arising from axisymmetry,
large inflation of hyperelastic membranes have been studied for circular
(Adkins and Rivlin, 1952; Hart-Smith and Crisp, 1967; Yang and Feng, 1970;
Saxena et al., 2019), cylindrical (Khayat et al., 1992; Guo, 2001; Pamplona et
al., 2006; Reddy and Saxena, 2018), spherical (Akkas, 1978; Verron et al.,
1999; Xie et al., 2016), and toroidal (Tamadapu and DasGupta, 2013; Reddy and
Saxena, 2017; Venkata and Saxena, 2019) membranes.

An overview of classical shell models with both analytical and computational
approaches is provided, for example, in (Simo and Fox, 1989; Cirak et al.,
2000; Bischoff et al., 2004; Kiendl et al., 2009), and more comprehensively in
the textbooks (Basar and Krätzig, 1985; Wempner and Talaslidis, 2002;
Blaauwendraad and Hoefakker, 2014; Radwańska et al., 2017). Various numerical
techniques have been proposed to simulate nonlinear deformation of
hyperelastic thin shells. Hughes et al. (2005) pioneered the use of NURBS
basis functions in finite elements which was applied to modelling thin shells
by Kiendl et al. (2009). Kiendl et al. (2015) further proposed an isogeometric
thin shell formulation which was extended by Liu et al. (2023) to model
instabilities in large deformation of thin shells. Additionally, Tepole et al.
(2015) and Roohbakhshan and Sauer (2017) developed isogeometric formulations
for modelling thin soft biological tissues. A complementary approach is based
on the definition of a finite number of restriction operators from the full
three-dimensional domain representing the slender structure and the lower
dimensional manifold, together with their corresponding extension operators
(Heltai and Zunino, 2023). This approach uses the original (non-linear) model
on the three-dimensional domain, while restricting the space of allowed
solutions and test functions to a smaller subspace, isomorphic to a space on
the lower dimensional manifold. This effectively reduces the dimensionality of
the problem, while naturally retaining the coupling with the surrounding
ambient (Alzetta and Heltai, 2020; Heltai et al., 2021).

In this work, a theory for modelling thin shells is derived from three-
dimensional incompressible hyperelasticity using the Kirchhoff–Love hypothesis
(Niordson, 1985). The effect of external loads is systematically captured in
the formalism as the two-dimensional system is derived through suitable
approximations. This approach proves indispensible when analyzing a
magnetoelastic shell, where careful consideration of changes in the
surrounding magnetic field in the embedding space due to shell deformation is
crucial (Ghosh et al., 2023). This contrasts with recent works on thin
magnetoelastic structures that need to impose a “weak magnetization”
assumption to simplify the problem by neglecting any perturbations in the
magnetic field in the surrounding free space (Barham et al., 2008; Reddy and
Saxena, 2017; Yan et al., 2020).

The salient features of the present work are the following:

  1. 1.

Distinction between top and bottom surface pressures: Traditional shell models
typically consider externally applied loads at the mid-surface of the
structure. In contrast, our approach distinguishes between loads applied to
the top and bottom surfaces, accurately capturing these boundary conditions.
The implications of this deviation from conventional practice are explored
through numerical examples.

  2. 2.

Comprehensive characterisation of nonlinear deformation: Nonlinear deformation
of a thin Kirchhoff–Love shell is completely characterised by the deformed
mid-surface position vector $(\bm{r})$, thickness stretch $(\lambda)$, and the
unit normal vector $(\bm{n})$ to the deformed surface. Variational
formulations of conventional shell models are directly based on $\bm{r}$ and
therefore only need to consider perturbations $\delta\bm{r}$ in the analysis.
Our approach, derived systematically from three to two dimensions, naturally
incorporates perturbations $\delta\lambda$ and $\delta\bm{n}$. This enrichment
adds complexity to the derivation of the shell system of equations and
necessitates a unique application of Green’s theorem. This study addresses
these complexities, providing a generalised system of partial differential
equations with boundary terms that encompass these effects.

  3. 3.

Surface-based partial differential equations: The partial differential
equations governing shell equilibrium derived in this work are entirely based
on the surface fields. Series expansion of all fields along the thickness
coordinate ensures all information from the top and bottom surfaces is
captured and an accuracy up to the linear order of the through-thickness
parameter is maintained. This is in contrast to conventional approaches (Cirak
and Ortiz, 2001; Kiendl et al., 2015) where stress resultants are integrated
through the thickness to formulate a virtual work principle based on effective
normal force and bending moment.

###  1.1 Outline

The structure of this paper is as follows: The preliminaries of nonlinear
elasticity are presented in Section 2. Section 3 presents the necessary
geometric definitions and the kinematics of the Kirchhoff–Love thin shell
formulation. Through a variational approach, equations of 3D elasticity are
systematically reduced to the mid-surface of the shell in Section 4. Section 5
provides the solution of two example boundary value problems to illustrate the
novel features of the derived shell equations. Section 6 concludes with
closing remarks and an outline for future research in this area. Equations
corresponding to the geometrical description of the shell are presented in
Appendix A, details about the Green’s theorem are presented in Appendix B, and
the equations corresponding to the variation of relevant physical variables
are presented in Appendix C.

###  1.2 Notation

A variable typeset in a normal weight font represents a scalar. A bold weight
font denotes a first or second-order tensor. Latin indices $i,j,k,\dots$ vary
from $1$ to $3$ while Greek indices $\alpha,\beta,\gamma,\dots$, used for
surface variable components, vary from $1$ to $2$. Einstein summation
convention is used throughout. $\mathbf{e}_{i}$ represent the basis vectors of
an orthonormal $(x_{1},x_{2},x_{3})$ coordinate system in three-dimensional
Euclidean space. The three covariant basis vectors for a surface point are
denoted as $\mathbf{a}_{i}$, where $\mathbf{a}_{\alpha}$ are the two
tangential vectors and $\mathbf{a}_{3}$ as the normal vector with
$\theta^{\alpha}$ and $\eta$ as the respective coordinates.

The comma symbol in a subscript represents partial derivative with respect to
the surface coordinates, for example, $A_{,\beta}$ is the partial derivative
of $A$ with respect to $\theta^{\beta}$. The scalar product of two vectors
$\bm{p}$ and $\bm{q}$ is denoted $\bm{p}\cdot\bm{q}$, and the dyadic product
of these vectors is a second order tensor $\bm{H}=\bm{p}\otimes\bm{q}$.
Operation of a second order tensor $\bm{H}$ on a vector $\bm{p}$ is given by
$\bm{H}\bm{p}$. The scalar product of two tensors $\bm{H}$ and $\bm{G}$ is
denoted $\bm{H}\mathrel{\mathop{\ordinarycolon}}\bm{G}$. The notation
$\mathinner{\\!\left\lVert\cdot\right\rVert}$ represents the Euclidean norm.
For a second order tensor in its component form
$\bm{H}=H^{ij}\mathbf{a}_{i}\otimes\mathbf{a}_{j}$, the matrix is denoted
$\left[H^{ij}\right]$. Circular brackets $(\,)$ are used to denote the
parameters of a function and square brackets $[\,]$ are used to group
mathematical expressions.

##  2 Non-linear elastostatics preliminaries

Consider a homogenous body occupying regions of space
${\mathcal{B}}_{\mathrm{0}}$ and $\mathcal{B}$ in ${\mathbb{R}}^{3}$ in its
reference and deformed configuration, respectively, as depicted in Figure 1.
The boundaries of the regions are $\partial{\mathcal{B}}_{\mathrm{0}}$ and
$\partial\mathcal{B}$. A point $\bm{X}\in\mathcal{B}_{0}$ is connected to a
point ${\bm{x}}\in{\mathcal{B}}$ through a bijective mapping
${\bm{\chi}}\mathrel{\mathop{\ordinarycolon}}{\mathcal{B}}_{\mathrm{0}}\rightarrow\mathcal{B}$.
The deformation gradient $\bm{F}$ is defined as
$\bm{F}=\dfrac{\partial\bm{\chi}}{\partial\bm{X}}$ with the Jacobian,
$J=\mathrm{det}\bm{F}>0$. The right Cauchy-Green tensor is then defined by
$\bm{C}={\bm{F}}^{\mathrm{T}}\bm{F}$.

The hyperelastic material model assumes that the free energy density function
per unit reference volume is of the form $\Omega=\Omega\left(\bm{F}\right).$
Objectivity and isotropy require that the free energy takes the form

$\Omega=\breve{{\Omega}}\left(\bm{C}\right)=\widetilde{{\Omega}}\left(I_{1},I_{2},I_{3}\right),$
(2.1)

where $I_{1},I_{2},I_{3}$ are scalar invariants of $\bm{C}$, given by

$I_{1}=\mathrm{tr}\bm{C},\quad
I_{2}=\dfrac{1}{2}\left[{\left[\mathrm{tr}\bm{C}\right]}^{2}-\mathrm{tr}{\bm{C}^{2}}\right],\quad\text{and}\quad
I_{3}=\mathrm{det}\bm{C}=J^{2}.$ (2.2)

Under the constraint of incompressibility ($J\equiv 1$), the energy density
function can be further simplified to

$\widetilde{{\Omega}}\left(I_{1},I_{2},I_{3}\right)=\grave{{\Omega}}\left(I_{1},I_{2}\right).$
(2.3)

For incompressible hyperelastic solids, the Piola stress tensor is given as

$\bm{P}=\dfrac{\partial{\Omega}}{\partial\bm{F}}-p{\bm{F}}^{-\mathrm{T}},$
(2.4)

where $p$ is a Lagrange multiplier due to the incompressibility constraint and
it is identified as the hydrostatic pressure.

##  3 Kirchhoff–Love thin shell

In this section, we present the necessary relations to model thin shells using
the classical Kirchhoff-Love assumption (see, e.g. Niordson, 1985).

###  3.1 Shell geometry and kinematics

$\bm{N}$
---
$\bm{n}$
---
$\bm{A}_{1}$
---
$\bm{A}_{2}$
---
$\bm{a}_{1}$
---
$\bm{a}_{2}$
---
$\mathcal{B}_{0}$
---
$\mathcal{B}$
---
$\bm{R}\left(\theta^{1},\theta^{2}\right)$
---
$\bm{r}\left(\theta^{1},\theta^{2}\right)$
---
$\left(\theta^{1},\theta^{2}\right)$
---
$\theta^{1}$
---
$\theta^{2}$
---
$\bm{e}_{1}$
---
$\bm{e}_{2}$
---
$\bm{e}_{3}$
---
$T$
---
$t$
---
O
---
O
---
O
---
O
---
Figure 1: The position vectors at the mid-surface of a shell in the two
configurations along with the local triads. The parametric coordinates
$\\{\theta^{\alpha}\\}$ are also shown.

Consider the three-dimensional body introduced in Section 2 to be a thin shell
as shown in Figure 1. Each point $\bm{X}\in{\mathcal{B}}_{\mathrm{0}}$ is
mapped from the parametric domain defined by the coordinate system
$\\{\theta^{1},\theta^{2},\eta\\}$. The Kirchhoff-Love hypothesis states that
for thin shell structures, lines perpendicular to the mid-surface of the shell
remain straight and perpendicular to the mid-surface after deformation. Hence,
assuming the shell has a thickness $T\left(\theta^{\alpha}\right)$ in the
reference configuration, the point $\bm{X}$ can be defined using a point on
the mid-surface ${S}_{\mathrm{m}}$ of the shell,
${\bm{R}}\in{S}_{\mathrm{m}}$, and the associated unit normal vector
${\bm{N}}$ by

$\bm{X}={\bm{R}}+\eta{\bm{N}},$ (3.1)

where $\eta\in[-T/{2},T/{2}]$. The same point on the mid-surface in the
deformed configuration $\bm{r}$ is related to $\bm{R}$ by

$\bm{r}=\bm{R}+\bm{u},$ (3.2)

where $\bm{u}$ denotes the mid-surface displacement vector. In a similar
fashion, each point $\bm{x}=\bm{\chi}(\bm{X})\in\mathcal{B}$ can be expressed
as

$\bm{x}=\bm{r}+\eta\bm{d},$ (3.3)

where the director $\bm{d}=\lambda\bm{n}$, $\bm{n}$ being the unit normal
vector in the deformed configuration and $\lambda$ the through-thickness
stretch defined by

$\lambda=\frac{t}{T}.$ (3.4)

Here $t(\theta^{\alpha})$ is the shell thickness after deformation. The
ensuing analysis is simplified by employing the long-wave assumption for thin
structures ($\lambda_{,\alpha}\approx 0$) that states that $\lambda$ changes
slowly along the mid-surface of the shell (Kiendl et al., 2015; Liu et al.,
2023). This assumption helps simplify the mathematical derivations in
obtaining the shell system of equations, since the shifter relations in the
deformed configuration no longer contain a term involving the parametric
derivative of $\lambda$. Mathematical details on this are provided in Appendix
A.3. For thin shells modelled with the Kirchhoff-Love assumption, transverse
shear and localised necking are unlikely to occur and therefore the long-wave
assumption is acceptable. Indeed, the long-wave assumption works well for most
problems concerning nonlinear deformation of thin shells modelled using the
Kirchhoff-Love model. However, it is recognised that it should be discarded
for problems that require modelling of localised necking and transverse shear
deformation.

Table 1 lists the surface parameters and Table 2 presents the surface and
volume elements of the shell. The expressions and associated derivations are
elaborated on in Appendix A. The boundaries,
$\partial{\mathcal{B}}_{\mathrm{0}}$ and $\partial\mathcal{B}$ can be written
as $\partial{\mathcal{B}}_{\mathrm{0}}=S_{\mathrm{t}}\cup S_{\mathrm{b}}\cup
S_{\mathrm{l}}$ and $\partial\mathcal{B}=s_{\mathrm{t}}\cup s_{\mathrm{b}}\cup
s_{\mathrm{l}}$, where the subscripts, t, b, and l, represent the top, bottom,
and lateral surfaces, respectively. Here, the top surface is the boundary that
is reached along the unit outward normal vector.

Table 1:  A list of parameters used to describe the geometry of the thin shell Surface Parameters | Reference Configuration | Deformed Configuration
---|---|---
Covariant basis vectors at the mid-surface | $\bm{A}_{\alpha}$ | $\bm{a}_{\alpha}$
Covariant metric tensor at the mid-surface | $A_{\alpha\beta}$ | $a_{\alpha\beta}$
Determinent of the covariant metric tensor at the mid-surface | $A$ | $a$
Stretch at the mid-surface | – – | $\hat{a}=a/A$
Contrvariant basis vectors at a mid-surface | $\bm{A}^{\alpha}$ | $\bm{a}^{\alpha}$
Contravariant metric tensor at the mid-surface | $A^{\alpha\beta}$ | $a^{\alpha\beta}$
Determinent of the contravariant metric tensor at the mid-surface | $A^{-1}$ | $a^{-1}$
Christoffel symbol at the mid-surface | ${\Gamma}^{\alpha}_{\beta\gamma}$ | ${\gamma}^{\alpha}_{\beta\gamma}$
Parametric derivative of the metric at the mid-surface | ${A}_{,\zeta}$ | ${a}_{,\zeta}$
Normal at the mid-surface | $\bm{N}$ | $\bm{n}$
Tangent on the bounding curve of the mid-surface | $\bm{\tau}$ | $--$
In-plane normal on the bounding curve of the mid-surface | $\bm{\nu}=\nu_{\alpha}\bm{A}^{\alpha}$ | $\overline{\bm{\nu}}=\nu_{\alpha}\bm{a}^{\alpha}$
Projection tensor at the mid-surface | $\bm{I}$ | $\bm{i}$
Curvature tensor at the mid-surface | $\bm{K}$ | $\bm{\kappa}$
Mean curvature at the mid-surface | $H$ | $h$
Gaussian curvature at the mid-surface | $K$ | $\kappa$
Covariant basis vectors at a shell-point | $\bm{G}_{\alpha}=\bm{M}\bm{A}_{\alpha}\ \text{and}\ \bm{M}=\bm{I}-\eta\bm{K}$ | $\bm{g}_{\alpha}=\bm{\mu}\bm{a}_{\alpha}\ \text{and}\ \bm{\mu}=\bm{i}-\eta\lambda\bm{\kappa}$
Covariant metric tensor at a shell-point | $G_{\alpha\beta}$ | $g_{\alpha\beta}$
Contravariant basis vectors at a shell-point | ${\bm{G}}^{\alpha}={\bm{M}}^{-\mathrm{T}}{\bm{A}}^{\alpha}$ | ${\bm{g}}^{\alpha}={\bm{\mu}}^{-\mathrm{T}}{\bm{a}}^{\alpha}$
Contravariant metric tensor at a shell-point | $G^{\alpha\beta}$ | $g^{\alpha\beta}$
Tangent at a shell-point on the lateral surface | $\bm{\tau}_{\mathrm{l}}=\bm{M}\bm{\tau}/c\ \text{and}\ c=\mathinner{\\!\left\lVert\bm{M}\bm{\tau}\right\rVert}$ | $--$
In-plane normal at a shell-point on the lateral surface | $\bm{\nu}_{\mathrm{l}}={c}^{-1}\left[\bm{I}+\eta\left[\bm{K}-2H\bm{I}\right]\right]\bm{\nu}$ | $--$

Table 2:  Surface and volume elemental quantities Volume/Surface elements | Reference Configuration | Deformed Configuration
---|---|---
Area for the convected coordinates | $dP$ | $dP$
Elemental area at the curved mid-surface | $dS_{\mathrm{m}}$ | $ds_{\mathrm{m}}$
Area at a shell-point | $dS=Md{S}_{\mathrm{m}},\,M=\mathrm{det}\bm{M}$ | $ds=\mu{\hat{a}}^{1/2}dS_{\mathrm{m}},\ \mu=\mathrm{det}\bm{\mu}$
Area at the top surface | $dS_{\mathrm{t}}$ | $ds_{\mathrm{t}}$
Area at the bottom surface | $dS_{\mathrm{b}}$ | $ds_{\mathrm{b}}$
Area on the lateral surface | $dS_{\mathrm{l}}=c\,dl\,d\eta$ | $--$
Volume at a shell-point | $dV=dS\,d\eta$ | $--$

The deformation gradient can be expanded so as to separate the thickness
variable $\eta$ from the surface parameters.

$\displaystyle\bm{F}={F}^{\alpha}_{0_{\beta}}\
\bm{a}_{\alpha}\otimes{\bm{A}}^{\beta}+\eta{F}^{\alpha}_{1_{\beta}}\
\bm{a}_{\alpha}\otimes{\bm{A}}^{\beta}+\eta^{2}{F}^{\alpha}_{2_{\beta}}\
\bm{a}_{\alpha}\otimes{\bm{A}}^{\beta}+\bm{\mathcal{O}}(\eta^{3})+\lambda\bm{n}\otimes\bm{N},$
(3.5)

where the individual scalar components are written using the geometric
parameters provided in Table 1 as

$\displaystyle{F}^{\alpha}_{0_{\beta}}=\delta^{\alpha}_{\beta},\quad{F}^{\alpha}_{1_{\beta}}=-\lambda{b}_{\beta}^{\
\alpha}+{B}_{\beta}^{\
\alpha},\quad\text{and}\quad{F}^{\alpha}_{2_{\beta}}={B}_{\delta}^{\
\alpha}{B}_{\beta}^{\ \delta}-\lambda{b}_{\delta}^{\ \alpha}{B}_{\beta}^{\
\delta}.$ (3.6)

Here ${B}_{\beta}^{\ \alpha}$ and ${b}_{\beta}^{\ \alpha}$ are the components
of the curvature tensors $\bm{K}$ and $\bm{\kappa}$, respectively. Applying
the same decomposition to the inverse of the deformation gradient tensor, we
obtain

$\displaystyle\bm{F}^{-1}={F}^{-1\alpha}_{0_{\beta}}\
{\bm{A}}_{\alpha}\otimes\bm{a}^{\beta}+\eta{F}^{-1\alpha}_{1_{\beta}}\
{\bm{A}}_{\alpha}\otimes\bm{a}^{\beta}+\eta^{2}{F}^{-1\alpha}_{2_{\beta}}\
{\bm{A}}_{\alpha}\otimes\bm{a}^{\beta}+\bm{\mathcal{O}}(\eta^{3})+\frac{1}{\lambda}\bm{N}\otimes\bm{n},$
(3.7)

with

${F}^{-1\alpha}_{0_{\beta}}=\delta^{\alpha}_{\beta},\quad{F}^{-1\alpha}_{1_{\beta}}=\lambda{b}_{\beta}^{\
\alpha}-{B}_{\beta}^{\
\alpha},\quad\text{and}\quad{F}^{-1\alpha}_{2_{\beta}}=\lambda^{2}{b}_{\delta}^{\
\alpha}{b}_{\beta}^{\ \delta}-\lambda{B}_{\delta}^{\ \alpha}{b}_{\beta}^{\
\delta}.$ (3.8)

Upon recognising that $\bm{a}_{\alpha}\cdot\bm{n}=0$, the right Cauchy–Green
tensor $\bm{C}$ can be expanded as

$\bm{C}={C}_{0_{\alpha\beta}}{\bm{A}}^{\alpha}\otimes{\bm{A}}^{\beta}+\eta{C}_{1_{\alpha\beta}}{\bm{A}}^{\alpha}\otimes{\bm{A}}^{\beta}+\eta^{2}{C}_{2_{\alpha\beta}}{\bm{A}}^{\alpha}\otimes{\bm{A}}^{\beta}+\bm{\mathcal{O}}(\eta^{3})+\lambda^{2}\bm{N}\otimes\bm{N},$
(3.9)

where the individual components are

$\displaystyle{C}_{0_{\alpha\beta}}$ $\displaystyle=$ $\displaystyle
a_{\alpha\beta},$ $\displaystyle{C}_{1_{\alpha\beta}}$ $\displaystyle=$
$\displaystyle-\lambda{b}_{\beta}^{\
\gamma}a_{\gamma\alpha}-\lambda{b}_{\alpha}^{\
\gamma}a_{\gamma\beta}+{B}_{\alpha}^{\ \gamma}a_{\gamma\beta}+{B}_{\beta}^{\
\gamma}a_{\gamma\alpha},$ $\displaystyle{C}_{2_{\alpha\beta}}$
$\displaystyle=$ $\displaystyle{B}_{\alpha}^{\ \gamma}{B}_{\beta}^{\
\delta}a_{\delta\gamma}+\lambda^{2}{b}_{\alpha}^{\ \delta}{b}_{\beta}^{\
\gamma}a_{\delta\gamma}+{B}_{\gamma}^{\ \delta}{B}_{\alpha}^{\
\gamma}a_{\delta\beta}+{B}_{\gamma}^{\ \delta}{B}_{\beta}^{\
\gamma}a_{\delta\alpha}-\lambda{b}_{\gamma}^{\ \delta}{B}_{\beta}^{\
\gamma}a_{\delta\alpha}-\lambda{b}_{\gamma}^{\ \delta}{B}_{\alpha}^{\
\gamma}a_{\delta\beta}$ (3.10) $\displaystyle-\lambda{B}_{\beta}^{\
\gamma}{b}_{\alpha}^{\ \delta}a_{\delta\gamma}-\lambda{B}_{\alpha}^{\
\gamma}{b}_{\beta}^{\ \delta}a_{\gamma\delta}.$

The decompositions in equations allow us to write the relevant quantities at
the mid-surface for later use. The mid-surface right Cauchy-Green tensor
$\bm{C}_{\mathrm{m}}$ is given by

$\bm{C}_{\mathrm{m}}=\bm{C}\Big{|}_{\eta=0}={C}_{0_{\alpha\beta}}{\bm{A}}^{\alpha}\otimes{\bm{A}}^{\beta}+\lambda^{2}\bm{N}\otimes\bm{N}.$
(3.11)

Its determinant squared is expressed as

$J_{0}^{2}=\mathrm{det}\bm{C}_{\mathrm{m}}=\mathrm{det}\left[{C_{\mathrm{m}_{j}}}^{i}\right]=\mathrm{det}\left[{C}_{\mathrm{m}_{jk}}{A}^{ki}\right]=\frac{\mathrm{det}\left[{C}_{\mathrm{m}_{ij}}\right]}{\mathrm{det}\left[{A}_{ij}\right]}=\frac{\mathrm{det}\left[a_{\alpha\beta}\right]\lambda^{2}}{\mathrm{det}\left[{A}_{\alpha\beta}\right]}=\frac{a\lambda^{2}}{A},$
(3.12)

where ${A}_{ij}={\bm{A}}_{i}\cdot{\bm{A}}_{j}$ and
${A}^{ij}={\bm{A}}^{i}\cdot{\bm{A}}^{j}$ are the components of the three-
dimensional covariant and contravariant metric tensors , respectively, on the
mid-surface with ${\bm{A}}_{3}={\bm{A}}^{3}={\bm{N}}$. The incompressibility
constraint $J_{0}=1$ implies

$\lambda={\hat{a}}^{-1/2},$ (3.13)

where $\hat{a}=a/A$.

###  3.2 Divergence of stress and Green’s theorem

The divergence of the Piola stress tensor is required for derivation of the
governing equations for the Kirchhoff-Love thin shell in later sections. It is
given by

$\displaystyle\mathrm{Div}\bm{P}$ $\displaystyle=$
$\displaystyle\bm{P}_{,\alpha}{\bm{G}}^{\alpha}+\frac{\partial\bm{P}}{\partial\eta}{\bm{N}},$
(3.14) $\displaystyle=$
$\displaystyle\bm{A}_{0}+\eta\bm{A}_{1}+\eta^{2}\bm{A}_{2}+\bm{\mathcal{O}}(\eta^{3}),$

where
${\bm{G}}^{\alpha}={\bm{M}}^{-\mathrm{T}}{\bm{A}}^{\alpha}+\bm{\mathcal{O}}(\eta^{2})$.
$\bm{A}_{0},\bm{A}_{1},\bm{A}_{2}$ are the zeroth; first; and second; order
terms independent of the thickness variable $\eta$ and are given by

$\displaystyle\bm{A}_{0}$ $\displaystyle=$
$\displaystyle{\bm{P}_{0}}_{,\alpha}{\bm{A}}^{\alpha}+\bm{P}_{1}\bm{N},$
$\displaystyle\bm{A}_{1}$ $\displaystyle=$ $\displaystyle{B}_{\delta}^{\
\alpha}{\bm{P}_{0}}_{,\alpha}{\bm{A}}^{\delta}+{\bm{P}_{1}}_{,\alpha}{\bm{A}}^{\alpha}+\bm{P}_{2}\bm{N},$
$\displaystyle\bm{A}_{2}$ $\displaystyle=$ $\displaystyle{B}_{\zeta}^{\
\alpha}{B}_{\delta}^{\
\zeta}{\bm{P}_{0}}_{,\alpha}{\bm{A}}^{\delta}+{B}_{\delta}^{\
\alpha}{\bm{P}_{1}}_{,\alpha}{\bm{A}}^{\delta}+\frac{1}{2}{\bm{P}_{2}}_{,\alpha}{\bm{A}}^{\alpha}+\frac{1}{2}\bm{P}_{3}\bm{N},$
(3.15)

with

$\displaystyle\bm{P}$ $\displaystyle=$
$\displaystyle\bm{P}_{0}+\eta\bm{P}_{1}+\frac{\eta^{2}}{2}\bm{P}_{2}+\frac{\eta^{3}}{6}\bm{P}_{3}+\bm{\mathcal{O}}(\eta^{4}).$
(3.16)

The stress at the top and bottom boundaries can be obtained by setting
$\eta=T/2$ and $-T/2$, respectively. We will require the zeroth. and the
first. order terms of the Piola stress explicitly given as

$\bm{P}_{0}=\bm{\mathcal{P}}_{0}-p_{0}\bm{F}^{-\mathrm{T}}_{0},\quad\text{and}\quad\bm{P}_{1}=\bm{\mathcal{P}}_{1}-p_{0}\bm{F}^{-\mathrm{T}}_{1}-p_{1}\bm{F}^{-\mathrm{T}}_{0},$
(3.17)

with

$\frac{\partial\Omega}{\partial\bm{F}}=\bm{\mathcal{P}}(\eta,\theta^{\alpha})=\bm{\mathcal{P}}_{0}+\eta\bm{\mathcal{P}}_{1}+\bm{\mathcal{O}}(\eta^{2}),\quad\text{and}\quad\bm{F^{-\mathrm{T}}}=\bm{F^{-\mathrm{T}}}_{0}+\eta\bm{F^{-\mathrm{T}}}_{1}+\bm{\mathcal{O}}(\eta^{2}).$
(3.18)

Mathematical analysis for shell mechanics often requires change of integration
from the parametric domain to the shell boundary with expressions that involve
parametric derivatives (see, for example, Section 4.2). This requires a
modified form of Green’s theorem which states that for a vector with
components ${T}^{\alpha}$ the following relation holds

$\int\limits_{P}\left[{A}^{1/2}{T}^{\alpha}\right]_{,\alpha}dP=\int\limits_{\mathcal{C}_{\mathrm{m}}}{T}^{\alpha}{\nu}_{\alpha}dl.$
(3.19)

Here $\mathcal{C}_{\mathrm{m}}$ is the boundary of the curved mid-surface
$S_{\mathrm{m}}$ . A complete derivation of this result is provided in
Appendix A.

##  4 Derivation of shell equations from 3D elasticity

We now systematically derive the necessary shell equations starting from a
classical 3D potential energy functional for a hyperelastic body.

###  4.1 Variational formulation in three dimensions

The total potential energy of an incompressible elastic body undergoing a
deformation $\bm{\chi}$ is given by

$\displaystyle\Pi[\bm{\chi},p]=\int\limits_{{\mathcal{B}}_{0}}\Omega\left(\bm{F}\right)dV-\int\limits_{\mathcal{B}_{0}}p[J-1]dV-\int\limits_{\mathcal{B}_{0}}\bm{\mathfrak{B}}\cdot\bm{\chi}dV-\int\limits_{{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}^{\mathrm{u}}_{\mathrm{m}}}\bm{t}_{\ell}\cdot\bm{\chi}dl-\int\limits_{s_{\mathrm{t}}}\bm{p}_{\mathrm{t}}\cdot\bm{\chi}ds_{\mathrm{t}}-\int\limits_{s_{\mathrm{b}}}\bm{p}_{\mathrm{b}}\cdot\bm{\chi}ds_{\mathrm{b}},$
(4.1)

where $\left(\bm{\chi},p\right)$ is the set of generalised solution variables.
The body force field per unit volume is $\bm{\mathfrak{B}}$ and
$\bm{t}_{\ell}$ is the applied traction due to dead load at
${C}_{\mathrm{m}}$. The Dirichlet part of the boundary is denoted by
${C}_{\mathrm{m}}^{\mathrm{u}}$. We defined
${p}_{\mathrm{t}}\left(\theta^{\alpha}\right)$ and
${p}_{\mathrm{b}}\left(\theta^{\alpha}\right)$ as the magnitudes of external
pressure at the top and bottom surfaces of the shell, respectively, such that

$\bm{p}_{\mathrm{t}}=-p_{\mathrm{t}}\bm{n}_{\mathrm{t}},\quad\text{and}\quad\bm{p}_{\mathrm{b}}=-p_{\mathrm{b}}\bm{n}_{\mathrm{b}},$
(4.2)

with $\bm{n}_{\mathrm{t}}=\bm{n}$ and $\bm{n}_{\mathrm{b}}=-\bm{n}$ as the
outward unit normals at these two surfaces in the current configuration. Upon
arbitrary variations $\delta\bm{\chi}$ and $\delta p$ of the solution
variables, the first variation of the energy functional becomes

$\displaystyle\delta\Pi[\bm{\chi},p;\delta\bm{\chi},\delta p]$
$\displaystyle=$
$\displaystyle\underbrace{\delta\left[\int\limits_{{\mathcal{B}}_{0}}\Omega\left(\bm{F}\right)dV\right]-\int\limits_{{\mathcal{B}}_{0}}p\delta
JdV}_{I}-\int\limits_{{\mathcal{B}}_{0}}\delta p[J-1]dV$ (4.3)
$\displaystyle-\int\limits_{{\mathcal{B}}_{0}}\bm{\mathfrak{B}}\cdot\delta\bm{\chi}dV-\int\limits_{{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}_{\mathrm{m}}^{\mathrm{u}}}\bm{t}_{\ell}\cdot\delta\bm{\chi}dl-\int\limits_{s_{\mathrm{t}}}\bm{p}_{\mathrm{t}}\cdot\delta\bm{\chi}ds_{\mathrm{t}}-\int\limits_{s_{\mathrm{b}}}\bm{p}_{\mathrm{b}}\cdot\delta\bm{\chi}ds_{\mathrm{b}}.$

Upon using Equation (2.4) and an application of the divergence theorem, the
first and second terms of the above Equation (4.3) can be written as

$\displaystyle
I=-\int\limits_{{\mathcal{B}}_{0}}\mathrm{Div}\bm{P}\cdot\delta\bm{\chi}dV+\int\limits_{S_{\mathrm{t}}}\bm{P}{\bm{N}}_{\mathrm{t}}\cdot\delta\bm{\chi}dS_{\mathrm{t}}+\int\limits_{S_{\mathrm{b}}}\bm{P}{\bm{N}}_{\mathrm{b}}\cdot\delta\bm{\chi}dS_{\mathrm{b}}+\int\limits_{S_{\mathrm{l}}}\bm{P}{\bm{\nu}}_{\mathrm{l}}\cdot\delta\bm{\chi}dS_{\mathrm{l}},$
(4.4)

The terms involving the Piola stress tensor are addressed further with in
Section 4.2.1. Note that, ${\bm{N}}_{\mathrm{t}}=\bm{N}$ and
${\bm{N}}_{\mathrm{b}}=-\bm{N}$ are the unit outward normals at the top and
bottom surfaces of the shell in the reference configuration, respectively. The
remaining terms in Equation (4.3), that is, the virtual work done by the dead
load traction, body force and pressures are addressed in Section 4.2.2, where
their contributions to a modified variational form for a Kirchoff-Love thin
shell are discussed.

###  4.2 Derivation of the shell model

The subsequent discussion introduces the key steps in deriving the Kirchhoff-
Love shell equations. It will shown that the derived equations can achieve an
accuracy of $\mathcal{O}\left(T\right)$. The solution variables are
alternatively expressed as quantities defined over the mid-surface of the
shell, $\bm{r}$ and $p_{0}$, such that

$\delta\Pi[\bm{\chi},p;\delta\bm{\chi},\delta
p]=\delta\widetilde{\Pi}[\bm{r},p_{0};\delta\bm{r},\delta p_{0}].$ (4.5)

In the following sections, all integrals in Equation (4.3) are systematically
converted from a dependence on $(\bm{\chi},p)$ to a dependence on
$(\bm{r},p_{0})$. This allows for the determination of the governing equations
over the mid-surface of the shell.

Remark: Shell models are usually simplified by incorporating the plane-stress
assumption due to a shell being a thin structure. In the case of shells made
of incompressible hyperelastic materials, application of the plane-stress
condition leads to the explicit determination of the Lagrange multiplier as a
function of deformation variables (Kiendl et al., 2015; Liu et al., 2023).
Consequently, the Lagrange multiplier within the equations changes in nature
from being an independent field to a constitutive expression. In the
subsequent derivations, the Lagrange multiplier is maintained as an
independent field until the plane-stress condition is applied in Section 4.3.

####  4.2.1 Terms involving the Piola stress tensor ($I$)

Upon using Equations (A.43) and (A.44), where the latter relates the
undeformed elemental area at a shell-point to the undeformed elemental area on
the midsurface, the first term of $I$ in Equation (4.4) can be written as

$\displaystyle-\int\limits_{{\mathcal{B}}_{0}}\mathrm{Div}\bm{P}\cdot\delta\bm{\chi}dV$
$\displaystyle=-\int\limits_{S}\int\limits_{\eta}\mathrm{Div}\bm{P}\cdot\delta\bm{\chi}d\eta
dS=-\int\limits_{S_{\mathrm{m}}}\int\limits_{\eta}\mathrm{Div}\bm{P}\cdot\delta\bm{\chi}Md\eta
d{\mathcal{S}}_{\mathrm{m}},$ (4.6)
$\displaystyle\quad\text{substituting}\,\bm{\chi}\,\text{from Equation
\eqref{eq:def_map_body} and }\,\mathrm{Div}\bm{P}\,\text{from Equation
\eqref{eq:div_total_stress} gives}$
$\displaystyle=-\int\limits_{S_{\mathrm{m}}}\int\limits_{\eta}\left[\bm{A}_{0}\cdot\delta\bm{r}+\eta\delta\lambda\bm{A}_{0}\cdot\bm{n}+\eta\lambda\bm{A}_{0}\cdot\delta\bm{n}+\eta\bm{A}_{1}\cdot\delta\bm{r}+\mathcal{O}(\eta^{2})\right]Md\eta
dS_{\mathrm{m}},$ (4.7) $\displaystyle\quad\text{substituting}\,M\,\text{from
Equation \eqref{eq:shift_ref} gives}$
$\displaystyle=-\int\limits_{S_{\mathrm{m}}}\int\limits_{\eta}\left[\bm{A}_{0}\cdot\delta\bm{r}+\eta\delta\lambda\bm{A}_{0}\cdot\bm{n}+\eta\lambda\bm{A}_{0}\cdot\delta\bm{n}+\eta\bm{A}_{1}\cdot\delta\bm{r}-2\eta
H\bm{A}_{0}\cdot\delta\bm{r}+\mathcal{O}\left(\eta^{2}\right)\right]d\eta
dS_{\mathrm{m}}.$ (4.8)

As $\eta\in[-T/{2},T/{2}]$ and we note that across the thickness of the shell

$\displaystyle{\int\limits_{\eta}d\eta}=T\quad\text{and}\quad{\int\limits_{\eta}\eta
d\eta}=0.$ (4.9)

This allows one to simplify Equation (4.8) as

$-\int\limits_{{\mathcal{B}}_{0}}\mathrm{Div}\bm{P}\cdot\delta\bm{\chi}dV=\int\limits_{S_{\mathrm{m}}}\left[-T\big{[}\widehat{\text{Div}}{\bm{P}_{0}}+\bm{T}_{1}\big{]}\cdot\delta\bm{r}+\mathcal{O}\left(T^{3}\right)\right]dS_{\mathrm{m}},$
(4.10)

where the mechanical traction is given by
$\bm{T}=\bm{P}\bm{N}=\bm{T}_{0}+\eta\bm{T}_{1}$ and the surface divergence is
defined as $\widehat{\text{Div}}\bm{W}=\bm{W}_{,\beta}\cdot\bm{A}^{\beta}$.

Integrals over the top and bottom boundaries in $I$:

The boundary terms in Equation (4.4) are now addressed. The second term in
Equation (4.4) over the top surface ($\eta=T/2,\bm{N}_{\mathrm{t}}=\bm{N}$)
can be written using Equation (3.16) as

$\displaystyle\int\limits_{S_{\mathrm{t}}}\bm{P}{\bm{N}}_{\mathrm{t}}\cdot\delta\bm{\chi}dS_{\mathrm{t}}$
$\displaystyle=\int\limits_{S_{\mathrm{t}}}\bm{P}\Big{|}_{\eta=T/2}{\bm{N}}_{\mathrm{t}}\cdot\delta\bm{\chi}\Big{|}_{\eta=T/2}dS_{\mathrm{t}},$
(4.11)
$\displaystyle=\int\limits_{S_{\mathrm{t}}}\left[\bm{T}_{0}\cdot\delta\bm{r}+\frac{T}{2}\delta\lambda\bm{T}_{0}\cdot\bm{n}+\frac{T}{2}\lambda\bm{T}_{0}\cdot\delta\bm{n}+\frac{T}{2}\bm{T}_{1}\cdot\delta\bm{r}+\mathcal{O}\left(T^{2}\right)\right]dS_{\mathrm{t}}.$
(4.12)

The above equation is evaluated over the top surface of the shell and should
be transformed such that the domain of integration changes to the mid-surface.
This is achieved using the shifter relation in Equation (A.46), giving

$\int\limits_{S_{\mathrm{t}}}\bm{T}_{0}\cdot\delta\bm{r}\,dS_{t}=\int\limits_{S_{\mathrm{m}}}\bm{T}_{0}\cdot\delta\bm{r}MdS_{\mathrm{m}}=\int\limits_{S_{\mathrm{m}}}\big{[}\bm{T}_{0}\cdot\delta\bm{r}-TH\bm{T}_{0}\cdot\delta\bm{r}+\mathcal{O}\left(T^{2}\right)\big{]}dS_{\mathrm{m}},$
(4.13)

and applying the same procedure to the second, third, and fourth terms in
Equation (4.12) retaining only terms of $\mathcal{O}(T)$ gives

$\int\limits_{S_{\mathrm{t}}}(\bullet)dS_{\mathrm{t}}=\int\limits_{S_{\mathrm{m}}}(\bullet)MdS_{\mathrm{m}}=\int\limits_{S_{\mathrm{m}}}\big{[}(\bullet)+\mathcal{O}\left(T^{2}\right)\big{]}dS_{\mathrm{m}}.$
(4.14)

Now that the integrals have been transformed to the mid-surface, it is
necessary to convert the variations of all the quantities
$(\delta\lambda,\delta\bm{n})$ to $\delta\bm{r}$ in order to derive the
Euler–Lagrange equations. This is achieved by transferring the integrals to
the parametric domain, applying the form of Green’s theorem stated in Equation
(3.19) and then transforming back to the shell domain. The second term of
Equation (4.12) with $\delta\lambda$ over $S_{\mathrm{m}}$ is transformed to
the parametric domain and written as

$\displaystyle\int\limits_{S_{\mathrm{m}}}T$
$\displaystyle\delta\lambda\bm{T}_{0}\cdot\bm{n}dS_{\mathrm{m}}=\int\limits_{P}T\delta\lambda\bm{T}_{0}\cdot\bm{n}{A}^{1/2}dP,$
(4.15) $\displaystyle\quad\text{substituting}\,\delta\lambda\,\text{using
Equation \eqref{eq:var_lambda}}$ (4.16)
$\displaystyle=-\int\limits_{P}T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\cdot\delta{\bm{a}}_{\alpha}{A}^{1/2}dP,$
(4.17)
$\displaystyle\quad\text{substituting}\,\delta\bm{a}_{\alpha}\,\text{from
Equation \eqref{eq: delta a lpha}}$
$\displaystyle=-\int\limits_{P}\left[{C}^{\alpha}{A}^{1/2}\right]_{,\alpha}dP+\int\limits_{S_{\mathrm{m}}}\left[T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\right]_{,\alpha}\cdot\delta\bm{r}dS_{\mathrm{m}}+\int\limits_{P}T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\left[A^{1/2}\right]_{,\alpha}\cdot\delta\bm{r}dP,$
(4.18) transforming back from parametric domain to shell domain
$\displaystyle=-\int\limits_{P}\left[{C}^{\alpha}{{A}}^{1/2}\right]_{,\alpha}dP+\int\limits_{S_{\mathrm{m}}}\left[T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\right]_{,\alpha}\cdot\delta\bm{r}dS_{\mathrm{m}}+\int\limits_{S_{\mathrm{m}}}T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}{\Gamma}^{\beta}_{\beta\alpha}\cdot\delta\bm{r}dS_{\mathrm{m}},$
(4.19)
$\displaystyle=-\int\limits_{P}\left[{C}^{\alpha}{{A}}^{1/2}\right]_{,\alpha}dP+\int\limits_{S_{\mathrm{m}}}\widehat{\text{Div}}\left(T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)\cdot\delta\bm{r}dS_{\mathrm{m}}.$
(4.20)

The first term containing the derivative over the parametric domain is
transformed using the modified Green’s theorem (3.19) as

$-\int\limits_{P}\left[{C}^{\alpha}{A}^{1/2}\right]_{,\alpha}dP=-\int\limits_{C_{\mathrm{m}}}{C}^{\alpha}{\nu}_{\alpha}dl=-\int\limits_{C_{\mathrm{m}}}T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\overline{\bm{\nu}}\cdot\delta\bm{r}dl=-\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\overline{\bm{\nu}}\cdot\delta\bm{r}dl.$
(4.21)

The third term of Equation (4.12) with $\delta\bm{n}$ is transformed by using
Equation (C.6) as

$\displaystyle\int\limits_{S_{\mathrm{m}}}T\lambda\bm{T}_{0}\cdot\delta\bm{n}dS_{\mathrm{m}}$
$\displaystyle=\int\limits_{P}T\lambda\bm{T}_{0}\cdot\delta\bm{n}{A}^{1/2}dP=-\int\limits_{P}T\lambda\big{[}\bm{T}_{0}\cdot\bm{a}^{\alpha}\big{]}\big{[}\bm{n}\cdot\delta\bm{a}_{\alpha}\big{]}{A}^{1/2}dP=0.$
(4.22)

The above term vanishes on account of the Kirchhoff-Love assumption of zero
shear stresses, i.e., $\bm{T}_{0}\cdot\bm{a}^{\alpha}=0$.

Application of similar mathematical manipulations to the second and third
terms in Equation (4.4) corresponding to the top and bottom boundaries and
neglecting the higher order terms results in

$\displaystyle\int\limits_{S_{\mathrm{t}}}\bm{P}{\bm{N}}_{\mathrm{t}}\cdot\delta\bm{\chi}dS_{\mathrm{t}}$
$\displaystyle+\int\limits_{S_{\mathrm{b}}}\bm{P}{\bm{N}}_{\mathrm{b}}\cdot\delta\bm{\chi}dS_{\mathrm{b}}$
$\displaystyle=$
$\displaystyle\int\limits_{{\mathcal{S}}_{\mathrm{m}}}\Bigg{[}\widehat{\text{Div}}\left(T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)+T\bm{T}_{1}-2TH\bm{T}_{0}\Bigg{]}\cdot\delta\bm{r}d{\mathcal{S}}_{\mathrm{m}}$
$\displaystyle-\int\limits_{{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}_{\mathrm{m}}^{\mathrm{u}}}T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\overline{\bm{\nu}}\cdot\delta\bm{r}dl.$
(4.23)

Integral over the lateral boundary in $I$:

The fourth term in Equation (4.4) is the contribution of Piola stress over the
lateral boundary. This integral can be expressed using Equations (4.9),
(A.62), and (A.65) as

$\int\limits_{S_{\mathrm{l}}}\bm{P}\bm{\nu}_{\text{l}}\cdot\delta\bm{\chi}dS_{\mathrm{l}}=\int\limits_{C_{\mathrm{m}}}\left[T\bm{P}_{0}\bm{\nu}\cdot\delta\bm{r}+\mathcal{O}\left(T^{3}\right)\right]dl=\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}\left[T\bm{P}_{0}\bm{\nu}\cdot\delta\bm{r}+\mathcal{O}\left(T^{3}\right)\right]dl.$
(4.24)

Total contribution:

Upon combining Equations (4.10), (4.23) and (4.24), and neglecting the higher
order terms in $T$, the entire contribution from the Piola-stress from
Equation (4.4) is given by

$\displaystyle
I=\int\limits_{{\mathcal{S}}_{\mathrm{m}}}\Bigg{[}-T\,\widehat{\text{Div}}\bm{P}_{0}+\widehat{\text{Div}}\left(T\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)-2TH\bm{T}_{0}\Bigg{]}\cdot\delta\bm{r}d{\mathcal{S}}_{\mathrm{m}}+\int\limits_{{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}_{\mathrm{m}}^{\mathrm{u}}}T\Big{[}{\bm{P}_{0}}\bm{\nu}-\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{\nu}\Big{]}\cdot\delta\bm{r}dl$
(4.25)

####  4.2.2 Terms involving the external loads

The terms involving the external pressure, traction, and body force in
Equation (4.3) involve the variation $\delta\bm{\chi}$. In order to derive the
shell equations, they need to be expressed in terms of the variation
$\delta\bm{r}$. It is noted that one of the main novelties in the present work
is the differentiation between the pressure load applied at the top and bottom
surfaces. Since the pressure is always applied in the current configuration,
the integrals need to be transformed to the reference and the parametric
domain for application of the Green’s theorem.

Consider the second to last term in Equation (4.3) involving pressure on the
top surface ($\eta=T/2,\bm{n}_{t}=\bm{n}$). It is written using Equations
(4.2), (A.52), and (A.53) as

$\displaystyle\int\limits_{s_{\mathrm{t}}}\bm{p}_{\mathrm{t}}\cdot\delta\bm{\chi}ds_{\mathrm{t}}=$
$\displaystyle-\int\limits_{S_{\mathrm{t}}}p_{\mathrm{t}}\bm{n}\cdot\delta\bm{\chi}\Big{|}_{\eta=T/2}\mu\Big{|}_{\eta=T/2}\frac{{\widehat{a}}^{1/2}}{M}dS_{\mathrm{t}}$
(4.26) using Equations (A.46) and (C.5) $\displaystyle=$
$\displaystyle-\int\limits_{S_{\mathrm{t}}}p_{\mathrm{t}}\bm{n}\cdot\left[\delta\bm{r}+\delta\lambda\frac{T}{2}\bm{n}+\frac{T}{2}\lambda\delta\bm{n}\right]\left[1-T\left[\lambda
h+H\right]+\mathcal{O}\left(T^{2}\right)\right]{\widehat{a}}^{1/2}dS_{\mathrm{t}},$
(4.27)
$\displaystyle\quad\text{using}\;\bm{n}\cdot\bm{n}=1,\;\bm{n}\cdot\delta\bm{n}=0,\;\lambda\widehat{a}^{1/2}=1$
$\displaystyle=$
$\displaystyle\int\limits_{S_{\mathrm{t}}}\left[-{\lambda}^{-1}p_{\mathrm{t}}\bm{n}\cdot\delta\bm{r}+\bm{a}^{\alpha}\cdot\delta\bm{a}_{\alpha}\frac{T}{2}p_{\mathrm{t}}+T\left[h+\lambda^{-1}H\right]p_{\mathrm{t}}\bm{n}\cdot\delta\bm{r}+\mathcal{O}\left(T^{2}\right)\right]dS_{\mathrm{t}}.$
(4.28)

Upon defining the average pressure
$\overline{p}=[p_{\mathrm{t}}+p_{\mathrm{b}}]/2$, the total contribution of
pressure to the top and bottom surfaces in Equation (4.3) can be written as

$\displaystyle\int\limits_{s_{\mathrm{t}}}\bm{p}_{\mathrm{t}}\cdot\delta\bm{\chi}ds_{\mathrm{t}}+\int\limits_{s_{\mathrm{b}}}\bm{p}_{\mathrm{b}}\cdot\delta\bm{\chi}ds_{\mathrm{b}}$
$\displaystyle=$
$\displaystyle\int\limits_{\mathcal{S}_{\mathrm{t}}}\left[-{\lambda}^{-1}p_{t}\bm{n}\cdot\delta\bm{r}+T\left[h+\lambda^{-1}H\right]p_{t}\bm{n}\cdot\delta\bm{r}+\mathcal{O}\left(T^{2}\right)\right]dS_{\mathrm{t}}$
(4.29) $\displaystyle+$
$\displaystyle\int\limits_{S_{\mathrm{b}}}\left[{\lambda}^{-1}p_{\mathrm{b}}\bm{n}\cdot\delta\bm{r}+T\left[h+\lambda^{-1}H\right]p_{\mathrm{b}}\bm{n}\cdot\delta\bm{r}+\mathcal{O}\left(T^{2}\right)\right]dS_{\mathrm{b}}$
$\displaystyle+$
$\displaystyle\int\limits_{S_{\mathrm{m}}}\left[\bm{a}^{\alpha}\cdot\delta\bm{a}_{\alpha}T\overline{p}+\mathcal{O}\left(T^{2}\right)\right]dS_{\mathrm{m}}.$

The only term above that needs further modification is the one with
$\delta\bm{a}_{\alpha}$. This is achieved by pulling it back to the parametric
domain and applying the Green’s theorem to get

$\displaystyle\int\limits_{S_{\mathrm{m}}}\bm{a}_{\alpha}\cdot\delta\bm{a}^{\alpha}T\overline{p}dS_{\mathrm{m}}$
$\displaystyle=$
$\displaystyle\int\limits_{P}T\overline{p}\bm{a}^{\alpha}\cdot\delta\bm{a}_{\alpha}{A}^{1/2}dP,$
(4.30) $\displaystyle=$
$\displaystyle\int\limits_{P}\left[{{E}}^{\alpha}A^{1/2}\right]_{,\alpha}dP-\int\limits_{S_{\mathrm{m}}}\Big{[}\left[T\overline{p}\bm{a}^{\alpha}\right]_{,\alpha}+T\overline{p}\bm{a}^{\alpha}{\Gamma}^{\beta}_{\beta\alpha}\Big{]}\cdot\delta\bm{r}dS_{\mathrm{m}},$
$\displaystyle=$
$\displaystyle\int\limits_{P}\left[{{E}}^{\alpha}A^{1/2}\right]_{,\alpha}dP-\int\limits_{S_{\mathrm{m}}}\widehat{\text{Div}}\left(T\overline{p}\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)\cdot\delta\bm{r}dS_{\mathrm{m}}$

with ${E}^{\alpha}=T\overline{p}\bm{a}^{\alpha}\cdot\delta\bm{r}$. An
application of the modified Green’s theorem (3.19) to the first term results
in

$\int\limits_{P}\left[{{E}}^{\alpha}A^{1/2}\right]_{,\alpha}dP=\int\limits_{C_{\mathrm{m}}}{E}^{\alpha}{\nu}_{\alpha}dl=\int\limits_{C_{\mathrm{m}}}T\overline{p}\overline{\bm{\nu}}\cdot\delta\bm{r}dl=\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}T\overline{p}\overline{\bm{\nu}}\cdot\delta\bm{r}dl.$
(4.31)

Upon substituting the Equations (4.30) and (4.31) back into Equation (4.29),
using the shifter relations (A.52) to transform all integrals to the mid-
surface, and neglecting the $\mathcal{O}(T^{2})$ terms, the total contribution
due to the pressure loads can be written as

$\displaystyle\int\limits_{s_{\mathrm{t}}}\bm{p}_{\mathrm{t}}\cdot\delta\bm{\chi}ds_{\mathrm{t}}+\int\limits_{s_{\mathrm{b}}}\bm{p}_{\mathrm{b}}\cdot\delta\bm{\chi}ds_{\mathrm{b}}=\int\limits_{{\mathcal{S}}_{\mathrm{m}}}\Big{[}-{\widehat{a}}^{1/2}\widetilde{p}\bm{n}\cdot\delta\bm{r}+2TH\overline{p}\bm{n}\cdot\delta\bm{r}-\widehat{\text{Div}}\left(T\overline{p}\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)\cdot\delta\bm{r}\Big{]}d{\mathcal{S}}_{\mathrm{m}}$
$\displaystyle+\int\limits_{{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}_{\mathrm{m}}^{\mathrm{u}}}T\overline{p}\overline{\bm{\nu}}\cdot\delta\bm{r}dl.$
(4.32)

where the difference in the top and bottom pressure is defined by
$\widetilde{p}=p_{\mathrm{t}}-p_{\mathrm{b}}$.

The contribution due to the dead load traction on the Neumann boundary of the
shell’s mid-surface is written as

$\displaystyle\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}\bm{t}_{\ell}\cdot\delta\bm{\chi}dl=\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}\bm{t}_{\ell}\cdot\delta\bm{\chi}\Big{|}_{\eta=0}dl=\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}\bm{t}_{\ell}\cdot\delta\bm{r}dl.$ (4.33)

The virtual work due to the body force per unit volume can be further
simplified as

$\displaystyle\int\limits_{{\mathcal{B}}_{0}}\bm{\mathfrak{B}}\cdot\delta\bm{\chi}dV$
$\displaystyle=$
$\displaystyle\int\limits_{S_{\mathrm{m}}}\left[T\bm{\mathfrak{B}}_{0}\cdot\delta\bm{r}+\mathcal{O}\left(T^{3}\right)\right]dS_{\mathrm{m}},$
(4.34)

with
$\bm{\mathfrak{B}}(\eta,\theta^{\alpha})=\bm{\mathfrak{B}}_{0}(\theta^{\alpha})+\eta\bm{\mathfrak{B}}_{1}(\theta^{\alpha})+\bm{\mathcal{O}}(\eta^{2})$
in ${\mathcal{B}}_{0}$.

Total contribution:

Considering Equations (4.32), (4.33), and (4.34), and neglecting the higher
order terms, the total contribution due to the external loads can be written
in terms of the variation $\delta\bm{r}$ as

$\displaystyle-$
$\displaystyle\int\limits_{s_{\mathrm{t}}}\bm{p}_{\mathrm{t}}\cdot\delta\bm{\chi}ds_{\mathrm{t}}-\int\limits_{s_{\mathrm{b}}}\bm{p}_{\mathrm{b}}\cdot\delta\bm{\chi}ds_{\mathrm{b}}-\int\limits_{{\mathcal{B}}_{0}}\bm{\mathfrak{B}}\cdot\delta\bm{\chi}dV-\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}\bm{t}_{\ell}\cdot\delta\bm{\chi}dl$
$\displaystyle=$
$\displaystyle\int\limits_{{\mathcal{S}}_{\mathrm{m}}}\left[\big{[}{\widehat{a}}^{1/2}\widetilde{p}-2TH\overline{p}\big{]}\bm{n}+\widehat{\text{Div}}\left(T\overline{p}\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)+T\bm{\mathfrak{B}}_{0}\right]\cdot\delta\bm{r}d{\mathcal{S}}_{\mathrm{m}}-\int\limits_{C_{\mathrm{m}}\setminus
C_{\mathrm{m}}^{\mathrm{u}}}\left[T\overline{p}\overline{\bm{\nu}}+\bm{t}_{\ell}\right]\cdot\delta\bm{r}dl.$
(4.35)

####  4.2.3 Shell equations in weak form

Starting with the first variation of the total energy functional in Equation
(4.3), the total contribution due to the Piola stress terms (internal virtual
work) were derived in Section 4.2.1 and the total contribution due to the
pressure, traction, and body force in Section 4.2.2. Combining the total
contributions in Equations (4.25) and (4.35), the first variation of energy
can be written in terms of the shell variables as

$\displaystyle\delta\widetilde{\Pi}[\bm{r},p_{0};\delta\bm{r},\delta p_{0}]=$
$\displaystyle\int\limits_{S_{\mathrm{m}}}\Bigg{[}-T\,\widehat{\text{Div}}\bm{P}_{0}+\widehat{\text{Div}}\left(T\widetilde{P}\bm{a}^{\alpha}\otimes\bm{A}_{\alpha}\right)-2TH\bm{T}_{0}$
$\displaystyle+\big{[}{\widehat{a}}^{1/2}\widetilde{p}-2TH\overline{p}\big{]}\bm{n}+T\bm{\mathfrak{B}}_{0}\Bigg{]}\cdot\delta\bm{r}d{\mathcal{S}}_{\mathrm{m}}-\int\limits_{S_{\mathrm{m}}}T\left[J_{0}-1\right]\delta
p_{0}dS_{\mathrm{m}}$
$\displaystyle+\int\limits_{{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}_{\mathrm{m}}^{\mathrm{u}}}\bigg{[}T\Big{[}\bm{P}_{0}{\bm{\nu}}-\widetilde{P}\overline{\bm{\nu}}\Big{]}-\bm{t}_{\ell}\bigg{]}\cdot\delta\bm{r}dl.$
(4.36)

Here $\widetilde{P}=\overline{p}+\lambda[\bm{T}_{0}\cdot\bm{n}]$ can be
physically interpreted as the “squeezing” pressure on the shell wall. The
necessary governing equations and boundary conditions are obtained by setting
$\delta\widetilde{\Pi}=0$.

The formulation is expressed in terms of two unknowns, the displacement vector
$\bm{r}$ of the shell mid-surface and the zeroth order term of the Lagrange
multiplier $p_{0}$ on the shell mid-surface. Since the top and bottom surfaces
were treated separately to make the model “embedding-aware” by accounting for
the appropriate boundary conditions, several non-standard terms arise in
comparison to the classical thin shell models (Cirak and Ortiz, 2001; Kiendl
et al., 2015; Liu et al., 2023). In particular, the pressure on the top and
bottom surface can be distinguished due to the terms $\widetilde{p}$ and
$\overline{p}$ as opposed to a single variable for pressure at the mid-surface
as appears in classical models. The effect of Piola stress is directly
accounted for by the terms involving $\bm{P}_{0}$ and its derivatives as
opposed to taking a stress resultant by integrating $\bm{P}$ through the
thickness in classical models. The influence of these terms in this novel
model are discussed with the aid of examples in Section 5.

###  4.3 Shell governing equations and boundary conditions

The necessary partial differential equations for the shell model are obtained
as the Euler–Lagrange equations by vanishing of the first variation
$\delta\widetilde{\Pi}=0$ in Equation (4.36). Since the variations
$\delta\bm{r}$ and $\delta p_{0}$ are arbitrary, the following PDE is obtained
over the mid-surface of the shell:

$\displaystyle T\,\widehat{\text{Div}}\bm{P}_{0}$
$\displaystyle-\widehat{\text{Div}}\left(T\Big{[}\lambda\left[\bm{T}_{0}\cdot\bm{n}\right]\bm{a}^{\alpha}+\overline{p}\bm{a}^{\alpha}\Big{]}\otimes\bm{A}_{\alpha}\right)$
$\displaystyle+2TH\bm{T}_{0}-\big{[}{\widehat{a}}^{1/2}\widetilde{p}-2TH\overline{p}\big{]}\bm{n}+T\bm{\mathfrak{B}}_{0}=\bm{0}\quad\forall\bm{X}\in
S_{\mathrm{m}},$ (4.37)

along with the incompressibility condition

$J_{0}-1=0\quad\forall\bm{X}\in S_{\mathrm{m}},$ (4.38)

and the boundary condition on the Neumann boundary of the mid-surface

$T\Big{[}\bm{P}_{0}{\bm{\nu}}-\widetilde{P}\overline{\bm{\nu}}\Big{]}-\bm{t}_{\ell}=\bm{0}\quad\forall\bm{X}\in{\mathcal{C}}_{\mathrm{m}}\setminus{\mathcal{C}}_{\mathrm{m}}^{\mathrm{u}}.$
(4.39)

The governing Equation (4.37) can be simplified by assuming uniform thickness
in the reference configuration and employing the plane-stress approximation
$\bm{T}_{0}\cdot\bm{n}=0$. Application of the plane-stress condition yields an
explicit expression for the Lagrange multiplier $p_{0}$ in terms of the
deformation gradient and the chosen hyperelastic energy density function. By
recognising that the Piola stress can be decomposed as
$\bm{P}_{0}={{P}_{0}}^{\alpha\beta}\bm{a}_{\alpha}\otimes\bm{A}_{\beta}+{P_{0}}_{\perp}\bm{n}\otimes\bm{N}$
where ${P_{0}}_{\perp}=0$ due to the plane-stress approximation, the equations
for in-plane and normal directions can be separated and when written in index
notation take the form

$\displaystyle{{P}_{0}}^{\gamma\delta}_{,\delta}+{{P}_{0}}^{\gamma\delta}{\Gamma}^{\alpha}_{\delta\alpha}-{{P}_{0}}^{\zeta\delta}\bm{a}_{\zeta}\cdot\bm{a}^{\gamma}_{,\delta}-\bigg{[}\overline{p}_{,\alpha}+\Big{[}{\Gamma}^{\beta}_{\beta\alpha}-{\gamma}^{\beta}_{\beta\alpha}\Big{]}\overline{p}\bigg{]}a^{\alpha\gamma}$
$\displaystyle=0\quad\forall\bm{X}\in S_{\mathrm{m}},$ (4.40)
$\displaystyle{{P}_{0}}^{\gamma\delta}a_{\gamma\beta}b^{\
\beta}_{\delta}+\widehat{a}^{1/2}T^{-1}\widetilde{p}+2[h-H]\overline{p}$
$\displaystyle=0\quad\forall\bm{X}\in S_{\mathrm{m}}.$ (4.41)

Here we note the departure from classical thin shell formulations where the
PDEs are presented in a homogenised form by integrating over the thickness and
defining average variables (Kiendl et al., 2015). There equations are
expressed in terms of bending moment and normal forces, with their conjugates
representing the change in curvature of the deforming body and the virtual
strain.

The equations derived here are used in the subsequent section to study the
response of nonlinear thin shells. The systems under consideration reflect the
cases of a thin plate, without any reference curvature of the mid-surface, and
an infinite cylinder with one nonzero component of the curvature tensor when
expressed in the local triads.

##  5 Illustrative examples

The primary aim of this section is to illustrate the theory of hyperelastic
thin shells derived in Section 4. The Euler–Lagrange equations presented in
Section 4.3 are used to derive the response equations and solve the boundary
value problems for the deformation of a flat plate and a cylindrical shell.
For analytical derivations, a simple form of an incompressible neo-Hookean
strain energy denstity function is used:

$W=\mu[I_{1}-3],$ (5.1)

where $\mu$ is the shear modulus. For this model, the Piola stress is
calculated as

$\bm{P}=\mu\bm{F}-p\bm{F}^{-\mathrm{T}},$ (5.2)

and its zeroth order component from Equation (3.16) is written in terms of the
local triad of the shell mid-surface as

$\displaystyle\bm{P}_{0}=\left[\mu
A^{\alpha\beta}-p_{0}a^{\alpha\beta}\right]\bm{a}_{\alpha}\otimes\bm{A}_{\beta}+\left[\mu\lambda-\frac{p_{0}}{\lambda}\right]\bm{n}\otimes\bm{N}.$
(5.3)

The plane stress condition $\bm{T}_{0}\cdot\bm{n}=0$ leads to the explicit
determination of the Lagrange multiplier as

$\displaystyle p_{0}=\mu\lambda^{2}.$ (5.4)

###  5.1 Thin plate under boundary traction and pressure across surfaces

Consider the problem of the deformation of a thin hyperelastic plate with
sides of length $L$ and a uniform thickness $T$ in its reference
configuration. Two types of loading scenarios are analysed as shown in Figure
2(a,b). In the first case, a uniformly distributed traction per unit area
$t_{\ell}$ is applied at the boundary $X_{1}=L$ along the $X_{1}$ axis; while
the points on the boundary $X_{1}=0$ are free to move along the $X_{2}$ axis
but their motion along the $X_{1}$ direction is constrained. In the second
case, the plate is squeezed by a uniform pressure $p$ on its two surfaces at
$X_{3}=\pm T/2$.

The homogeneous deformation of the plate can be represented by the
3-dimensional deformation gradient
$[\bm{F}]=\text{diag}(\lambda_{1},\lambda_{2},\lambda)$. Considering the mid-
surface of the shell, the covariant metric tensors in the reference and the
deformed configurations are given by

$\displaystyle\left[{A}_{\alpha\beta}\right]=\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\quad\text{and}\quad\left[{a}_{\alpha\beta}\right]=\begin{bmatrix}\lambda^{2}_{1}&0\\\
0&\lambda^{2}_{2}\end{bmatrix},$ (5.5)

with the determinants $A=1$ and $a=\lambda_{1}^{2}\lambda_{2}^{2}$. The
curvature tensor vanishes for the flat plate. The incompressibility condition
(4.38) gives the thickness stretch as $\lambda=\sqrt{A/a}$ and imposes the
constraint $\lambda\lambda_{1}\lambda_{2}=1$.

$L$
---
$\mathrm{midsurface}$
---
$t_{\ell}$
---
$x_{2},\ X_{2}$
---
$x_{3},\ X_{3}$
---
$x_{1},\ X_{1}$
---

$x_{3},\ X_{3}$
---
$x_{1},\ X_{1}$
---
$x_{2},\ X_{2}$
---
$p$
---
$p$
---

0
---
0.1
---
0.2
---
0.3
---
0.4
---
0
---
0.5
---
1
---
1.5
---
2
---
2.5
---
3
---

$\mathrm{Stretch}$
---

$t_{\ell}/\mu$
---
$\grave{T}=1/20,\ \lambda_{1}$
---
$\grave{T}=1/20,\ \lambda_{2}$
---
$\grave{T}=1/30,\ \lambda_{1}$
---
$\grave{T}=1/30,\ \lambda_{2}$
---

0
---
2
---
4
---
6
---
8
---
1
---
1.4
---
1.8
---
2.2
---
2.6
---
3
---
$p/\mu$
---
$\lambda_{1}$
---

Figure 2:  (a) Schematic of a thin plate with roller support at $X_{1}=0$ and
a uniform traction applied at $X_{1}=L$. (b) Schematic of a thin plate with
equal and opposite pressure $p$ applied on the two surfaces at $X_{3}=\pm
T/2$. (c) Response of the plate under traction governed by Equation (5.7) for
two thickness values ($\grave{T}=\\{1/20,1/30\\}$). (d) Response of the plate
squeezed by equal and opposite pressure governed by Equation (5.9).

For the first problem with traction, Equations (4.40) and (4.41) are trivially
satisfied due to homogeneous deformation. The boundary condition (4.39)
applied on the traction-free boundary $X_{2}=\\{0,L\\},\bm{\nu}=\pm\bm{e}_{2}$
along with Equation (5.4) results in the relation

$\lambda_{2}=\lambda.$ (5.6)

Application of Equation (4.39) on the boundary $X_{1}=L,\bm{\nu}=\bm{e}_{1}$
results in the response equation

$\displaystyle\frac{t_{\ell}}{\mu
L}=\grave{T}\left[\lambda_{1}-\lambda_{1}^{-2}\right]=\grave{T}\left[\lambda_{2}^{-2}-\lambda_{2}^{4}\right],$
(5.7)

where the dimensionless thickness $\grave{T}=\dfrac{T}{L}$ is defined in order
to incorporate the effect of thickness on the behaviour of the system. This
response is plotted in Fig. 2(c). It is observed that as traction $t_{\ell}$
is applied on one boundary while the other boundary remains free, the thin
plate deforms exhibiting a monotonic increase in $\lambda_{1}$ and a monotonic
decrease in $\lambda_{2}$. The effect of the thickness of the thin plate is
also depicted, and it is observed that for a thinner case, the system is
easier to deform.

For the second problem where an equal and opposite pressure $p$ is applied on
the two surfaces of the plate in the deformed configuration, symmetry and
incompressibility lead to the condition

$\lambda_{1}=\lambda_{2}=\lambda^{-2}.$ (5.8)

The equilibrium Equations (4.40) and (4.41) are trivially satisfied due to
homogeneity of the deformation. Substituting $\bm{t}_{\ell}=\bm{0}$ and
$\overline{p}=p$ in the boundary condition (4.39) results in

$\frac{p}{\mu}=\lambda_{1}^{2}-\lambda_{1}^{-4}=\lambda_{2}^{2}-\lambda_{2}^{-4}.$
(5.9)

This response is plotted in Figure 2(d), and it is evident that as the applied
pressure increases, the stretch $\lambda_{1}$ increases monotonically. It is
noteworthy that the shell model is able to capture this effect even with the
plane-stress assumption in place. Several classical thin shell models (Cirak
and Ortiz, 2001; Kiendl et al., 2015; Liu et al., 2023) are unable to capture
this loading condition because they consider an average pressure acting
directly on the mid-surface. It is noted that the exact same response
Equations (5.7) and (5.9) are obtained from the equations of 3D elasticity
(Ogden, 1997) if the plate is considered as a 3D block, thereby validating the
current model.

###  5.2 Inflation of a thin cylindrical shell

Consider an infinitely long thin hyperelastic cylindrical shell as shown in
Figure 3(a). Different magnitudes of pressure are applied at the internal and
the external surfaces of the shell while a unit stretch is maintained along
the $z$ axis resulting in an axisymmetric inflation.

0
---
0.01
---
0.02
---
0.03
---
1
---
2
---
3
---
4
---
5
---
6
---
7
---
8
---
9
---
10
---
$\lambda_{\theta}$
---
$\Delta p$
---
$\tilde{T}$
---
$p_{t}/\mu$
---
$\tilde{T}=1/20,p_{t}/\mu=0$
---
$\tilde{T}=1/30,p_{t}/\mu=0$
---
$p_{t}/\mu=0.08$
---
$\mathrm{3D}\ \mathrm{case}$
---
$\mathrm{3D}\ \mathrm{case}$
---
$p_{t}/\mu=0.08$
---

Figure 3:  (a) A schematic showing the deformed configuration of an infinite
cylindrical shell with thickness $t$ under the influence of pressures, $p_{b}$
and $p_{t}$, at the inner and outer surfaces respectively. (b) Variation of
the mid-surface azimuthal stretch ($\lambda_{\theta}$) with the dimensionless
pressure difference $\Delta p=[p_{b}-p_{t}]/\mu$ for the cylinder modelled as
a thin shell and a 3D solid. The shell response can be controlled by changing
the pressure on the outer surface $p_{t}/\mu$ and the response for the 3D
model lies within the region $p_{t}/\mu\in(0,0.08)$.

The axisymmetric deformation of a cylindrical shell can be written as

$\theta=\Theta,\quad\quad
z=Z,\quad\quad\bm{r}=\bm{R}+\bm{u}=\left[R+u\right]\bm{e}_{\rho}+Z\bm{e}_{z}=r\bm{e}_{\rho}+Z\bm{e}_{z}$
(5.10)

Here, $\theta$ and $z$ denote the deformed coordinates corresponding to their
azimuthal and axial counterparts in the reference configuration ($\Theta$ and
$Z$, respectively). The unit vectors along the axial and radial directions are
represented by $\bm{e}_{z}$ and $\bm{e}_{\rho}$, respectively. Moreover, $R$
and $r$ identify the radius at the mid-surface of the cylindrical shell in the
two configurations. The displacement vector is expressed as
$\bm{u}=u(\rho)\bm{e}_{\rho}$. The covariant and contravariant vectors at the
mid-surface in the two configurations, along with the reference and deformed
normals, are specified by

$\displaystyle{\color[rgb]{1,1,1}=}$
$\displaystyle{\bm{A}}_{1}=R\bm{e}_{\theta},\quad{\bm{A}}_{2}=\bm{e}_{z},\quad{\bm{A}}^{1}=\frac{1}{R}\bm{e}_{\theta},\quad{\bm{A}}^{2}=\bm{e}_{z},$
(5.11) $\displaystyle{\color[rgb]{1,1,1}=}$
$\displaystyle{\bm{a}}_{1}=r\bm{e}_{\theta},\quad{\bm{a}}_{2}=\bm{e}_{z},\quad{\bm{a}}^{1}=\frac{1}{r}\bm{e}_{\theta},\quad{\bm{a}}^{2}=\bm{e}_{z},$
$\displaystyle{\color[rgb]{1,1,1}=}$
$\displaystyle\bm{n}=\bm{N}=\bm{e}_{\rho}.$

The azimuthal unit vector is denoted by $\bm{e}_{\theta}$. The components of
the covariant metric tensors at the mid-surface in the reference and deformed
configurations are:

$\displaystyle\left[{A}_{\alpha\beta}\right]=\begin{bmatrix}{R}^{2}&0\\\
0&1\end{bmatrix}\quad\text{and}\quad\left[{a}_{\alpha\beta}\right]=\begin{bmatrix}{r}^{2}&0\\\
0&1\end{bmatrix},$ (5.12)

along with the determinant of the covariant metric tensors at the mid-surface
given by

$A={R}^{2},\quad\text{and}\quad a={r}^{2}.$ (5.13)

The incompressibility condition (4.38) results in

$J_{0}-1=0\Rightarrow\lambda=\sqrt{\frac{A}{a}}=\lambda_{\theta}^{-1},$ (5.14)

with $\lambda_{\theta}=r/R$ as the azimuthal stretch. The non-zero components
of the curvature tensor at the mid-surface are

${B}_{1}^{\ 1}=-\frac{1}{R},\quad\text{and}\quad{b}_{1}^{\ 1}=-\frac{1}{r}.$
(5.15)

In this case, the in-plane governing Equation (4.40) is trivially satisfied
due to axisymmetry of the deformation. Upon substituting
$\overline{p}=[p_{\mathrm{t}}+p_{\mathrm{b}}]/2$,
$\widetilde{p}=p_{\mathrm{t}}-p_{\mathrm{b}}$ and the expression for stress
from Equation (5.3) in the out-of-plane balance Equation (4.41), one obtains

$\displaystyle\dfrac{p_{\mathrm{b}}-p_{\mathrm{t}}}{\mu}=\dfrac{\widetilde{T}\left[\lambda-1\right]\dfrac{p_{t}}{\mu}+\widetilde{T}\left[1-\lambda^{4}\right]}{\dfrac{\widetilde{T}}{2}\left[\lambda-1\right]+\lambda^{-1}}$
(5.16)

where $\widetilde{T}=T/R$ is the dimensionless thickness of the cylindrical
shell. Pressures at the inner and outer surfaces of the cylinder are denoted
by $p_{b}$ and $p_{t}$, respectively. When the response of the system is
derived from the perspective of considering the cylinder as a three
dimensional hyperelastic body, the response expression is given by (Haughton
and Ogden, 1979)

$\displaystyle\dfrac{p_{\mathrm{b}}-p_{\mathrm{t}}}{\mu}=6\left[\lambda_{\mathrm{b}}^{-4}-\lambda_{\mathrm{a}}^{-4}\right]-2\left[\ln{\lambda_{\mathrm{b}}}-\ln{\lambda_{\mathrm{a}}}\right],$
(5.17)

where $\lambda_{\mathrm{a}}$ and $\lambda_{\mathrm{b}}$ are the stretches at
the inner and outer surfaces of the infinite cylinder and are related by the
incompressibility constraint

$\displaystyle\lambda_{\mathrm{b}}^{2}-1={\left[\frac{R_{\mathrm{out}}-R_{\mathrm{in}}}{R_{\mathrm{in}}}\right]}^{2}\left[\lambda_{\mathrm{a}}^{2}-1\right],$
(5.18)

where $R_{\mathrm{in}}$ and $R_{\mathrm{out}}$ ane the radius at the inner and
outer surfaces of the cylindrical shell. Moreover, let $\lambda_{\mathrm{c}}$
represent the stretch at the mid-surface of the cylinder, which can be
similarly connected to the stretches at the inner and outer surfaces from the
incompressibility relation to eventually obtain a response function for the
difference in pressures corresponding to $\lambda_{\mathrm{c}}$. These 3D
relations allow for direct comparison with the shell system of equations.

The comparison between the 3D and shell theory is presented in Figure 3(b) for
various shell thickness values and external pressure values ($p_{t}/\mu$). As
the pressure difference ($\Delta p=[p_{b}-p_{t}]/\mu$) between the inner and
external shell surface increases, the stretch $\lambda_{\theta}$ monotonically
increases until reaching a critical value $\Delta p_{\mathrm{c}}$
corresponding to a limit point instability. At this juncture, minor changes in
applied pressure lead to significant inflationary shifts. Similar limit point
instabilities have been observed during inflation of thin hyperelastic shells
and soft cylindrical cavities (Kiendl et al., 2015; Cheewaruangroj et al.,
2019; Mehta et al., 2022). The critical limit point pressure decreases with
decreasing shell thickness. This example demonstrates the distinction between
considering the pressure on the top and bottom surfaces of the shell
separately, as opposed to the common practice of considering a pressure
difference on the mid-surface. Varying the pressure on the external surface
$p_{t}/\mu$ significantly alters the shell’s response to the applied pressure
difference $\Delta p$. Reduction in shell thickness brings these response
curves closer together, as observed for $\widetilde{T}=1/30$ and $1/20$.

##  6 Concluding Remarks

This study rigorously derived the governing equations for the embedding aware
large deformation of Kirchhoff-Love hyperelastic thin shells, using a
variational approach and dimension reduction. The general deformation map
within the shell was redefined to introduce terms for through-thickness
stretch and the deformed normal. Deformation function for the shell’s mid-
surface and the Lagrange multiplier associated with the incompressibility
constraint were employed as the generalised set of solution variables for a
two-field variational formulation, leading to the derivation of equilibrium
equations for the shell.

The formulation preserves boundary conditions at the top and bottom surfaces
of the shell thereby accounting for the fields in the embedding space. In the
current work, this manifests as a distinction between hydrostatic pressures
applied to the top and bottom surfaces. This lays a foundation for multi-
physics shell models of electro- and magnetoelasticity that must account for
perturbations in the electromagnetic fields in the embedding media due to
shell deformation. This enriched approach effectively models large
deformations and limit point instabilities, as demonstrated by scenarios such
as the stretching and squashing of a soft thin plate and the inflation of an
infinitely long cylindrical hyperelastic shell.

This work provides a fresh perspective on shell system of equations,
emphasising geometric precision and accuracy by capturing arbitrary geometric
and constitutive nonlinearities. The derived equations, rooted in a
variational formulation, are well-prepared for future numerical implementation
via the finite element method.

## Acknowledgements

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) grants EP/V030833/1 and EP/R008531/1, and a Royal Society
grant IES/R1/201122.

##  Appendix A Geometry of a Kirchoff-Love thin shell

###  A.1 The natural basis at the mid-surface

The covariant basis vectors for the mid-surface in the reference and deformed
configurations, respectively, can be computed as

$\bm{A}_{\alpha}=\frac{\partial{\bm{x}}}{\partial\theta^{\alpha}},\quad\mathrm{and}\quad\bm{a}_{\alpha}=\frac{\partial\bm{x}}{\partial\theta^{\alpha}}.$
(A.1)

Thus, the unit normal vectors in the two configurations are defined by

$\bm{N}=\frac{{\bm{A}}_{1}\times{\bm{A}}_{2}}{{A}^{1/2}},\quad\text{and}\quad\bm{n}=\frac{\bm{a}_{1}\times\bm{a}_{2}}{{a^{1/2}}},$
(A.2)

where $A$ and $a$ are

$A=\mathinner{\\!\left\lVert{\bm{A}}_{1}\times{\bm{A}}_{2}\right\rVert}^{2},\quad\text{and}\quad
a=\mathinner{\\!\left\lVert\bm{a}_{1}\times\bm{a}_{2}\right\rVert}^{2}.$ (A.3)

Further, it can be shown that

$A=\mathrm{det}\left[{A}_{\alpha\beta}\right],\quad\text{and}\quad
a=\mathrm{det}[{a}_{\alpha\beta}].$ (A.4)

The covariant components of the metric tensor for the mid-surface points
$\bm{R}$ and $\bm{r}$ are respectively given by

${A}_{\alpha\beta}={\bm{A}}_{\alpha}\cdot{\bm{A}}_{\beta},\quad\text{and}\quad{a}_{\alpha\beta}=\bm{{a}}_{\alpha}\cdot\bm{{a}}_{\beta}.$
(A.5)

Also, the contravariant metric tensor components for the mid-surface are

${A}^{\alpha\gamma}{A}_{\gamma\beta}=\delta^{\alpha}_{\beta},\quad\text{and}\quad{a}^{\alpha\gamma}{a}_{\gamma\beta}=\delta^{\alpha}_{\beta},$
(A.6)

where $\delta^{\alpha}_{\beta}$ denotes the Kronecker Delta. Again,
${\bm{A}}^{\alpha}$ and $\bm{a}^{\alpha}$ denotes the contravariant basis
vectors for the mid-surface in the two configurations, and they are defined as

${\bm{A}}^{\alpha}\cdot{\bm{A}}_{\beta}=\delta^{\alpha}_{\beta},\quad\text{and}\quad\bm{a}^{\alpha}\cdot\bm{a}_{\beta}=\delta^{\alpha}_{\beta}.$
(A.7)

###  A.2 The unit alternator and permutation symbol

In general, for a surface tensor ${\bm{Q}}$, the surface inverse
${\bm{Q}}^{-1}$ defined from

${\bm{Q}}^{-1}{\bm{Q}}={\bm{I}},$ (A.8)

with
$\bm{I}={\bm{A}}^{\beta}\otimes{\bm{A}}_{\beta}={\bm{A}}_{\beta}\otimes{\bm{A}}^{\beta}$
($\bm{I}$ denotes the projection onto the tangent plane of ${S}_{\mathrm{m}}$)
has the contravariant components as

${Q}_{\mathrm{inv}}^{\alpha\beta}=\frac{1}{Q}{e}^{\alpha\gamma}{Q}_{\delta\gamma}{e}^{\beta\delta},$
(A.9)

where $Q=\mathrm{det}\left[{Q}_{\alpha\beta}\right]$, and the so-called unit
alternator given as

$\left[e^{\alpha\gamma}\right]=\begin{bmatrix}0&1\\\ -1&0\end{bmatrix}.$
(A.10)

Further, the permutation tensor is defined as

$\bm{E}={E}^{\alpha\beta}\
{\bm{A}}_{\alpha}\otimes{\bm{A}}_{\beta}=\frac{1}{{A}^{1/2}}{e}^{\alpha\beta}\
{\bm{A}}_{\alpha}\otimes{\bm{A}}_{\beta},\quad\text{and}\quad\bm{\varepsilon}={\varepsilon}^{\alpha\beta}\
\bm{a}_{\alpha}\otimes\bm{a}_{\beta}=\frac{1}{{a}^{1/2}}e^{\alpha\beta}\
\bm{a}_{\alpha}\otimes\bm{a}_{\beta},$ (A.11)

in the two configurations. In particular, Equation (A.9) yields

${A}^{\alpha\beta}=\frac{1}{{A}}{e}^{\alpha\gamma}{A}_{\gamma\delta}{e}^{\beta\delta},$
(A.12)

and similarly,

$a^{\alpha\beta}=\frac{1}{a}e^{\alpha\gamma}a_{\gamma\delta}e^{\beta\delta}.$
(A.13)

From the above, and using the relation,
$e^{\gamma\alpha}e_{\gamma\beta}=\delta^{\alpha}_{\beta}$, we get

${A}^{\alpha\beta}{e}_{\beta\gamma}{A}={e}^{\alpha\beta}{A}_{\beta\gamma}\quad\text{and}\quad{a}^{\alpha\beta}{e}_{\beta\gamma}a={e}^{\alpha\beta}{a}_{\beta\gamma},$
(A.14)

which can be further used to rewrite the permutation tensors as

$\bm{E}={E}_{\alpha\beta}\
{\bm{A}}^{\alpha}\otimes{\bm{A}}^{\beta}={A}^{1/2}{e}_{\alpha\beta}\
{\bm{A}}^{\alpha}\otimes{\bm{A}}^{\beta},\quad\text{and}\quad\bm{\varepsilon}={\varepsilon}_{\alpha\beta}\
\bm{a}^{\alpha}\otimes\bm{a}^{\beta}={a}^{1/2}e_{\alpha\beta}\
\bm{a}^{\alpha}\otimes\bm{a}^{\beta}.$ (A.15)

Again, multiplying Equations (A.12) and (A.13) by ${A}_{\alpha\beta}$ and
$a_{\alpha\beta}$, respectively, we obtain

$A=\frac{1}{2}{e}^{\alpha\gamma}{e}^{\beta\delta}{A}_{\alpha\beta}{A}_{\gamma\delta},\quad\text{and}\quad
a=\frac{1}{2}e^{\alpha\gamma}e^{\beta\delta}a_{\alpha\beta}a_{\gamma\delta}.$
(A.16)

From the above, we can write

${A}_{,\zeta}={e}^{\alpha\gamma}{e}^{\beta\delta}{A}_{\alpha\beta,\zeta}{A}_{\gamma\delta},\quad\text{and}\quad
a_{,\zeta}=e^{\alpha\gamma}e^{\beta\delta}a_{\alpha\beta,\zeta}a_{\gamma\delta},$
(A.17)

which can be further rewritten as

${A}_{,\zeta}={A}{A}^{\alpha\beta}{A}_{\alpha\beta,\zeta},\quad\text{and}\quad
a_{,\zeta}=aa^{\alpha\beta}a_{\alpha\beta,\zeta},$ (A.18)

by using Equations (A.12) and (A.13). Now,

${A}_{\alpha\beta,\zeta}={\bm{A}}_{\alpha,\zeta}\cdot{\bm{A}}_{\beta}+{\bm{A}}_{\alpha}\cdot{\bm{A}}_{\beta,\zeta},\quad\text{and}\quad{a}_{\alpha\beta,\zeta}={\bm{a}}_{\alpha,\zeta}\cdot{\bm{a}}_{\beta}+{\bm{a}}_{\alpha}\cdot{\bm{a}}_{\beta,\zeta}.$
(A.19)

Therefore,

${A}_{,\zeta}=2{A}{\Gamma}^{\alpha}_{\alpha\zeta},\quad\text{and}\quad{a}_{,\zeta}=2{a}{\gamma}^{\alpha}_{\alpha\zeta},$
(A.20)

with the Christoffel symbols of the second kind in the two configurations
defined by

${\Gamma}^{\alpha}_{\zeta\gamma}={\bm{A}}^{\alpha}\cdot{\bm{A}}_{\zeta,\gamma},\quad\text{and}\quad{\gamma}^{\alpha}_{\zeta\gamma}={\bm{a}}^{\alpha}\cdot{\bm{a}}_{\zeta,\gamma}.$
(A.21)

###  A.3 The natural basis at a shell-point

A point $\bm{x}\in\mathcal{B}$ can be written as

$\bm{x}=\bm{r}+\eta\bm{d},$ (A.22)

where $\bm{d}=\lambda\bm{n}$ and $\lambda=\dfrac{t}{T}$. The covariant basis
vectors at a point $\bm{x}$ in the shell are

$\displaystyle\bm{g}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{\partial\bm{x}}{\partial\theta^{\alpha}},$ (A.23)
$\displaystyle=$
$\displaystyle\frac{\partial\bm{r}}{\partial\theta^{\alpha}}+\eta\bm{n}\frac{\partial\lambda}{\partial\theta^{\alpha}}+\eta\lambda\frac{\partial\bm{n}}{\partial\bm{r}}\frac{\partial\bm{r}}{\partial\theta^{\alpha}},$
$\displaystyle=$ $\displaystyle\bm{\mu}\bm{a}_{\alpha}.$

In the above, neglecting the variation of the through-thickness stretch along
the mid-surface of the shell-structure (Kiendl et al., 2015; Liu et al.,
2023), we arrive at

$\bm{\mu}=\bm{i}-\eta\lambda\bm{\kappa},$ (A.24)

where
$\bm{i}=\bm{a}^{\beta}\otimes\bm{a}_{\beta}=\bm{a}_{\beta}\otimes\bm{a}^{\beta}$
denotes the projection onto the tangent plane of $s_{\mathrm{m}}$, the
deformed counterpart of ${S}_{\mathrm{m}}$. Also,

$\bm{\kappa}=-\dfrac{\partial\bm{n}}{\partial\bm{r}}=-{\bm{n}}_{,\beta}\otimes{\bm{a}}^{\beta}.$
(A.25)

The covariant basis vectors at a point $\bm{X}$ in the shell are given by

$\displaystyle{\bm{G}}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{\partial\bm{X}}{\partial\theta^{\alpha}},$ (A.26)
$\displaystyle=$
$\displaystyle\frac{\partial\bm{R}}{\partial\theta^{\alpha}}+\eta\frac{\partial\bm{N}}{\partial\bm{R}}\frac{\partial\bm{R}}{\partial\theta^{\alpha}},$
$\displaystyle=$ $\displaystyle\bm{M}{\bm{A}}_{\alpha},$

where

$\bm{M}=\bm{I}-\eta\bm{K},$ (A.27)

with

$\bm{K}=-\dfrac{\partial\bm{N}}{\partial\bm{R}}=-{\bm{N}}_{,\beta}\otimes{\bm{a}}^{\beta}.$
(A.28)

Again, ${\bm{N}}_{,\beta}$ and ${\bm{n}}_{,\beta}$ appears in the formula of
Weingarten as

${\bm{N}}_{,\beta}=-{B}_{\beta}^{\ \gamma}\
{\bm{A}}_{\gamma},\quad\text{and}\quad{\bm{n}}_{,\beta}=-{{b}}_{\beta}^{\
\gamma}\ {\bm{a}}_{\gamma},$ (A.29)

with the surface curvature tensors in the two configurations defined as

$\bm{B}={B}_{\beta\delta}\
{\bm{A}}^{\beta}\otimes{\bm{A}}^{\delta}\quad\text{and}\quad{\bm{b}}={b}_{\beta\delta}\
{\bm{a}}^{\beta}\otimes{\bm{a}}^{\delta},$ (A.30)

where

${B}_{\beta\delta}=\bm{N}\cdot{\bm{A}}_{\beta,\delta},\quad\text{and}\quad{b}_{\beta\delta}={\bm{n}}\cdot{\bm{a}}_{\beta,\delta},$
(A.31)

and further,

${{B}_{\beta}}^{\gamma}={B}_{\beta\delta}{A}^{\delta\gamma},\quad\text{and}\quad{{b}_{\beta}}^{\gamma}={b}_{\beta\delta}{a}^{\delta\gamma}.$
(A.32)

Therefore,

$\bm{K}={B}_{\beta}^{\ \gamma}\
{\bm{A}}_{\gamma}\otimes{\bm{A}}^{\beta},\quad\text{and}\quad\bm{\kappa}={{b}}_{\beta}^{\
\gamma}\ {\bm{a}}_{\gamma}\otimes{\bm{a}}^{\beta}.$ (A.33)

For a point in the shell, the components of the covariant and contravariant
metric tensors in the reference configuration are defined by

${G}_{\alpha\beta}={\bm{G}}_{\alpha}\cdot{\bm{G}}_{\beta}\quad\text{and}\quad{G}^{\alpha\beta}={\bm{G}}^{\alpha}\cdot{\bm{G}}^{\beta},$
(A.34)

with the deformed counterparts as

$g_{\alpha\beta}=\bm{g}_{\alpha}\cdot\bm{g}_{\beta}\quad\text{and}\quad
g^{\alpha\beta}=\bm{g}^{\alpha}\cdot\bm{g}^{\beta},$ (A.35)

where ${\bm{G}}^{\alpha}$ and $\bm{g}^{\alpha}$ denotes the contravariant
basis vectors in the two configurations defined by

${\bm{G}}^{\alpha}\cdot{\bm{G}}_{\beta}=\delta^{\alpha}_{\beta},\quad\text{and}\quad\bm{g}^{\alpha}\cdot\bm{g}_{\beta}=\delta^{\alpha}_{\beta}.$
(A.36)

It can be shown that

${\bm{G}}^{\alpha}={\bm{M}}^{-\mathrm{T}}{\bm{A}}^{\alpha},$ (A.37)

and ${\bm{M}}^{-1}$ can be expanded as

$\displaystyle{\bm{M}}^{-1}$ $\displaystyle=$
$\displaystyle{M}^{-1\gamma}_{0_{\beta}}\
{\bm{A}}_{\gamma}\otimes{\bm{A}}^{\beta}+\eta{M}^{-1\gamma}_{1_{\beta}}\
{\bm{A}}_{\gamma}\otimes{\bm{A}}^{\beta}+\eta^{2}{M}^{-1\gamma}_{2_{\beta}}\
{\bm{A}}_{\gamma}\otimes{\bm{A}}^{\beta}+\bm{\mathcal{O}}(\eta^{3}),$ (A.38)

with

${M}^{-1\gamma}_{0_{\beta}}=\delta^{\gamma}_{\beta},\quad{M}^{-1\gamma}_{1_{\beta}}={B}_{\beta}^{\
\gamma},\quad\text{and}\quad{M}^{-1\gamma}_{2_{\beta}}={B}_{\delta}^{\
\gamma}{B}_{\beta}^{\ \delta}.$ (A.39)

Similarly,

$\bm{g}^{\alpha}=\bm{\mu}^{-\mathrm{T}}\bm{a}^{\alpha},$ (A.40)

and $\bm{\mu}^{-1}$ can be expanded as

$\bm{\mu}^{-1}={\mu}^{-1\gamma}_{0_{\beta}}\
{\bm{a}}_{\gamma}\otimes{\bm{a}}^{\beta}+\eta{\mu}^{-1\gamma}_{1_{\beta}}\
{\bm{a}}_{\gamma}\otimes{\bm{a}}^{\beta}+\eta^{2}{\mu}^{-1\gamma}_{2_{\beta}}\
{\bm{a}}_{\gamma}\otimes{\bm{a}}^{\beta}+\bm{\mathcal{O}}(\eta^{3}),$ (A.41)

with

${\mu}^{-1\gamma}_{0_{\beta}}=\delta^{\gamma}_{\beta},\quad{\mu}^{-1\gamma}_{1_{\beta}}=\lambda{b}_{\beta}^{\
\gamma},\quad\text{and}\quad{\mu}^{-1\gamma}_{2_{\beta}}=\lambda^{2}{b}_{\delta}^{\
\gamma}\ {b}_{\beta}^{\ \delta}.$ (A.42)

####  A.3.1 The volume and surface elements

The volume element in the reference configuration can be expressed as

$\displaystyle dV$ $\displaystyle=$
$\displaystyle\left[{\bm{G}}_{1}\times{\bm{G}}_{2}\right]\cdot\bm{N}\
d\theta^{1}d\theta^{2}d\eta,$ (A.43) $\displaystyle=$ $\displaystyle\
\left[{\bm{A}}_{1}\times{\bm{A}}_{2}\right]\cdot\bm{N}Md\theta^{1}d\theta^{2}d\eta,$
$\displaystyle=$ $\displaystyle dSd\eta,$

where the undeformed elemental area $dS$ is given by

$dS=Md{S}_{\mathrm{m}},$ (A.44)

with $d{S}_{\mathrm{m}}$ as the area element on ${S}_{\mathrm{m}}$ written as

$dS_{\mathrm{m}}={A}^{1/2}dP,$ (A.45)

and the area element for the convected coordinates is
$dP=d\theta^{1}d\theta^{2}$. Also,

$\displaystyle M$ $\displaystyle=$ $\displaystyle\mathrm{det}{\bm{M}},$ (A.46)
$\displaystyle=$ $\displaystyle 1-2\eta H+\eta^{2}K,$

where $H$ and $K$ are the mean and Gaussian curvatures of the undeformed mid-
surface, respectively, and are expressed as

$H=\frac{1}{2}\mathrm{tr}{\bm{K}}=\frac{1}{2}\frac{\partial{\bm{N}}}{\partial{\bm{R}}}\mathrel{\mathop{\ordinarycolon}}{\bm{I}}=\frac{1}{2}{B}_{\alpha\beta}{A}^{\alpha\beta}=\frac{1}{2}{B}_{\alpha}^{\alpha},$
(A.47)

and

$K=\mathrm{det}{\bm{K}}=\mathrm{det}\left[{B}_{\alpha}^{\
\beta}\right]=\mathrm{det}\left[{B}_{\alpha\gamma}{A}^{\gamma\beta}\right]=\frac{B}{A},$
(A.48)

with $B=\mathrm{det}\left[{B}_{\alpha\beta}\right]$. Further, an elemental
area in the deformed configuration is given by

$ds=\mu{\hat{a}}^{1/2}dS_{\mathrm{m}},$ (A.49)

with the surface stretch $\hat{a}=\dfrac{a}{A}$, and

$\mu=1-2\eta\lambda h+\eta^{2}\lambda^{2}\kappa,$ (A.50)

where the mean and Gaussian curvatures of the deformed mid-surface are

$h=\frac{1}{2}{b}_{\alpha}^{\alpha},\quad\text{and}\quad\kappa=\frac{b}{a},$
(A.51)

with $b=\mathrm{det}[b_{\alpha\beta}]$. The boundaries,
$\partial{\mathcal{B}}_{\mathrm{0}}$ and $\partial\mathcal{B}$ can be written
as $\partial{\mathcal{B}}_{\mathrm{0}}=S_{\mathrm{t}}\cup S_{\mathrm{b}}\cup
S_{\mathrm{l}}$, and $\partial\mathcal{B}=s_{\mathrm{t}}\cup
s_{\mathrm{b}}\cup s_{\mathrm{l}}$, where, the subscripts, t, b, and l,
represents the top, bottom, and lateral surfaces in the two configurations,
and the top surface is the side of the boundary that is reached along the unit
outward normal vector. Therefore,

$dS_{\mathrm{t}}={M}\Big{|}_{\eta=T/2}dS_{\mathrm{m}},\quad\text{and}\quad
dS_{\mathrm{b}}={M}\Big{|}_{\eta=-T/2}d{S}_{\mathrm{m}}.$ (A.52)

Also,

$ds_{\mathrm{t}}=\mu\Big{|}_{\eta=T/2}{\hat{a}}^{1/2}dS_{\mathrm{m}},\quad\text{and}\quad
ds_{\mathrm{b}}=\mu\Big{|}_{\eta=-T/2}{\hat{a}}^{1/2}dS_{\mathrm{m}}.$ (A.53)

If the bounding curve $C_{\mathrm{m}}$ of the mid-surface $S_{\mathrm{m}}$ is
characterized by the arc-length parameter $l$, then the infinitesimal length
$dl$ between two points ${\bm{R}}(\theta^{1},\theta^{2})$ and
${\bm{R}}(\theta^{1}+d\theta^{1},\theta^{2}+d\theta^{2})$ is given by

$\displaystyle dl$ $\displaystyle=$
$\displaystyle\mathinner{\\!\left\lVert{\bm{R}}(\theta^{1}+d\theta^{1},\theta^{2}+d\theta^{2})-{\bm{R}}(\theta^{1},\theta^{2})\right\rVert},$
(A.54) $\displaystyle=$
$\displaystyle\mathinner{\\!\left\lVert{\bm{R}}(\theta^{1},\theta^{2})+\frac{\partial\bm{R}}{\partial\theta^{\alpha}}d\theta^{\alpha}-{\bm{R}}(\theta^{1},\theta^{2})\right\rVert},$
$\displaystyle=$
$\displaystyle\mathinner{\\!\left\lVert{\bm{A}}_{\alpha}d\theta^{\alpha}\right\rVert},$
$\displaystyle=$
$\displaystyle\sqrt{{\bm{A}}_{\alpha}d\theta^{\alpha}\cdot{\bm{A}}_{\beta}d\theta^{\beta}},$
$\displaystyle=$
$\displaystyle\sqrt{A_{\alpha\beta}d\theta^{\alpha}d\theta^{\beta}}.$

The tangent vector at a point $\bm{R}$ on $C_{\mathrm{m}}$ is defined as

$\bm{\tau}=\frac{d\bm{R}}{dl}=\frac{\partial\bm{R}}{\partial\theta^{\beta}}\frac{d\theta^{\beta}}{dl}={\bm{A}}_{\beta}\frac{d\theta^{\beta}}{dl},$
(A.55)

and using Equation (A.54), we get

$\displaystyle\bm{\tau}\cdot\bm{\tau}$ $\displaystyle=$
$\displaystyle{\bm{A}}_{\alpha}\cdot{\bm{A}}_{\beta}\frac{d\theta^{\alpha}}{dl}\frac{d\theta^{\beta}}{dl},$
(A.56) $\displaystyle=$
$\displaystyle{A}_{\alpha\beta}\frac{d\theta^{\alpha}}{dl}\frac{d\theta^{\beta}}{dl},$
$\displaystyle=$ $\displaystyle 1,$

implying that $\bm{\tau}$ is a unit tangent vector. Further, we define

$\bm{\nu}=\bm{E}\bm{\tau}={E}^{\alpha\beta}{A}_{\beta\gamma}\frac{d\theta^{\gamma}}{dl}{\bm{A}}_{\alpha}={E}_{\eta\delta}\frac{d\theta^{\delta}}{dl}{\bm{A}}^{\eta},$
(A.57)

such that,

$\displaystyle\bm{\nu}\cdot\bm{\tau}$ $\displaystyle=$
$\displaystyle{E}_{\alpha\beta}\frac{d\theta^{\alpha}}{dl}\frac{d\theta^{\beta}}{dl},$
(A.58) $\displaystyle=$
$\displaystyle\frac{1}{2}\left[{E}_{\alpha\beta}+{E}_{\beta\alpha}\right]\frac{d\theta^{\alpha}}{dl}\frac{d\theta^{\beta}}{dl},$
$\displaystyle=$ $\displaystyle 0.$

Again,

$\displaystyle\bm{\nu}\cdot\bm{\nu}$ $\displaystyle=$
$\displaystyle{E}^{\alpha\beta}{A}_{\beta\gamma}\frac{d\theta^{\gamma}}{dl}{\bm{A}}_{\alpha}\cdot{E}_{\eta\delta}\frac{d\theta^{\delta}}{dl}{\bm{A}}^{\eta},$
(A.59) $\displaystyle=$
$\displaystyle{E}^{\alpha\beta}{E}_{\eta\delta}\delta^{\eta}_{\alpha}{A}_{\beta\gamma}\frac{d\theta^{\gamma}}{dl}\frac{d\theta^{\delta}}{dl},$

and following the relation,
${E}^{\alpha\beta}{E}_{\eta\delta}\delta^{\eta}_{\alpha}=\left[\delta^{\alpha}_{\eta}\delta^{\beta}_{\delta}-\delta^{\alpha}_{\delta}\delta^{\beta}_{\eta}\right]\delta^{\eta}_{\alpha}=\delta^{\beta}_{\delta}$,
we get

$\displaystyle{\bm{\nu}}\cdot\bm{\nu}$ $\displaystyle=$
$\displaystyle\delta^{\beta}_{\delta}{A}_{\beta\gamma}\frac{d\theta^{\gamma}}{dl}\frac{d\theta^{\delta}}{dl},$
(A.60) $\displaystyle=$
$\displaystyle{A}_{\delta\gamma}\frac{d\theta^{\gamma}}{dl}\frac{d\theta^{\delta}}{dl},$
$\displaystyle=$ $\displaystyle\bm{\tau}\cdot\bm{\tau},$ $\displaystyle=$
$\displaystyle 1,$

implying that $\bm{\nu}$ is the in-plane unit normal to $\bm{\tau}$ on
$\mathcal{C}_{\mathrm{m}}$, and

$\bm{\nu}=\bm{\tau}\times\bm{N}.$ (A.61)

An elemental area, $d{S}_{\mathrm{l}}$, at a point $\bm{X}$ on the lateral
surface is given by

$\displaystyle dS_{\mathrm{l}}$ $\displaystyle=$
$\displaystyle\mathinner{\\!\left\lVert\frac{\partial\bm{X}}{\partial
l}\times\frac{\partial\bm{X}}{\partial\eta}\right\rVert}dl\,d\eta,$ (A.62)
$\displaystyle=$
$\displaystyle\mathinner{\\!\left\lVert\bm{M}\bm{\tau}\times\bm{N}\right\rVert}dl\,d\eta,$
$\displaystyle=$ $\displaystyle
c\mathinner{\\!\left\lVert{\bm{\tau}}_{\mathrm{l}}\times\bm{\bm{N}}\right\rVert}dl\,d\eta,$
$\displaystyle=$ $\displaystyle c\,dl\,d\eta,$

with

$c=\mathinner{\\!\left\lVert\bm{M}\bm{\tau}\right\rVert}={\left[1-2\eta\bm{K}\bm{\tau}\cdot\bm{\tau}+\eta^{2}\bm{K}\bm{\tau}\cdot\bm{K}\bm{\tau}\right]}^{1/2},\quad\text{and}\quad\bm{\tau}_{\mathrm{l}}=\frac{\bm{M}\bm{\tau}}{c},$
(A.63)

where ${\bm{\tau}}_{\mathrm{l}}$ is the unit tangent vector at a point on the
lateral surface on the bounding curve $\mathcal{C}$ of the surface $S$. The
in-plane unit normal to ${\bm{\tau}}_{\mathrm{l}}$ is given by

${\bm{\nu}}_{\mathrm{l}}={\bm{\tau}}_{\mathrm{l}}\times\bm{\bm{N}}.$ (A.64)

The above can be written as

$c{\bm{\nu}}_{\mathrm{l}}=\left[\bm{I}+\eta\left[\bm{K}-2H\bm{I}\right]\right]\bm{\nu},$
(A.65)

by using the relation,

$\bm{K}\bm{\tau}\times\bm{N}=\left[2H\bm{I}-\bm{K}\right]\bm{\nu}.$ (A.66)

##  Appendix A Application of Green’s Theorem at the mid-surface of the shell

For a scalar ${T}^{\alpha}$, consider the following integral:

$\int\limits_{P}\left[{A}^{1/2}{T}^{\alpha}\right]_{,\alpha}dP,$ (B.1)

which can be rewritten by applying the Green’s theorem as

$\displaystyle\int\limits_{P}\left[{A}^{1/2}{T}^{\alpha}\right]_{,\alpha}dP$
$\displaystyle=$
$\displaystyle\int\limits_{P}\left[\left[{A}^{1/2}{T}^{1}\right]_{,1}+\left[{A}^{1/2}{T}^{2}\right]_{,2}\right]dP,$
(B.2) $\displaystyle=$
$\displaystyle\int\limits_{\mathcal{C}_{\mathrm{p}}}{A}^{1/2}{e}_{\alpha\beta}{T}^{\alpha}d\theta^{\beta},$

where $\mathcal{C}_{\mathrm{p}}$ is the boundary of the parametric domain $P$,
and the above boundary integral can be simplified as

$\displaystyle\int\limits_{\mathcal{C}_{\mathrm{p}}}{A}^{1/2}{e}_{\alpha\beta}{T}^{\alpha}d\theta^{\beta}$
$\displaystyle=$
$\displaystyle\int\limits_{\mathcal{C}_{\mathrm{p}}}{E}_{\alpha\beta}{T}^{\alpha}d\theta^{\beta},$
(B.3) $\displaystyle=$
$\displaystyle\int\limits_{\mathcal{C}_{\mathrm{m}}}\left[{E}_{\alpha\beta}\frac{d\theta^{\beta}}{dl}\right]{T}^{\alpha}dl,$

which on using Equation A.57 can be further wriiten as

$\int\limits_{\mathcal{C}_{\mathrm{p}}}{A}^{1/2}{e}_{\alpha\beta}{T}^{\alpha}d\theta^{\beta}=\int\limits_{\mathcal{C}_{\mathrm{m}}}{T}^{\alpha}{\nu}_{\alpha}dl.$
(B.4)

Therefore, we relate the integral over the parametric domain to the line-
integral along the boundary of the curved mid-surface.

##  Appendix A Variation of some relevant quantities

Here, we list the first variation of key kinematic variables (essential for
the calculations in Sec. 4), for example,

$\delta\bm{F}=\delta\frac{\partial\bm{\chi}}{\partial\bm{X}}=\frac{\partial\delta\bm{\chi}}{\partial\bm{X}},$
(C.1)

and following

$\bm{F}\bm{F}^{-1}=\mathbold{\mathbb{1}},$ (C.2)

we obtain

$\delta\bm{F}^{-1}=-\bm{F}^{-1}\delta\bm{F}\bm{F}^{-1}.$ (C.3)

Also,

$\delta J=J\bm{F}^{-T}\mathrel{\mathop{\ordinarycolon}}\delta\bm{F}.$ (C.4)

Again,

$\displaystyle\delta\bm{\chi}$ $\displaystyle=$
$\displaystyle\delta\bm{r}+\eta\delta\bm{d},$ (C.5) $\displaystyle=$
$\displaystyle\delta\bm{r}+\eta\delta\lambda\bm{n}+\eta\lambda\delta\bm{n},$

where $\delta\bm{n}$ can be obtained by using the relations,
$\bm{n}\cdot\bm{n}=1$ and $\bm{a}_{\alpha}\cdot\bm{n}=0$ as

$\delta\bm{n}=-\left[\bm{a}^{\alpha}\otimes\bm{n}\right]\delta\bm{a}_{\alpha}=-\bm{a}^{\alpha}\left[\bm{n}\cdot\delta\bm{a}_{\alpha}\right],$
(C.6)

with

$\delta\bm{a}_{\alpha}=\delta\frac{\partial\bm{r}}{\partial\theta^{\alpha}}=\left[\delta\bm{r}\right]_{,\alpha},$
(C.7)

and moreover, from Equation (3.13) follows

$\delta\lambda=-\frac{\lambda}{2}a^{-1}\delta a,$ (C.8)

where from Equation (A.16),

$\delta a=aa^{\alpha\beta}\delta a_{\alpha\beta},$ (C.9)

which can be rewritten by using,

$\delta
a_{\alpha\beta}=\delta{\bm{a}}_{\alpha}\cdot{\bm{a}}_{\beta}+{\bm{a}}_{\alpha}\cdot\delta{\bm{a}}_{\beta},$
(C.10)

as

$\delta a=2a\bm{a}^{\alpha}\cdot\delta{\bm{a}}_{\alpha}.$ (C.11)

Therefore, the variation in the through-thickness stretch can be rewritten as

$\delta\lambda=-\lambda\bm{a}^{\alpha}\cdot\delta{\bm{a}}_{\alpha}.$ (C.12)

Also,

$\delta p=\delta p_{0}+\eta\delta p_{1}+\mathcal{O}\left(\eta^{2}\right),$
(C.13)

and $p=p_{0}+\eta p_{1}+\mathcal{O}\left(\eta^{2}\right)$.

## References

  * Adkins and Rivlin (1952) Adkins J.E. and Rivlin R.S.  “Large elastic deformations of isotropic materials ix. the deformation of thin shells”.  _Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences_ , 244(888):505–531 (1952)
  * Akkas (1978) Akkas N.  “On the dynamic snap-out instability of inflated non-linear spherical membranes”.  _International Journal of Non-Linear Mechanics_ , 13(3):177–183 (1978)
  * Alzetta and Heltai (2020) Alzetta G. and Heltai L.  “Multiscale modeling of fiber reinforced materials via non-matching immersed methods”.  _Computers & Structures_, 239:106334 (2020)
  * Barham et al. (2008) Barham M., Steigmann D.J., McElfresh M., and Rudd R.E.  “Limit-point instability of a magnetoelastic membrane in a stationary magnetic field”.  _Smart Materials and Structures_ , 17(5):055003 (2008)
  * Basar and Krätzig (1985) Basar Y. and Krätzig W.  “Mechanik der flächentragwerke”.  _International Journal of Nonlinear Mechanics_ (1985)
  * Bischoff et al. (2004) Bischoff M., Bletzinger K., Wall W., and Ramm E.  “Models and finite elements for thin-walled structures”.  _Encyclopedia of Computational Mechanics_ (2004)
  * Blaauwendraad and Hoefakker (2014) Blaauwendraad J. and Hoefakker J.H.  “Structural shell analysis”.  _Solid Mechanics and Its Applications_ , 200 (2014)
  * Cheewaruangroj et al. (2019) Cheewaruangroj N., Leonavicius K., Srinivas S., and Biggins J.S.  “Peristaltic elastic instability in an inflated cylindrical channel”.  _Phys. Rev. Lett._ , 122:068003 (2019)
  * Cirak and Ortiz (2001) Cirak F. and Ortiz M.  “Fully ${C}^{1}$-conforming subdivision elements for finite deformation thin-shell analysis”.  _International Journal for Numerical Methods in Engineering_ , 51:813–833 (2001)
  * Cirak et al. (2000) Cirak F., Ortiz M., and Schröder P.  “Subdivision surfaces: a new paradigm for thin-shell finite-element analysis”.  _International Journal for Numerical Methods in Engineering_ , 47(12):2039–2072 (2000)
  * Dorfmann and Ogden (2014) Dorfmann L. and Ogden R.W.  _Nonlinear Theory of Electroelastic and Magnetoelastic Interactions._ Springer (2014)
  * Galley et al. (2019) Galley A., Knopf G.K., and Kashkoush M.  “Pneumatic hyperelastic actuators for grasping curved organic objects”.  _In_ “Actuators”, volume 8, page 76. MDPI (2019)
  * Ghosh et al. (2023) Ghosh A., McBride A., Liu Z., Heltai L., Steinmann P., and Saxena P.  “A fully-coupled nonlinear magnetoelastic thin shell formulation”.  _arXiv:2308.12300_ (2023)
  * Guo (2001) Guo X.  “Large deformation analysis for a cylindrical hyperelastic membrane of rubber-like material under internal pressure”.  _Rubber chemistry and technology_ , 74(1):100–115 (2001)
  * Hao et al. (2017) Hao Y., Wang T., Ren Z., Gong Z., Wang H., Yang X., Guan S., and Wen L.  “Modeling and experiments of a soft robotic gripper in amphibious environments”.  _International Journal of Advanced Robotic Systems_ , 14(3):1729881417707148 (2017)
  * Hart-Smith and Crisp (1967) Hart-Smith L. and Crisp J.  “Large elastic deformations of thin rubber membranes”.  _International Journal of Engineering Science_ , 5(1):1–24 (1967)
  * Haughton and Ogden (1979) Haughton D. and Ogden R.  “Bifurcation of inflated circular cylinders of elastic material under axial loading—ii. exact theory for thick-walled tubes”.  _Journal of the Mechanics and Physics of Solids_ , 27(5):489–512 (1979)
  * Heltai et al. (2021) Heltai L., Caiazzo A., and Müller L.O.  “Multiscale coupling of one-dimensional vascular models and elastic tissues”.  _Annals of Biomedical Engineering_ , 49:3243–3254 (2021)
  * Heltai and Zunino (2023) Heltai L. and Zunino P.  “Reduced lagrange multiplier approach for non-matching coupling of mixed-dimensional domains”.  _Mathematical Models and Methods in Applied Sciences_ , 33(12):2425–2462 (2023)
  * Hughes et al. (2005) Hughes T.J.R., Cottrell J.A., and Bazilevs Y.  “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement”.  _Computer Methods in Applied Mechanics and Engineering_ , 194(39):4135–4195 (2005)
  * Jolly et al. (1996) Jolly M.R., Carlson J.D., and Muñoz B.C.  “A model of the behaviour of magnetorheological materials”.  _Smart Materials and Structures_ , 5(5):607 (1996)
  * Khayat et al. (1992) Khayat R.E., Derdorri A., and García-Rejón A.  “Inflation of an elastic cylindrical membrane: non-linear deformation and instability”.  _International journal of solids and structures_ , 29(1):69–87 (1992)
  * Kiendl et al. (2009) Kiendl J., Bletzinger K.U., Linhard J., and Wüchner R.  “Isogeometric shell analysis with Kirchhoff–Love elements”.  _Computer Methods in Applied Mechanics and Engineering_ , 198(49):3902–3914 (2009)
  * Kiendl et al. (2015) Kiendl J., Hsu M.C., Wu M.C., and Reali A.  “Isogeometric kirchhoff–love shell formulations for general hyperelastic materials”.  _Computer Methods in Applied Mechanics and Engineering_ , 291:280–303 (2015)
  * Kim et al. (2012) Kim J., Hanna J.A., Byun M., Santangelo C.D., and Hayward R.C.  “Designing responsive buckled surfaces by halftone gel lithography”.  _Science_ (2012)
  * Liu et al. (2023) Liu Z., McBride A., Ghosh A., Heltai L., Huang W., Yu T., Steinmann P., and Saxena P.  “Computational instability analysis of inflated hyperelastic thin shells using subdivision surfaces”.  _Computational Mechanics_ (2023)
  * McKay et al. (2010) McKay T., O’Brien B., Calius E., and Anderson I.  “Self-priming dielectric elastomer generators”.  _Smart Materials and Structures_ , 19(5):055025 (2010)
  * Mehta et al. (2022) Mehta S., Raju G., Kumar S., and Saxena P.  “Instabilities in a compressible hyperelastic cylindrical channel under internal pressure and external constraints”.  _International Journal of Non-Linear Mechanics_ , 144:104031 (2022)
  * Niordson (1985) Niordson F.  “Shell theory”.  _North-Holland Series in Applied Mathematics and Mechanics: Elsevier Science_ (1985)
  * Ogden (1997) Ogden R.W.  _Non-linear elastic deformations_.  Dover Publications (1997)
  * Pamplona et al. (2006) Pamplona D., Goncalves P., and Lopes S.  “Finite deformations of cylindrical membrane under internal pressure”.  _International journal of mechanical sciences_ , 48(6):683–696 (2006)
  * Pelteret and Steinmann (2020) Pelteret J.P. and Steinmann P.  _Magneto-Active Polymers: Fabrication, Characterisation, Modelling and Simulation at the Micro- and Macro-Scale_.  DeGruyter (2020)
  * Radwańska et al. (2017) Radwańska M., Stankiewicz A., Wosatko A., and Pamin J.  “Plate and shell structures”.  _John Wiley and Sons_ (2017)
  * Reddy and Saxena (2017) Reddy N.H. and Saxena P.  “Limit points in the free inflation of a magnetoelastic toroidal membrane”.  _International Journal of Non-Linear Mechanics_ , 95:248–263 (2017)
  * Reddy and Saxena (2018) Reddy N.H. and Saxena P.  “Instabilities in the axisymmetric magnetoelastic deformation of a cylindrical membrane”.  _International Journal of Solids and Structures_ , 136-137:203–219 (2018)
  * Roohbakhshan and Sauer (2017) Roohbakhshan F. and Sauer R.A.  “Efficient isogeometric thin shell formulations for soft biological materials”.  _Biomechanics and Modeling in Mechanobiology_ , 16(5):1569–1597 (2017)
  * Saxena et al. (2019) Saxena P., Reddy N.H., and Pradhan S.P.  “Magnetoelastic deformation of a circular membrane: Wrinkling and limit point instabilities”.  _International Journal of Non-Linear Mechanics_ , 116:250–261 (2019)
  * Sheng et al. (2012) Sheng J., Chen H., Qiang J., Li B., and Wang Y.  “Thermal, mechanical, and dielectric properties of a dielectric elastomer for actuator applications”.  _Journal of Macromolecular Science, Part B_ , 51(10):2093–2104 (2012)
  * Simo and Fox (1989) Simo J.C. and Fox D.  “On a stress resultant geometrically exact shell model. Part I: formulation and optimal parametrization”.  _Applied Mechanics and Engineering_ , 72:267–304 (1989)
  * Steinmann (2015) Steinmann P.  “Geometrical foundations of continuum mechanics”.  _LAMM, volume 2_ (2015)
  * Tamadapu and DasGupta (2013) Tamadapu G. and DasGupta A.  “Finite inflation analysis of a hyperelastic toroidal membrane of initially circular cross-section”.  _International Journal of Non-Linear Mechanics_ , 49:31–39 (2013)
  * Tepole et al. (2015) Tepole A.B., Kabaria H., Bletzinger K.U., and Kuhl E.  “Isogeometric kirchhoff–love shell formulations for biological membranes”.  _Computer methods in applied mechanics and engineering_ , 293:328–347 (2015)
  * Venkata and Saxena (2019) Venkata S.P. and Saxena P.  “Instabilities in the free inflation of a nonlinear hyperelastic toroidal membrane”.  _Journal of Mechanics of Materials and Structures_ , 14(4):473–496 (2019)
  * Verron et al. (1999) Verron E., Khayat R., Derdouri A., and Peseux B.  “Dynamic inflation of hyperelastic spherical membranes”.  _Journal of Rheology_ , 43(5):1083–1097 (1999)
  * Wempner and Talaslidis (2002) Wempner G. and Talaslidis D.  “Mechanics of solids and shells: Theories and approximations”.  _CRC Press_ (2002)
  * Xie et al. (2016) Xie Y.X., Liu J.C., and Fu Y.  “Bifurcation of a dielectric elastomer balloon under pressurized inflation and electric actuation”.  _International Journal of Solids and Structures_ , 78-79:182–188 (2016)
  * Yan et al. (2020) Yan D., Pezzulla M., Cruveiller L., Abbasi1 A., and Reis P.M.  “Magneto-active elastic shells with tunable buckling strength”.  _Nature Communications_ , 12:2831 (2020)
  * Yang and Feng (1970) Yang W. and Feng W.  “On axisymmetrical deformations of nonlinear membranes”.  _Journal of Applied Mechanics_ , 37(4):1002–1011 (1970)