A Geometrically Exact Nonlinear Kirchhoff–Love Shell Theory in the Language of Fibre Bundles, Moving Frames and Differential Forms

1. Introduction

Shell theory occupies an awkward position between three-dimensional continuum mechanics and the geometry of surfaces. The classical engineering literature [16,30,34,52] reaches the field equations by a sequence of coordinate manipulations in which the geometric content—what is a tensor, what is a connection, what is merely a coordinate artefact—is progressively obscured. Two costs follow. First, the compatibility conditions of Gauss and Codazzi–Mainardi have to be imposed by hand as extra equations, when they are in fact identities. Second, and more seriously for computation, the strain measures produced by displacement-based derivations are frequently only approximately invariant under rigid rotation, so that a rotating but unstrained shell element registers a fictitious strain. In an explicit dynamic computation, or in a buckling analysis where the physical membrane strains are themselves of order ${10}^{-3}$, this is not a small matter: Section 10.1 shows the effect is measured in hundreds of megapascals.

Both costs disappear if one takes the geometry seriously from the start. The observation organising this paper is that a shell is a fibre bundle. The mid-surface $\mathcal{M}$ is the base, the through-thickness coordinate parameterises the fibre $I=[-h/2,h/2]$, and the shell body $\mathcal{S}$ is the total space. Everything then falls into place:

• Kinematics is the pull-back of the ambient Euclidean metric along the deformation. Because a pull-back of a metric is manifestly invariant under isometries of the target, objectivity is automatic rather than a property to be checked.
• Geometry of the reference state is encoded in an orthonormal moving coframe $\left\{{\Theta }^I\right\}$ and its connection 1-forms ${{\omega }^I}_J$. Cartan’s structure equations [14] then deliver the Gauss and Codazzi–Mainardi relations as the statement that a certain curvature 2-form vanishes, because the reference configuration is embedded in flat ${\mathbb{R}}^3$.
• Stress is a vector-bundle-valued 2-form, its exterior covariant derivative is the divergence, and Cauchy’s theorem is Stokes’ theorem.
• Dimensional reduction is integration over the fibre. The resultants $N^{ab}$ and $M^{ab}$ are the zeroth and first fibre moments, and no assumption beyond the Kirchhoff–Love constraint is required to obtain them.

• What is new here.

The bundle picture is not itself new. Intrinsic, coordinate-free formulations of shell theory go back to Koiter [25,26], Naghdi [34] and Sanders [39]; the geometrically exact stress-resultant models of Simo, Fox and Rifai placed the theory squarely in the language of configuration manifolds, with the configuration space taken to be the manifold of maps into ${\mathbb{R}}^3\times S^2$ and the rotation of the director realised as motion on the unit sphere: Part I [41,48] formulates the model and its parametrisation—optimal in the precise sense that no smaller set of parameters describes the director field without introducing singularities ([41], §1.2(i))— ([42], Part II) develops the linear theory and the mixed interpolation, and Part III [43] carries both into the finite-deformation range; Libai and Simmonds [28] and Antman [2] give the definitive treatments of the nonlinear theory; and Ciarlet [9] and Marsden and Hughes [33] supply the differential-geometric apparatus. Rigorous justifications of two-dimensional models as limits of three-dimensional elasticity are given by Fox, Raoult and Simo [15] and Steigmann [47].

The contribution of this paper is to carry it through in complete detail, to correct a family of errors that recur in presentations of this kind, and to close the loop with a verified computation. Specifically, we draw attention to the following points, each of which is a place where published derivations are commonly wrong or incomplete.

(i)

The variation of the bending strain is not ${\partial }_a\delta \mathit{\phi }·{\partial }_b\mathit{n}+{\partial }_a\mathit{\phi }·{\partial }_b\delta \mathit{n}$ with $\delta \mathit{n}$ free; under the Kirchhoff–Love constraint $\mathit{n}$ is a function of $\mathit{\phi }$, and the correct result (Proposition 4) is $\delta K_{ab}=\mathit{n}·{\stackrel{\circ }{\nabla }}_a{\stackrel{\circ }{\nabla }}_b\delta \mathit{\phi }$ with the current-surface connection.

(ii)

The linear-momentum balance is not ${\tilde{\nabla }}_b{N}^{ab}+{\overline{n}}^a={\rho }_{0}h{\ddot{\phi }}^a$. The correct tangential equation carries a shear–curvature coupling term and the effective resultant ${\widehat{N}}^{ab}=N^{ab}+M^{ac}{{\tilde{b}}_c}^b$ (Section 5.4). The necessity of an effective resultant is classical. Simo and Fox ([41], §1.2(viii), Eqs. (4.24)–(4.27)) obtain the companion result for the one-director Cosserat surface: balance of angular momentum forces the combination ${\tilde{n}}^{\beta \alpha }=n^{\beta \alpha }-{\lambda }_{\mu }^{\beta }{\tilde{m}}^{\alpha \mu }$ to be symmetric while $n^{\beta \alpha }$ itself is not, and it is this effective resultant, not $n^{\beta \alpha }$, that is conjugate to the membrane strain rate in the reduced stress power ([41], Cor. 5.2). Part II states the same conclusion in the form used below: the constitutive restriction imposed by angular-momentum balance ([42], Equation (2.4)) defines the effective resultants ([42], Equation (2.5)), and only the symmetric effective membrane resultant survives in the weak form ([42], Equation (2.12)). The resultants themselves cannot be symmetric in the presence of initial curvature [34,35], ([42], §1(i)). The term is nevertheless routinely dropped in derivations that proceed by analogy with plate theory.

(iii)

Because the shell theory is of second differential order in $\mathit{\phi }$, the boundary term does not separate into independent shear and twist; it must be integrated by parts along $\partial \mathcal{M}$, producing the Kirchhoff effective shear and concentrated corner forces (Section 6).

(iv)

Fibre integration of the exact volume form ${\mu }_0\left(z\right)\mathrm{d}z\wedge \mathrm{d}A$ with ${\mu }_0=1+2Hz+K{z}^2$ produces a membrane–bending coupling stiffness $\propto H{h}^3$ and corrections $\propto K{h}^2$ to the extensional and flexural stiffnesses, which vanish only in the limit $h/R\to 0$ (Section 5). Such curvature-induced terms are not new—they are given by Simo, Fox and Rifai ([42], Table 2) as a “higher-order constitutive model” obtained by asymptotic expansion—but the route here is different, and comparing the two exposes a consistency requirement that is easy to miss (Section 7.2).

(v)

Elementary but consequential: the resultant constitutive law is $N^{ab}=h{H}^{abcd}{E}_{cd}^{\left(0\right)}$, not $C{H}^{abcd}{E}_{cd}^{\left(0\right)}$ with $C=Eh/(1-{\nu }^2)$ and $H^{abcd}$ already containing $E/(1-{\nu }^2)$; the latter double-counts the modulus.

• Relation to computational practice.

The $C^1$ regularity that the Kirchhoff–Love constraint demands (Remark 4) has driven much of the recent finite-element literature: subdivision-surface elements [10,11] and isogeometric analysis [19,23] both exist largely to supply it. The alternative is to give up $C^1$ altogether by carrying an independent director with its own rotational degrees of freedom—the route taken by the degenerated-solid elements descending from Ahmad, Irons and Zienkiewicz [1] and by the one-director Cosserat models of [41,42,43]. The price is a shear-flexible kinematics with five degrees of freedom per node and, with it, the problem of updating a rotation exactly. Part III [43] settles that problem in closed form: since the director is a point of $S^2$, its update is the exponential map along a great circle, which is well conditioned for increments of arbitrary magnitude and needs none of the quaternion machinery required for rods [45]. The discrete version of that map is already in Part II, where the nodal director interpolation is built so that the interpolated director has unit length and the incremental rotation stays orthogonal to it throughout the element ([42], Props. 4.1–4.3); the alternative of interpolating a nodal rotation vector is rejected there precisely because extracting it requires careful use of quaternions to avoid singularities ([42], Rem. 4.4(2)).

It is worth being explicit that this is the opposite of the choice made here. Simo and Fox parametrise the weak form so that the Riemannian connection of the mid-surface never appears, on the grounds that Christoffel symbols and the second fundamental form are not readily available in a computational setting ([41], §1.2(vi)); Part II makes the point explicitly, noting that Christoffel symbols, the second fundamental form and covariant derivatives never appear in the field equations as they cast them ([42], §1(iii)). The present paper does the reverse, making the connection 1-forms the carriers of all metric information (Section 2). Neither choice is more correct; they optimise for different things—theirs for an element that can be coded without differential geometry, ours for a formulation in which compatibility is an identity rather than an imposed constraint. Section 10.10 compares the two where they overlap; the difference in kinematic hypothesis—shear-rigid here, shear-flexible there—is what fixes the sign of the residual discrepancy we find. Surveys of the competing discretisation strategies are given by Bischoff, Ramm and Irslinger [4] and Bathe [3]; standard benchmark problems for geometrically nonlinear shells are collected by Sze, Liu and Lo [49]. The present paper is not a finite-element paper, but Section 10 shows that the formulation discretises cleanly, that its theoretical guarantees survive discretisation, and that the consistent tangent it supplies delivers the asymptotically quadratic Newton convergence that such a tangent is supposed to deliver (Section 10.11).

• Outline.

Section 2 sets up the bundle, the moving frame and the structure equations. Section 3 develops the finite-deformation kinematics and the strain measures, with the objectivity proof and the exact truncation estimate. Section 4 introduces the bundle-valued stress forms and the three-dimensional virtual-work statement. Section 5 performs the fibre integration and derives the strong form; Section 6 the boundary conditions. Section 7 closes the system constitutively, Section 8 verifies that the classical theories are recovered as limits, and Section 9 gives the consistent tangent operator required for Newton solution. Section 10 is the numerical case study.

2. Geometric Foundations

2.1. The Shell as a Fibre Bundle

Definition 1 

(Shell bundle). Let $\mathcal{M}$ be a compact orientable 2-manifold with piecewise-smooth boundary, embedded in ${\mathbb{R}}^3$ by ${\mathit{\phi }}_0:\mathcal{M}\to {\mathbb{R}}^3$ with unit normal field $\mathit{N}$. Let $I=[-h/2,h/2]$. The shell bundle is the trivial fibre bundle

$$\pi :\mathcal{S}=\mathcal{M}\times I⟶\mathcal{M},\pi (\mathit{x},z)=\mathit{x},$$

(1)

(Figure 1) with reference embedding

$${\mathbf{\Phi }}_0:\mathcal{S}\to {\mathbb{R}}^3,{\mathbf{\Phi }}_0(\mathit{x},z)={\mathit{\phi }}_0\left(\mathit{x}\right)+z\mathit{N}\left(\mathit{x}\right).$$

(2)

Figure 1. The shell as a fibre bundle. Left: the mid-surface $\mathcal{M}$ is the base, the fibres (red) are the through-thickness segments $I=[-h/2,h/2]$ carried along the unit normal, and the total space $\mathcal{S}=\mathcal{M}\times I$ is the shell body. Right: the orthonormal moving frame $\{{\mathit{E}}_1,{\mathit{E}}_2,{\mathit{E}}_3=\mathit{N}\}$ dual to the coframe (5), in which the pulled-back metric is $\mathit{G}={\delta }_{IJ}{\Theta }^I\otimes {\Theta }^J$ and all geometric information is carried by the connection 1-forms ${{\omega }^I}_J$.

Figure 1. The shell as a fibre bundle. Left: the mid-surface $\mathcal{M}$ is the base, the fibres (red) are the through-thickness segments $I=[-h/2,h/2]$ carried along the unit normal, and the total space $\mathcal{S}=\mathcal{M}\times I$ is the shell body. Right: the orthonormal moving frame $\{{\mathit{E}}_1,{\mathit{E}}_2,{\mathit{E}}_3=\mathit{N}\}$ dual to the coframe (5), in which the pulled-back metric is $\mathit{G}={\delta }_{IJ}{\Theta }^I\otimes {\Theta }^J$ and all geometric information is carried by the connection 1-forms ${{\omega }^I}_J$.

Remark 1. 

${\mathbf{\Phi }}_0$ is an embedding provided $\left|z\right|max\left\{\right|{\kappa }_1|,|{\kappa }_2\left|\right\}<1$ pointwise, i.e., $h<2{min}_{\mathcal{M}}\left\{\right|{\kappa }_1{|}^{-1},|{\kappa }_2{|}^{-1}\}$: the thickness must not reach the focal surfaces. All statements below are made under this hypothesis, which is the precise form of the “thin shell” assumption at the level of geometry. It is weaker than the constitutive thinness assumption $h/R\ll 1$ used later.

Triviality of the bundle is a convenience, not a restriction of substance: for a non-orientable mid-surface one works with the associated ${\mathbb{Z}}_2$-bundle, and every formula below is local and therefore unaffected.

2.2. Moving Coframe and Connection Forms

Let $({\alpha }^1,{\alpha }^2)$ be lines-of-curvature coordinates on $\mathcal{M}$ with Lamé parameters $A_1,A_2$ and principal radii $R_1,R_2$, so that the mid-surface first and second fundamental forms are

$$\mathit{a}=\sum _a{A}_a^2\mathrm{d}{\alpha }^a\otimes \mathrm{d}{\alpha }^a,\mathit{d}=\sum _a{\displaystyle \frac{A_a^2}{R_a}}\mathrm{d}{\alpha }^a\otimes \mathrm{d}{\alpha }^a.$$

(3)

Remark 2 

(Sign convention). We fix once and for all

$$d_{ab}:={\partial }_a{\mathit{\phi }}_0·{\partial }_b\mathit{N}=-\mathit{N}·{\partial }_a{\partial }_b{\mathit{\phi }}_0,$$

(4)

so that $\mathit{d}={\sum }_a(A_a^2/R_a)\mathrm{d}{\alpha }^a\otimes \mathrm{d}{\alpha }^a$ with $R_a>0$ when the centre of curvature lies on the $-\mathit{N}$ side. This is the convention forced on us by the metric expansion (18), and it is the opposite of the more common $b_{ab}=\mathit{N}·{\partial }_a{\partial }_b{\mathit{\phi }}_0$. Mixing the two is a frequent source of sign errors in the bending strain; we therefore never use the symbol b for a reference quantity, and write ${\tilde{b}}_{ab}:=\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi }=-d_{t,ab}$ only where the outward-normal convention is genuinely wanted, namely in the normal equilibrium equation (40).

Define on $\mathcal{S}$ the orthonormal coframe

$${\Theta }^1=A_1\left(1+{\displaystyle \frac{z}{R_1}}\right)\mathrm{d}{\alpha }^1,{\Theta }^2=A_2\left(1+{\displaystyle \frac{z}{R_2}}\right)\mathrm{d}{\alpha }^2,{\Theta }^3=\mathrm{d}z,$$

(5)

so that the pull-back of the Euclidean metric along ${\mathbf{\Phi }}_0$ is

$$\mathit{G}={\mathbf{\Phi }}_0^{*}\mathbf{\delta }={\delta }_{IJ}{\Theta }^I\otimes {\Theta }^J,I,J\in \{1,2,3\},$$

(6)

and the Riemannian volume form is

$$\mathrm{d}V={\Theta }^1\wedge {\Theta }^2\wedge {\Theta }^3={\mu }_0\left(z\right)\mathrm{d}A\wedge \mathrm{d}z,\begin{array}{|c|}\hline {\mu }_0\left(z\right)=\left(1+{\displaystyle \frac{z}{R_1}}\right)\left(1+{\displaystyle \frac{z}{R_2}}\right)=1+2Hz+K{z}^2\\ \hline\end{array}$$

(7)

with $H=\frac{1}{2}(R_1^{-1}+R_2^{-1})$ the mean and $K={\left(R_1{R}_2\right)}^{-1}$ the Gaussian curvature. The function ${\mu }_0$ is the determinant of the shifter $\mathit{Z}=\mathbf{1}+z{\mathit{d}}^{♯}$ that maps the mid-surface tangent plane onto the tangent plane of the parallel surface at height z; it is the exact Jacobian of the fibre coordinates and is retained throughout Section 5. Shifter tensors of this kind are standard in the shell literature [16,34] but are usually set to unity at the first opportunity.

Proposition 1 

(Connection forms). The unique torsion-free metric connection 1-forms ${{\omega }^I}_J=-{{\omega }^J}_I$ satisfying the first structure equation $\mathrm{d}{\Theta }^I+{{\omega }^I}_J\wedge {\Theta }^J=0$ are

$${{\omega }^a}_3={\displaystyle \frac{{\Theta }^a}{R_a+z}}(a=1,2),{{\omega }^1}_2={\displaystyle \frac{1}{{\overline{A}}_1{\overline{A}}_2}}\left({\displaystyle \frac{\partial {\overline{A}}_1}{\partial {\alpha }^2}}{\Theta }^1-{\displaystyle \frac{\partial {\overline{A}}_2}{\partial {\alpha }^1}}{\Theta }^2\right),$$

(8)

where ${\overline{A}}_a:=A_a(1+z/R_a)$ and no summation is implied in ${{\omega }^a}_3$.

Proof. 

Antisymmetry is metric compatibility of (6). For $I=3$: $\mathrm{d}{\Theta }^3=0$ and ${{\omega }^3}_a\wedge {\Theta }^a=-{(R_a+z)}^{-1}{\Theta }^a\wedge {\Theta }^a=0$, so the third structure equation holds identically. For $I=1$,

$$\mathrm{d}{\Theta }^1={\displaystyle \frac{\partial {\overline{A}}_1}{\partial {\alpha }^2}}\mathrm{d}{\alpha }^2\wedge \mathrm{d}{\alpha }^1+{\displaystyle \frac{A_1}{R_1}}\mathrm{d}z\wedge \mathrm{d}{\alpha }^1,$$

while ${{\omega }^1}_3\wedge {\Theta }^3={\displaystyle \frac{{\overline{A}}_1}{R_1+z}}\mathrm{d}{\alpha }^1\wedge \mathrm{d}z={\displaystyle \frac{A_1}{R_1}}\mathrm{d}{\alpha }^1\wedge \mathrm{d}z$, which cancels the second term exactly. The remaining requirement ${\displaystyle \frac{\partial {\overline{A}}_1}{\partial {\alpha }^2}}\mathrm{d}{\alpha }^2\wedge \mathrm{d}{\alpha }^1+{{\omega }^1}_2\wedge {\Theta }^2=0$ fixes the ${\Theta }^1$-component of ${{\omega }^1}_2$ as stated, and the companion equation for ${\Theta }^2$ fixes the ${\Theta }^2$-component. Uniqueness is the fundamental lemma of Riemannian geometry.    □

2.3. Curvature 2-Form; Gauss and Codazzi–Mainardi

Theorem 1 

(Compatibility as flatness). Because ${\mathbf{\Phi }}_0\left(\mathcal{S}\right)\subset {\mathbb{R}}^3$ and ${\mathbb{R}}^3$ is flat, the curvature 2-form

$${{\mathsf{\Omega }}^I}_J=\mathrm{d}{{\omega }^I}_J+{{\omega }^I}_L\wedge {{\omega }^L}_J$$

(9)

vanishes identically. Its components are equivalent to the classical compatibility conditions of surface theory:

$$\begin{array}{cc}\hfill {{\mathsf{\Omega }}^1}_3=0& \iff {\displaystyle \frac{\partial }{\partial {\alpha }^2}}\left({\displaystyle \frac{A_1}{R_1}}\right)={\displaystyle \frac{1}{R_2}}{\displaystyle \frac{\partial A_1}{\partial {\alpha }^2}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {{\mathsf{\Omega }}^2}_3=0& \iff {\displaystyle \frac{\partial }{\partial {\alpha }^1}}\left({\displaystyle \frac{A_2}{R_2}}\right)={\displaystyle \frac{1}{R_1}}{\displaystyle \frac{\partial A_2}{\partial {\alpha }^1}},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {{\mathsf{\Omega }}^1}_2=0& \iff {\displaystyle \frac{\partial }{\partial {\alpha }^1}}\left({\displaystyle \frac{1}{A_1}}{\displaystyle \frac{\partial A_2}{\partial {\alpha }^1}}\right)+{\displaystyle \frac{\partial }{\partial {\alpha }^2}}\left({\displaystyle \frac{1}{A_2}}{\displaystyle \frac{\partial A_1}{\partial {\alpha }^2}}\right)=-{\displaystyle \frac{A_1{A}_2}{R_1{R}_2}}.\hfill \end{array}$$

(12)

Equations (10)–(11) are the Codazzi–Mainardi equations and (12) is Gauss’ theorema egregium [13], $K={\left(R_1{R}_2\right)}^{-1}$.

Proof. 

${\mathbf{\Phi }}_0$ is an isometric embedding into flat space, so the induced Levi-Civita connection has vanishing curvature; (9) is the second Cartan structure equation. Substituting (8) and collecting the coefficients of ${\Theta }^a\wedge {\Theta }^3$ and ${\Theta }^1\wedge {\Theta }^2$ gives the stated component identities; note that the z-dependence cancels identically, which is precisely the statement that the parallel surfaces of a surface satisfying (10)–(12) again satisfy them.    □

Remark 3. 

This is the first concrete payoff of the formalism. In a coordinate treatment, (10)–(12) are three additional partial differential equations that must be appended to the field equations and then shown to be preserved. Here they are consequences of the ambient space being flat, and they are automatically satisfied by any configuration that is actually realised as a surface in ${\mathbb{R}}^3$—which is exactly the class of configurations a displacement-based formulation produces.

2.4. Exterior Covariant Derivative

For a ${\mathbb{R}}^3$-valued p-form $\mathit{\eta }={\eta }^I{\mathit{e}}_I$ on $\mathcal{S}$ (components in the moving frame ${\mathit{e}}_I$ dual to ${\Theta }^I$) the exterior covariant derivative is

$${\mathrm{d}}^{\nabla }\mathit{\eta }=\left(\mathrm{d}{\eta }^I+{{\omega }^I}_J\wedge {\eta }^J\right){\mathit{e}}_I.$$

(13)

Theorem 1 gives ${\mathrm{d}}^{\nabla }\circ {\mathrm{d}}^{\nabla }=0$ on $\mathcal{S}$, so the complex of bundle-valued forms is exact and Stokes’ theorem holds in the form

$${\int }_{\mathcal{S}}{\mathrm{d}}^{\nabla }\mathit{\eta }={\int }_{\partial \mathcal{S}}\mathit{\eta }.$$

(14)

This identity is what replaces the divergence theorem in Section 4, and it is the reason no Christoffel symbols appear explicitly anywhere in the derivation of the balance laws.

3. Finite-Deformation Kinematics

3.1. Deformation and the Kirchhoff–Love Constraint

Let $\mathit{\phi }:\mathcal{M}\to {\mathbb{R}}^3$ be the deformed mid-surface placement, assumed an immersion, and let

$$\mathit{n}={\displaystyle \frac{{\partial }_1\mathit{\phi }\times {\partial }_2\mathit{\phi }}{|{\partial }_1\mathit{\phi }\times {\partial }_2\mathit{\phi }|}}$$

(15)

be the deformed unit normal. The shell deformation is

$$\mathbf{\Phi }:\mathcal{S}\to {\mathbb{R}}^3,\mathbf{\Phi }(\mathit{x},z)=\mathit{\phi }\left(\mathit{x}\right)+z\mathit{n}\left(\mathit{x}\right).$$

(16)

Remark 4 

(On the status of $\mathit{n}$). It is common to see $\mathit{n}$ described as “an independent kinematic variable subject to the orthogonality constraint $\mathit{n}·{\partial }_a\mathit{\phi }=0$”. These two statements are incompatible: once the constraint is imposed, $\mathit{n}$ is determined by $\mathit{\phi }$ through (15) and is not independent. Two consistent formulations exist [2,41]:

(a) 

the constrained (Kirchhoff–Love) theory used here, with the single field $\mathit{\phi }\in H^2(\mathcal{M};{\mathbb{R}}^3)$ and $\mathit{n}$ given by (15); or

(b) 

a mixed theory with independent director $\mathit{n}$ and Lagrange multipliers ${\lambda }^a$ enforcing $\mathit{n}·{\partial }_a\mathit{\phi }=0$, whose stationarity conditions recover (a) after elimination, and which relaxes to Reissner–Mindlin kinematics if the constraint is penalised [4]. The one-director Cosserat models of Simo, Fox and Rifai take the unconstrained limit of (b) as their starting point: the director is a free point of $S^2$ with two rotational degrees of freedom, so the kinematics is shear-flexible from the outset and no constraint is imposed ([42], §1(v)).

The distinction is not pedantic: it dictates the regularity required of the trial space ($C^1$ versus $C^0$) and it changes the variation of the bending strain, Proposition 4.

3.2. Pull-Back Metric: Exact Expansion

Write ${\mathit{a}}_t={\mathit{\phi }}^{*}\mathit{\delta }$ for the deformed first fundamental form, $a_{t,ab}={\partial }_a\mathit{\phi }·{\partial }_b\mathit{\phi }$, and, consistently with Remark 2,

$$d_{t,ab}:={\partial }_a\mathit{\phi }·{\partial }_b\mathit{n}=-\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi }.$$

(17)

Proposition 2 

(Exact pull-back). With Φ as in (16),

$${\mathbf{\Phi }}^{*}\mathit{\delta }={\mathit{a}}_t+2z{\mathit{d}}_t+z^2{\mathit{c}}_t+\mathrm{d}z\otimes \mathrm{d}z,c_{t,ab}:={\partial }_a\mathit{n}·{\partial }_b\mathit{n},$$

(18)

and correspondingly ${\mathbf{\Phi }}_0^{*}\mathit{\delta }=\mathit{a}+2z\mathit{d}+z^2\mathit{c}+\mathrm{d}z\otimes \mathrm{d}z$ with $c_{ab}={\partial }_a\mathit{N}·{\partial }_b\mathit{N}$. The tensor ${\mathit{c}}_t$ is the third fundamental form and satisfies the Cayley–Hamilton identity ${\mathit{c}}_t=2H_t{\mathit{d}}_t-K_t{\mathit{a}}_t$.

Proof. 

${\partial }_a\mathbf{\Phi }={\partial }_a\mathit{\phi }+z{\partial }_a\mathit{n}$ and ${\partial }_z\mathbf{\Phi }=\mathit{n}$. Taking inner products and using $\mathit{n}·{\partial }_a\mathit{\phi }=0$, $\mathit{n}·{\partial }_a\mathit{n}=0$, $\mathit{n}·\mathit{n}=1$ gives (18) with no remainder. The last identity is Cayley–Hamilton for the $2\times 2$ Weingarten map ${\mathit{d}}_t^{♯}$.    □

Note that (18) is exact: no expansion has been performed. This matters for the error estimate below.

3.3. Strain Measures

Definition 2 

(Membrane and bending strains).

$$\begin{array}{|c|}\hline {\mathit{E}}^{\left(0\right)}:=\frac{1}{2}\left({\mathit{a}}_t-\mathit{a}\right),{\mathit{K}}^{\left(1\right)}:={\mathit{d}}_t-\mathit{d}\\ \hline\end{array}$$

(19)

These are the finite-strain counterparts of Koiter’s membrane and bending measures [25,39]; see also [36,40] for equivalent finite-rotation forms. Both are covariant 2-tensor fields on $\mathcal{M}$, i.e., sections of ${\mathrm{Sym}}^2{T}^{*}\mathcal{M}$; in particular they are ordinary tensors, not “pseudo-strains”, and they may be pushed forward, pulled back and differentiated covariantly without qualification.

Proposition 3 

(Green–Lagrange strain of the shell body). The three-dimensional Green–Lagrange strain $\mathit{E}\left(z\right):=\frac{1}{2}({\mathbf{\Phi }}^{*}\mathit{\delta }-{\mathbf{\Phi }}_0^{*}\mathit{\delta })$ has components

$$E_{ab}\left(z\right)=E_{ab}^{\left(0\right)}+z{K}_{ab}^{\left(1\right)}+\frac{z^2}{2}\left(c_{t,ab}-c_{ab}\right),E_{a3}=E_{33}=0.$$

(20)

The vanishing of the transverse strains is the Kirchhoff–Love hypothesis; it is a consequence of (16), not a separate assumption.

Corollary 1 

(Truncation error). The thin-shell strain field used in Section 5 is $\mathit{E}\simeq {\mathit{E}}^{\left(0\right)}+z{\mathit{K}}^{\left(1\right)}$, with pointwise error

$$\parallel \mathit{E}\left(z\right)-\left({\mathit{E}}^{\left(0\right)}+z{\mathit{K}}^{\left(1\right)}\right)\parallel =\frac{z^2}{2}\parallel {\mathit{c}}_t-\mathit{c}\parallel =O\left(z^2\kappa \parallel {\mathit{K}}^{\left(1\right)}\parallel \right)=O\left(h^2/R^2\right)·\parallel \mathit{E}\parallel ,$$

(21)

uniformly on $\mathcal{S}$. Retaining ${\mathit{c}}_t$ is therefore consistent only if one also retains $O(h^2/R^2)$ constitutive terms; we do so in Section 5 through ${\mu }_0$ but discard the strain remainder, which is the standard—and here explicitly quantified—Kirchhoff–Love truncation [26,47]. Section 10.2 confirms the quadratic rate numerically.

3.4. Exact Objectivity

Theorem 2 

(Invariance under the Euclidean group). Let $\mathit{Q}\in \mathrm{SO}\left(3\right)$ and $\mathit{c}\in {\mathbb{R}}^3$ be constants and let ${\mathit{\phi }}^{†}=\mathit{Q}\mathit{\phi }+\mathit{c}$. Then

$${\mathit{E}}^{\left(0\right)}\left[{\mathit{\phi }}^{†}\right]={\mathit{E}}^{\left(0\right)}\left[\mathit{\phi }\right],{\mathit{K}}^{\left(1\right)}\left[{\mathit{\phi }}^{†}\right]={\mathit{K}}^{\left(1\right)}\left[\mathit{\phi }\right],$$

(22)

identically—not to any order in a rotation parameter.

Proof. 

${\partial }_a{\mathit{\phi }}^{†}=\mathit{Q}{\partial }_a\mathit{\phi }$, hence $a_{t,ab}^{†}=\mathit{Q}{\partial }_a\mathit{\phi }·\mathit{Q}{\partial }_b\mathit{\phi }={\partial }_a\mathit{\phi }·{\partial }_b\mathit{\phi }=a_{t,ab}$ because ${\mathit{Q}}^{\top }\mathit{Q}=\mathit{1}$. Since $det\mathit{Q}=+1$, (15) gives ${\mathit{n}}^{†}=\mathit{Q}\mathit{n}$, whence $d_{t,ab}^{†}=\mathit{Q}{\partial }_a\mathit{\phi }·\mathit{Q}{\partial }_b\mathit{n}=d_{t,ab}$. Both reference quantities are untouched.    □

Remark 5 

(Why this is not automatic). Linear and moderate-rotation shell theories [26,52] replace ${\mathit{E}}^{\left(0\right)}$ by its first variation at the reference configuration, ${\epsilon }_{ab}^{\mathrm{lin}}\left[\mathit{u}\right]=sym({\partial }_a\mathit{u}·{\partial }_b{\mathit{\phi }}_0)$. For a rigid rotation $\mathit{u}=(\mathit{Q}-\mathit{1}){\mathit{\phi }}_0+\mathit{c}$ this gives

$${\epsilon }_{ab}^{\mathrm{lin}}=\frac{1}{2}{\partial }_a{\mathit{\phi }}_0·\left(\mathit{Q}+{\mathit{Q}}^{\top }-2\mathit{1}\right){\partial }_b{\mathit{\phi }}_0=-(1-cos\theta ){\left({\mathit{P}}_{k^{\perp }}\right)}_{ab}\ne 0,$$

(23)

where θ is the rotation angle and ${\mathit{P}}_{k^{\perp }}$ the projector orthogonal to the rotation axis. The spurious strain therefore grows like ${\theta }^2/2$: $3.8\times {10}^{-5}$ at $0.5^{\circ }$, $3.0\times {10}^{-3}$ at $5^{\circ }$, $0.23$ at ${45}^{\circ }$. Section 10.1 confirms (23) numerically and converts it into fictitious stress.

3.5. Admissible Variations

Let $\delta \mathit{\phi }\in H^2(\mathcal{M};{\mathbb{R}}^3)$ be a variation. Because $\mathit{n}$ is slaved to $\mathit{\phi }$, its variation is determined.

Lemma 1 

(Variation of the normal).

$$\delta \mathit{n}=-\left({\partial }_c\delta \mathit{\phi }·\mathit{n}\right)a_t^{cd}{\partial }_d\mathit{\phi }.$$

(24)

Proof. 

$\left|\mathit{n}\right|=1$ gives $\mathit{n}·\delta \mathit{n}=0$, so $\delta \mathit{n}$ is tangential, $\delta \mathit{n}={\beta }^d{\partial }_d\mathit{\phi }$. Varying $\mathit{n}·{\partial }_c\mathit{\phi }=0$ gives $\delta \mathit{n}·{\partial }_c\mathit{\phi }=-\mathit{n}·{\partial }_c\delta \mathit{\phi }$, i.e., ${\beta }^d{a}_{t,dc}=-\mathit{n}·{\partial }_c\delta \mathit{\phi }$. Invert.    □

Proposition 4 

(Strain variations). With $\stackrel{\circ }{\nabla }$ the Levi-Civita connection of the deformed metric ${\mathit{a}}_t$ (Christoffel symbols ${\Gamma }_{ab}^c$),

$$\begin{array}{cc}\hfill \delta E_{ab}^{\left(0\right)}& =sym\left({\partial }_a\delta \mathit{\phi }·{\partial }_b\mathit{\phi }\right),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \delta K_{ab}^{\left(1\right)}& =-\mathit{n}·\left({\partial }_a{\partial }_b\delta \mathit{\phi }-{\Gamma }_{ab}^c{\partial }_c\delta \mathit{\phi }\right)=-\mathit{n}·{\stackrel{\circ }{\nabla }}_a{\stackrel{\circ }{\nabla }}_b\delta \mathit{\phi }.\hfill \end{array}$$

(26)

Proof. (25) is immediate. For (26), using $K_{ab}^{\left(1\right)}=d_{t,ab}-d_{ab}$ and $d_{t,ab}=-\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi }$,

$$\delta K_{ab}^{\left(1\right)}=-\delta \mathit{n}·{\partial }_a{\partial }_b\mathit{\phi }-\mathit{n}·{\partial }_a{\partial }_b\delta \mathit{\phi }.$$

Insert (24) and use $a_t^{cd}{\partial }_d\mathit{\phi }·{\partial }_a{\partial }_b\mathit{\phi }={\Gamma }_{ab}^c$, the Christoffel symbol of the deformed surface.    □

Remark 6 

(A common error). It is sometimes written that $\delta K_{ab}^{\left(1\right)}={\partial }_a\delta \mathit{\phi }·{\partial }_b\mathit{n}+{\partial }_a\mathit{\phi }·{\partial }_b\delta \mathit{n}-\delta ({\partial }_a\mathit{\phi }·{\partial }_b\mathit{n})$. The right-hand side is identically zero, since the third term is by the product rule the sum of the first two. The correct statement is (26). The presence of ${\Gamma }_{ab}^c$—the deformed, not reference, Christoffel symbol—is what makes $\delta {\mathit{K}}^{\left(1\right)}$ a tensor, and it is exactly the term that a naive derivation loses.

4. Stress as a Bundle-Valued Form; Virtual Work

4.1. The Piola Stress 2-Form

Let $\mathit{P}$ be the first Piola–Kirchhoff stress of the three-dimensional body, a two-point tensor. Its natural differential-geometric avatar is the ${\mathbb{R}}^3$-valued 2-form

$$\Xi =\frac{1}{2}{{\Xi }^I}_{JK}{\Theta }^J\wedge {\Theta }^K{\mathit{e}}_I,{{\Xi }^I}_{JK}={ϵ}_{JKL}{P}^{IL},$$

(27)

so that for any surface $\Sigma \subset \mathcal{S}$ with unit normal $\mathit{\nu }$ the traction resultant is ${\int }_{\Sigma }\Xi$. Equivalently $\Xi =★(\mathit{P}·\Theta )$, the Hodge dual taken in the reference metric (6).

Proposition 5 

(Cauchy’s theorem is Stokes’ theorem). Balance of linear momentum for the three-dimensional body,

$$Div\mathit{P}+{\rho }_0\mathit{b}={\rho }_0\ddot{\mathbf{\Phi }},$$

(28)

is equivalent to the statement

$${\mathrm{d}}^{\nabla }\Xi +\left({\rho }_0\mathit{b}-{\rho }_0\ddot{\mathbf{\Phi }}\right)\mathrm{d}V=0\mathrm{on}\mathcal{S},$$

(29)

and integrating (29) over any subbody $\mathcal{P}\subseteq \mathcal{S}$ and applying (14) returns Cauchy’s theorem ${\int }_{\partial \mathcal{P}}\Xi +{\int }_{\mathcal{P}}{\rho }_0\mathit{b}\mathrm{d}V={\int }_{\mathcal{P}}{\rho }_0\ddot{\mathbf{\Phi }}\mathrm{d}V$.

This is the shell-theoretic instance of the general geometric formulation of continuum mechanics in [33]. Balance of angular momentum is, in this language, the statement that the ${\mathbb{R}}^3$-valued 3-form $\mathbf{\Phi }\wedge {\mathrm{d}}^{\nabla }\Xi +\mathrm{d}\mathbf{\Phi }\wedge \Xi$ has vanishing skew part, which is equivalent to symmetry of the second Piola–Kirchhoff stress $\mathit{S}={\mathit{F}}^{-1}\mathit{P}$.

4.2. Three-dimensional virtual work

For admissible $\delta \mathbf{\Phi }$ vanishing where $\mathbf{\Phi }$ is prescribed, pairing (29) with $\delta \mathbf{\Phi }$ and applying (14) gives

$$\underset{\delta W_{\mathrm{int}}}{\underbrace{{\int }_{\mathcal{S}}\mathit{S}:\delta \mathit{E}\mathrm{d}V}}=\underset{\delta W_{\mathrm{ext}}}{\underbrace{{\int }_{\mathcal{S}}{\rho }_0\mathit{b}·\delta \mathbf{\Phi }\mathrm{d}V+{\int }_{{\partial }_t\mathcal{S}}\overline{\mathit{t}}·\delta \mathbf{\Phi }\mathrm{d}A}}-{\int }_{\mathcal{S}}{\rho }_0\ddot{\mathbf{\Phi }}·\delta \mathbf{\Phi }\mathrm{d}V.$$

(30)

Everything that follows is obtained from (30) by fibre integration; no further mechanical assumption is introduced.

5. Dimensional Reduction by Fibre Integration

5.1. Resultants

Using Proposition 3 with the truncation of Corollary 1, $\delta E_{ab}\left(z\right)=\delta E_{ab}^{\left(0\right)}+z\delta K_{ab}^{\left(1\right)}$, and the exact volume form (7),

$$\delta W_{\mathrm{int}}={\int }_{\mathcal{M}}{\int }_{-h/2}^{h/2}{S}^{ab}\left(\delta E_{ab}^{\left(0\right)}+z\delta K_{ab}^{\left(1\right)}\right){\mu }_0\left(z\right)\mathrm{d}z\mathrm{d}A={\int }_{\mathcal{M}}\left(N^{ab}\delta E_{ab}^{\left(0\right)}+M^{ab}\delta K_{ab}^{\left(1\right)}\right)\mathrm{d}A,$$

(31)

which defines the resultants as the weighted fibre moments

$$\begin{array}{|c|}\hline N^{ab}:={\int }_{-h/2}^{h/2}{S}^{ab}{\mu }_0\left(z\right)\mathrm{d}z,M^{ab}:={\int }_{-h/2}^{h/2}{S}^{ab}z{\mu }_0\left(z\right)\mathrm{d}z.\\ \hline\end{array}$$

(32)

The weight ${\mu }_0$ is the exact shifter determinant; setting ${\mu }_0\equiv 1$ is an additional approximation of relative order $h/R$ which we do not make.

5.2. Exact Through-Thickness Stiffness

With the linear-elastic law $S^{ab}=H^{abcd}(E_{cd}^{\left(0\right)}+z{K}_{cd}^{\left(1\right)})$ of Section 7 and ${\mu }_0=1+2Hz+K{z}^2$, the elementary integrals

$${\int }_{-h/2}^{h/2}{\mu }_0\mathrm{d}z=h\left(1+{\displaystyle \frac{K{h}^2}{12}}\right),{\int }_{-h/2}^{h/2}z{\mu }_0\mathrm{d}z={\displaystyle \frac{H{h}^3}{6}},{\int }_{-h/2}^{h/2}{z}^2{\mu }_0\mathrm{d}z={\displaystyle \frac{h^3}{12}}+{\displaystyle \frac{K{h}^5}{80}},$$

(33)

give the exact resultant constitutive law

$$\left(\begin{array}{c}{N}^{ab}\\ M^{ab}\end{array}\right)=\left(\begin{array}{cc}\mathsf{A}& \mathsf{B}\\ \mathsf{B}& \mathsf{D}\end{array}\right)H^{abcd}\left(\begin{array}{c}{E}_{cd}^{\left(0\right)}\\ K_{cd}^{\left(1\right)}\end{array}\right),\begin{array}{cc}\hfill \mathsf{A}& =h\left(1+\frac{1}{12}K{h}^2\right),\hfill \\ \hfill \mathsf{B}& =\frac{1}{6}H{h}^3,\hfill \\ \hfill \mathsf{D}& =\frac{1}{12}{h}^3\left(1+\frac{3}{20}K{h}^2\right).\hfill \end{array}$$

(34)

Remark 7 

(Curvature-induced coupling). Even for a homogeneous isotropic shell, initial curvature produces a membrane–bending coupling $\mathsf{B}=H{h}^3/6$, of relative size $\frac{1}{2}h/R$ compared with $\sqrt{\mathsf{A}\mathsf{D}}$. This is a purely geometric effect—the outer fibres of a curved shell are longer than the inner ones—and it is invisible in formulations that set ${\mu }_0\equiv 1$. For the caps of Section 10, $h/R=2\times {10}^{-3}$ and the effect is negligible, as it should be; for thick pressure-vessel heads or laminated tubes with $h/R\sim 0.1$ it is not.

5.3. Inertia

Fibre integration of the inertial term in (30), using $\ddot{\mathbf{\Phi }}=\ddot{\mathit{\phi }}+z\ddot{\mathit{n}}$, gives

$${\int }_{\mathcal{S}}{\rho }_0\ddot{\mathbf{\Phi }}·\delta \mathbf{\Phi }\mathrm{d}V={\int }_{\mathcal{M}}\left[{\rho }_0\mathsf{A}\ddot{\mathit{\phi }}·\delta \mathit{\phi }+{\rho }_0\mathsf{B}\left(\ddot{\mathit{n}}·\delta \mathit{\phi }+\ddot{\mathit{\phi }}·\delta \mathit{n}\right)+{\rho }_0\mathsf{D}\ddot{\mathit{n}}·\delta \mathit{n}\right]\mathrm{d}A.$$

(35)

The last term is the rotary inertia, of relative size $O(h^2/L^2)$ where L is the shortest deformation wavelength; it is retained here for completeness and neglected in Section 10, where $h/L\lesssim 2\times {10}^{-2}$.

5.4. Strong Form

Theorem 3 

(Local balance laws). Define the difference tensor ${\Delta }_{ab}^c:={\Gamma }_{ab}^c-{\tilde{\Gamma }}_{ab}^c$ between the deformed and reference Christoffel symbols (a genuine $(1,2)$-tensor), the effective membrane resultant and the effective shear

$$\begin{array}{|c|}\hline {\widehat{N}}^{ab}:=N^{ab}+M^{ac}{{\tilde{b}}_c}^b,Q^a:=-\left({\tilde{\nabla }}_b{M}^{ab}+M^{cb}{\Delta }_{cb}^a\right),{\tilde{b}}_{ab}=\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi },\\ \hline\end{array}$$

(36)

indices on $\tilde{b}$ being raised with $a_t^{ab}$, and the stress vector

$${\mathit{T}}^a:={\widehat{N}}^{ab}{\partial }_b\mathit{\phi }+Q^a\mathit{n}.$$

(37)

Then (30) is equivalent to the vector balance law

$$\begin{array}{|c|}\hline {\tilde{\nabla }}_a{\mathit{T}}^a+\overline{\mathit{p}}={\rho }_0\mathsf{A}\ddot{\mathit{\phi }}+{\rho }_0\mathsf{B}\ddot{\mathit{n}}\\ \hline\end{array}$$

(38)

on $\mathcal{M}$, whose tangential and normal projections read

$$\begin{array}{cc}\hfill {\tilde{\nabla }}_a{\widehat{N}}^{ab}-Q^a{{\tilde{b}}_a}^b+{\overline{p}}^b& ={\rho }_0\mathsf{A}\left(\ddot{\mathit{\phi }}·{\partial }^b\mathit{\phi }\right)+\cdots ,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\tilde{\nabla }}_a{Q}^a+{\widehat{N}}^{ab}{\tilde{b}}_{ab}+{\overline{p}}^n& ={\rho }_0\mathsf{A}\left(\mathit{n}·\ddot{\mathit{\phi }}\right)+\cdots \hfill \end{array}$$

(40)

Proof. 

Insert (25)–(26) into (31). Since $\tilde{\nabla }$ is compatible with $\mathrm{d}A$ but $\stackrel{\circ }{\nabla }$ is not, write ${\stackrel{\circ }{\nabla }}_a{\stackrel{\circ }{\nabla }}_b\delta \mathit{\phi }={\tilde{\nabla }}_a{\tilde{\nabla }}_b\delta \mathit{\phi }-{\Delta }_{ab}^c{\tilde{\nabla }}_c\delta \mathit{\phi }$, so that

$$\delta W_{\mathrm{int}}={\int }_{\mathcal{M}}\left[N^{ab}{\partial }_b\mathit{\phi }·{\tilde{\nabla }}_a\delta \mathit{\phi }-M^{ab}\mathit{n}·{\tilde{\nabla }}_a{\tilde{\nabla }}_b\delta \mathit{\phi }+M^{ab}{\Delta }_{ab}^c\mathit{n}·{\tilde{\nabla }}_c\delta \mathit{\phi }\right]\mathrm{d}A.$$

Integrating the second term by parts once and using the Weingarten relation ${\partial }_b\mathit{n}=-{{\tilde{b}}_b}^c{\partial }_c\mathit{\phi }$ collects the tangential parts into ${\widehat{N}}^{ab}$ and the normal parts into $Q^a$, yielding $\delta W_{\mathrm{int}}={\int }_{\mathcal{M}}{\mathit{T}}^a·{\tilde{\nabla }}_a\delta \mathit{\phi }\mathrm{d}A+{\oint }_{\partial \mathcal{M}}(\cdots )$. A second integration by parts and the fundamental lemma give (38). Projecting on ${\partial }_b\mathit{\phi }$ and $\mathit{n}$, and using $\mathit{n}·{\tilde{\nabla }}_a\left({\widehat{N}}^{ab}{\partial }_b\mathit{\phi }\right)={\widehat{N}}^{ab}{\tilde{b}}_{ab}$, gives (39)–(40). The omitted terms “…” are the $\mathsf{B}$- and $\mathsf{D}$-weighted director inertias of (35).    □

Remark 8 

(What is usually missing). Three features of Theorem 3 are commonly absent from presentations of Kirchhoff–Love shell theory, though all three are present in the careful treatments [28,34,41]:

(a) 

the replacement of $N^{ab}$ by ${\widehat{N}}^{ab}=N^{ab}+M^{ac}{{\tilde{b}}_c}^b$, which is what makes ${\mathit{T}}^a$ the true traction vector;

(b) 

the shear–curvature coupling $-Q^a{{\tilde{b}}_a}^b$ in the tangential equation (39). Writing simply ${\tilde{\nabla }}_b{N}^{ab}+{\overline{n}}^a={\rho }_{0}h{\ddot{\phi }}^a$ is correct only for a flat plate;

(c) 

the term $M^{cb}{\Delta }_{cb}^a$ in $Q^a$, which is $O(\parallel \mathit{E}\parallel )$ and hence negligible in small-strain/large-rotation regimes but must be present for the formulation to be exactly variationally consistent (and therefore for Newton’s method to converge quadratically).

The term ${\widehat{N}}^{ab}{\tilde{b}}_{ab}$ in (40) is the geometric-stiffness density: it is the mechanism by which membrane compression destabilises a curved shell, and it is what produces the limit points computed in Section 10.8.

6. Boundary Conditions and Corner Forces

The boundary term generated in the proof of Theorem 3 is

$${\oint }_{\partial \mathcal{M}}\left[\left({\mathit{T}}^a{\nu }_a\right)·\delta \mathit{\phi }-M^{ab}{\nu }_a\mathit{n}·{\tilde{\nabla }}_b\delta \mathit{\phi }\right]\mathrm{d}s,$$

(41)

with $\mathit{\nu }$ the outward unit conormal and $\mathit{\tau }$ the unit tangent of $\partial \mathcal{M}$. Because $\delta \mathit{\phi }$ and its tangential derivative along $\partial \mathcal{M}$ are not independent, (41) must be reduced further. Decomposing ${\tilde{\nabla }}_b={\nu }_b{\partial }_{\nu }+{\tau }_b{\partial }_{\tau }$ and writing $M_{\nu \nu }=M^{ab}{\nu }_a{\nu }_b$, $M_{\nu \tau }=M^{ab}{\nu }_a{\tau }_b$, integration by parts along $\partial \mathcal{M}$ gives

Proposition 6 

(Kirchhoff boundary conditions).

$$\left(41\right)={\oint }_{\partial \mathcal{M}}\left[\left({\mathit{T}}^a{\nu }_a+{\partial }_{\tau }\left(M_{\nu \tau }\right)\mathit{n}\right)·\delta \mathit{\phi }-M_{\nu \nu }\left(\mathit{n}·{\partial }_{\nu }\delta \mathit{\phi }\right)\right]\mathrm{d}s-\sum _{\mathit{corners}i}{\left[\left[M_{\nu \tau }\right]\right]}_i\mathit{n}·\delta \mathit{\phi }{|}_i,$$

(42)

so that the natural (traction) boundary conditions are

$$\underset{\mathit{effective}\mathit{traction}/\mathit{Kirchhoff}\mathit{shear}}{\underbrace{{\mathit{T}}^a{\nu }_a+{\partial }_{\tau }\left(M_{\nu \tau }\right)\mathit{n}=\overline{\mathit{t}}}},\underset{\mathit{bending}\mathit{moment}}{\underbrace{M_{\nu \nu }=\overline{m}}},\underset{\mathit{corner}\mathit{force}}{\underbrace{{\left[\left[M_{\nu \tau }\right]\right]}_i={\overline{f}}_i}},$$

(43)

and the essential conditions are prescription of $\mathit{\phi }$ and of the normal slope $\mathit{n}·{\partial }_{\nu }\mathit{\phi }$.

Remark 9. 

The corner forces ${\left[\left[M_{\nu \tau }\right]\right]}_i$ are not an artefact. A Kirchhoff–Love shell has three boundary conditions per edge, not four, and the twisting moment cannot be prescribed independently of the shear; the concentrated corner reactions are what restores consistency. Any formulation that lists four independent edge conditions is over-determined. The work-conjugate structure of shell boundary conditions at finite rotation is analysed in detail by Makowski and Pietraszkiewicz [31]. The clamped conditions used in Section 10 ($\mathit{\phi }$ and ${\partial }_{\nu }\mathit{\phi }$ prescribed) are essential throughout, so no corner forces arise there; they would for a simply-supported polygonal panel.

7. Constitutive Closure

7.1. Plane-Stress Condensation

Take a Saint-Venant–Kirchhoff material with strain energy per unit reference volume $W_3\left(\mathit{E}\right)=\mu \mathit{E}:\mathit{E}+\frac{\lambda }{2}{(tr\mathit{E})}^2$, $\lambda ,\mu$ the Lamé parameters. The Kirchhoff–Love kinematics forces $E_{a3}=0$ but leaves $E_{33}$ undetermined; it is fixed by the plane-stress condition $S^{33}=0$, which gives

$$E_{33}=-{\displaystyle \frac{\lambda }{\lambda +2\mu }}\left(E_{11}+E_{22}\right),$$

(44)

and substituting back replaces $\lambda$ by the reduced modulus ${\lambda }^{*}={\displaystyle \frac{2\lambda \mu }{\lambda +2\mu }}$.

Proposition 7 

(Reduced elasticity tensor).

$$\begin{array}{|c|}\hline H^{abcd}={\lambda }^{*}{\delta }^{ab}{\delta }^{cd}+\mu \left({\delta }^{ac}{\delta }^{bd}+{\delta }^{ad}{\delta }^{bc}\right)={\displaystyle \frac{E}{1-{\nu }^2}}\left[\nu {\delta }^{ab}{\delta }^{cd}+{\displaystyle \frac{1-\nu }{2}}\left({\delta }^{ac}{\delta }^{bd}+{\delta }^{ad}{\delta }^{bc}\right)\right],\\ \hline\end{array}$$

(45)

components referred to the orthonormal coframe $\{{\Theta }^1,{\Theta }^2\}$.

Remark 10 

(Dimensional consistency). $H^{abcd}$ already carries the modulus and has units of stress. The resultant law is therefore

$$N^{ab}=\mathsf{A}{H}^{abcd}{E}_{cd}^{\left(0\right)}+\mathsf{B}{H}^{abcd}{K}_{cd}^{\left(1\right)},M^{ab}=\mathsf{B}{H}^{abcd}{E}_{cd}^{\left(0\right)}+\mathsf{D}{H}^{abcd}{K}_{cd}^{\left(1\right)},$$

(46)

with $\mathsf{A},\mathsf{B},\mathsf{D}$ from (34). Writing instead $N^{ab}=C{H}^{abcd}{E}_{cd}^{\left(0\right)}$ with $C=Eh/(1-{\nu }^2)$—as is occasionally seen—multiplies the modulus in twice and produces membrane stiffnesses too large by the factor $E/(1-{\nu }^2)$.

In the thin-shell limit $h/R\to 0$, $\mathsf{A}\to h$, $\mathsf{B}\to 0$, $\mathsf{D}\to h^3/12$, and (46) collapses to the familiar uncoupled form with extensional and flexural rigidities

$$C={\displaystyle \frac{Eh}{1-{\nu }^2}},B_f={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}.$$

(47)

7.2. Comparison with the Higher-Order Model of Simo, Fox and Rifai

Simo, Fox and Rifai ([42], Table 2) give a higher-order resultant constitutive model obtained by asymptotic expansion of the three-dimensional stored energy, in which the membrane and bending resultants read

$$\begin{array}{cc}\hfill {\tilde{n}}^{ab}& ={\displaystyle \frac{E}{1-{\nu }^2}}\left(h{H}_0^{abcd}{\overline{\epsilon }}_{cd}-\frac{h^3}{6}{H}_1^{abcd}{\overline{\rho }}_{cd}+\frac{h^3}{3}{H}_2^{abcd}{\overline{\epsilon }}_{cd}+\frac{h^3}{4}{H}_3^{abcd}{\overline{\epsilon }}_{cd}\right),\hfill \\ \hfill {\tilde{m}}^{\left(ab\right)}& ={\displaystyle \frac{E}{1-{\nu }^2}}\left(\frac{h^3}{12}{H}_0^{abcd}{\overline{\rho }}_{cd}-\frac{h^3}{6}{H}_1^{abcd}{\overline{\epsilon }}_{cd}\right),\hfill \end{array}$$

(48)

with $H_1$ linear in the second fundamental form $b_{ab}^0={\mathit{N}}_{,a}·{\mathit{\phi }}_{,b}^0$, $H_2$ quadratic in $b^0$, and $H_3$ linear in the third fundamental form $d_{ab}^0={\mathit{N}}_{,a}·{\mathit{N}}_{,b}$—which is exactly the tensor $\mathit{c}$ of Proposition 2. Three features of (48) should be recorded before any comparison is attempted.

(i)

The conventions coincide. Their second fundamental form $b_{ab}^0={\mathit{N}}_{,a}·{\mathit{\phi }}_{,b}^0$ is our $d_{ab}$ of (4), and their bending measure ${\overline{\rho }}_{ab}={\mathit{\phi }}_{,a}·{\mathit{t}}_{,b}-{\mathit{\phi }}_{,a}^0·{\mathit{t}}_{,b}^0$ ([42], Equation (2.11c)) is our ${\mathit{K}}^{\left(1\right)}$ (Definition 2) verbatim, the director $\mathit{t}$ replacing the normal $\mathit{n}$. No sign reconciliation is needed, and the discrepancies found below are therefore real and not artefacts of Remark 2.

(ii)

The model is hyperelastic. The cross-coefficient $-\frac{h^3}{6}{H}_1$ is shared by $\tilde{n}$ and $\tilde{m}$, and $H_1$ possesses the major symmetry $H_1^{abcd}=H_1^{cdab}$. Hence (48) derives from a quadratic stored energy of exactly the form (51), with $\mathsf{A}{H}_0\mapsto h{H}_0+\frac{h^3}{3}{H}_2+\frac{h^3}{4}{H}_3$, $\mathsf{B}{H}_0\mapsto -\frac{h^3}{6}{H}_1$ and $\mathsf{D}{H}_0\mapsto \frac{h^3}{12}{H}_0$. The two models are thus structurally identical; only the coefficients differ.

(iii)

The resultants are the effective ones. The quantities ${\tilde{n}}^{ab}$, ${\tilde{m}}^{ab}$, ${\tilde{q}}^a$ of Table 2 are the effective resultants defined by ${\tilde{n}}^{ba}=n^{ba}-{\lambda }_{\mu }^b{\tilde{m}}^{a\mu }$, ${\tilde{q}}^a=q^a-{\lambda }_{\mu }^3{\tilde{m}}^{a\mu }$ ([42], Equation (2.5)), not the bare resultants $N^{ab}$ of (34). Comparing the two requires converting one to the other, and at order $h^3/R$ the conversion is not negligible.

Specialising both models to a sphere of radius R, where $a_{ab}^0={\delta }_{ab}$, $b_{ab}^0={\delta }_{ab}/R$, $d_{ab}^0={\delta }_{ab}/R^2$, one finds $H_1=(2/R)H_0$, $H_2=(1/R^2)H_0$, $H_3=(2/R^2)H_0$, and ${\lambda }_a^b={\delta }_a^b/R$, so that the coefficient of $H_0$ multiplying $\overline{\epsilon }$ and $\overline{\rho }$ in the membrane resultant is

$$\begin{array}{ccc}\mathrm{membrane}\mathrm{resultant}\hfill & \mathrm{coefficient}\mathrm{of}\overline{\epsilon }& \mathrm{coefficient}\mathrm{of}\overline{\rho }\\ \text{([42], Table 2), effective} \tilde{n}\hfill & h\left(1+\frac{5}{6}{h}^2/R^2\right)& -h^3/3R\\ \text{([42], Table 2), bare} n=\tilde{n}+R^{-1}\tilde{m}\hfill & h\left(1+\frac{1}{2}{h}^2/R^2\right)& -h^3/4R\\ \text{present, (34)}\hfill & h\left(1+\frac{1}{12}{h}^2/R^2\right)& +h^3/6R\end{array}$$

(49)

The third row is $\mathsf{A}$ and $\mathsf{B}$ of (34) with $H=1/R$, $K=1/R^2$. All entries, together with the three tensor identities above, are verified symbolically in code/verify_table2.py.

Remark 11 

(An inconsistent truncation, and its repair). Even between the two directly comparable rows of (49) the coupling coefficients differ by a factor of $-\frac{3}{2}$ and the $O(h^2/R^2)$ membrane correction by a factor of 6; against the effective resultant as tabulated the factors are $-2$ and 10. Converting between effective and bare resultants therefore narrows the gap but does not close it, and by item (i) above the sign reversal cannot be attributed to conventions. The discrepancy is nevertheless not an error in either derivation; it is the truncation warned about in Corollary 1. Equation (34) retains the shifter ${\mu }_0$ exactly but discards the $O\left(z^2\right)$ strain remainder $\frac{z^2}{2}({\mathit{c}}_t-\mathit{c})$, and both are of order $h^2/R^2$. Retaining one without the other is inconsistent. Carrying both through the fibre integral gives, in place of (34),

$$N^{ab}=H^{abcd}\left[\mathsf{A}{E}_{cd}^{\left(0\right)}+\mathsf{B}{K}_{cd}^{\left(1\right)}+\frac{h^3}{24}{C}_{cd}\right],M^{ab}=H^{abcd}\left[\mathsf{B}{E}_{cd}^{\left(0\right)}+\mathsf{D}{K}_{cd}^{\left(1\right)}+O\left(h^5\right)\right],$$

(50)

with $\mathit{C}:={\mathit{c}}_t-\mathit{c}$, and the $H_3$-type term of (48) is recovered. That the two models should differ at this order is not a surprise on their side either: Simo and Fox derive the lowest-order relations by asymptotic expansion of the three-dimensional Saint-Venant–Kirchhoff model retaining terms through $O(h/R)$, and note explicitly that products of $a^{\alpha \beta }$ with the curvature parameters ${\lambda }_{\beta }^{\alpha }$ are of higher order and are conventionally discarded, so that the resultant law is left depending on the reference surface through its first fundamental form alone ([41], Rem. 5.5(2)); the general dependence of resultant constitutive laws on reference geometry is the subject of Carroll and Naghdi [8]. Equations (48) and (34) are two different partial retentions of the terms that expansion discards, which is why they disagree. The practical consequence is small: ([42], Remarks 2.2) report that the lowest-order and higher-order models give numerically indistinguishable results, and that through-thickness integration of three-dimensional laws “leads to no improved accuracy for most thin, or moderately thick, shell problems”. For the shells of Section 10, $h/R=2\times {10}^{-3}$ and every term in (49) beyond the leading h is below ${10}^{-6}$ relative, so the factors of $-\frac{3}{2}$ and 6 are without practical consequence here. The point is one of principle: (34) should be read as the ${\mu }_0$-exact model, not as a consistent $O(h^2/R^2)$ theory, and (50) is the consistent one.

Remark 12 

(The third fundamental form as the carrier of the $O(h^2/R^2)$ corrections). The same tensor $\mathit{c}={\mathit{d}}^0$ appears in the third of the Table 2 relations, the transverse shear law ${\tilde{q}}^a=\kappa G\left(h{a}^{0ab}+\frac{h^3}{4}{d}^{0ab}\right){\overline{\gamma }}_b$ ([42], Table 2). That the leading curvature correction to every resultant stiffness—membrane, coupling and shear—is carried by ${\mathit{d}}^0$ and by the shifter built from it is not a coincidence: by Proposition 2 the pull-back metric is exactly $\mathit{a}+2z\mathit{d}+z^2\mathit{c}$, so $\mathit{c}$ is the only new tensor that $O\left(z^2\right)$ accuracy can produce. The present theory is shear-rigid and has no counterpart of $\tilde{q}$, but the structural agreement is a useful check on (50).

7.3. Stored Energy and Hyperelasticity

The reduced two-dimensional stored-energy function is

$$\mathsf{\Psi }\left({\mathit{E}}^{\left(0\right)},{\mathit{K}}^{\left(1\right)}\right)=\frac{1}{2}\mathsf{A}{H}^{abcd}{E}_{ab}^{\left(0\right)}{E}_{cd}^{\left(0\right)}+\mathsf{B}{H}^{abcd}{E}_{ab}^{\left(0\right)}{K}_{cd}^{\left(1\right)}+\frac{1}{2}\mathsf{D}{H}^{abcd}{K}_{ab}^{\left(1\right)}{K}_{cd}^{\left(1\right)},$$

(51)

with $N^{ab}=\partial \mathsf{\Psi }/\partial E_{ab}^{\left(0\right)}$, $M^{ab}=\partial \mathsf{\Psi }/\partial K_{ab}^{\left(1\right)}$. Because ${\mathit{E}}^{\left(0\right)}$ and ${\mathit{K}}^{\left(1\right)}$ are exactly objective (Theorem 2), $\mathsf{\Psi }$ is frame-indifferent by construction, and the resulting boundary-value problem is variational: for conservative loading, equilibria are the critical points of

$$\mathsf{\Pi }\left[\mathit{\phi }\right]={\int }_{\mathcal{M}}\mathsf{\Psi }\mathrm{d}A-{\mathsf{\Pi }}_{\mathrm{ext}}\left[\mathit{\phi }\right].$$

(52)

This is what licenses the energy-based numerical treatment of Section 10, including the use of the second variation as a stability criterion.

8. Reduction to Classical Theories

A formulation of this generality is only trustworthy if the standard theories drop out of it. We record the four relevant limits; each is a statement about which terms are retained in ${\mathit{E}}^{\left(0\right)}$, ${\mathit{K}}^{\left(1\right)}$ and Theorem 3.

(i)

Koiter’s linear shell theory.

Write $\mathit{\phi }={\mathit{\phi }}_0+\mathit{u}$ and linearise (19) at $\mathit{u}=0$:

$$E_{ab}^{\left(0\right)}\to {\epsilon }_{ab}=sym({\partial }_a\mathit{u}·{\partial }_b{\mathit{\phi }}_0),K_{ab}^{\left(1\right)}\to {\rho }_{ab}=-\mathit{N}·\left({\partial }_a{\partial }_b\mathit{u}-{\tilde{\Gamma }}_{ab}^c{\partial }_c\mathit{u}\right).$$

(53)

With $\mathsf{B}\to 0$ this is exactly Koiter’s energy [9,26]. The price is (23): the theory is no longer rotation-invariant.

(ii)

von Kármán plate.

Take $\mathcal{M}$ flat, $\mathit{u}=(u_1,u_2,w)$, retain quadratic terms only in w:

$$E_{ab}^{\left(0\right)}=\frac{1}{2}({\partial }_a{u}_b+{\partial }_b{u}_a)+\frac{1}{2}{\partial }_{a}w{\partial }_{b}w,K_{ab}^{\left(1\right)}=-{\partial }_a{\partial }_{b}w,$$

(54)

and (40) becomes $B_f{\Delta }^{2}w-N^{ab}{\partial }_a{\partial }_{b}w={\overline{p}}^n$, the von Kármán equation [52], with $\widehat{N}\to N$ because $\tilde{b}=O\left({\nabla }^{2}w\right)$ is already second order.

(iii)

Föppl membrane.

Set $\mathsf{D}=\mathsf{B}=0$. Then $Q^a\equiv 0$, (40) reduces to $N^{ab}{\tilde{b}}_{ab}+{\overline{p}}^n=0$ — the Young–Laplace law—and (39) to ${\tilde{\nabla }}_a{N}^{ab}+{\overline{p}}^b=0$.

(iv)

Inextensional bending.

Impose ${\mathit{E}}^{\left(0\right)}\equiv 0$. By Gauss’ theorema egregium (Theorem 1) [13] the deformation preserves K, so the admissible motions are the isometries of $\mathcal{M}$; the energy reduces to $\frac{1}{2}\mathsf{D}{H}^{abcd}{K}_{ab}^{\left(1\right)}{K}_{cd}^{\left(1\right)}$. This limit is the correct description whenever $B_f/\left(C{L}^2\right)\ll 1$ and it is the reason that developable shells roll up at essentially zero membrane stress — behaviour that a formulation with spurious rotational strain cannot reproduce.

9. Linearisation and the Consistent Tangent

Newton solution of (52) requires the second variation. Writing $G[\mathit{\phi };\delta \mathit{\phi }]=\delta \mathsf{\Pi }$, the tangent operator splits in the usual way into material and geometric parts:

$$\begin{array}{cc}\hfill \mathrm{D}G·\Delta \mathit{\phi }=& \underset{\mathrm{material}}{\underbrace{{\int }_{\mathcal{M}}\left(\begin{array}{c}\delta {\mathit{E}}^{\left(0\right)}:{\mathbb{C}}_{NN}:\Delta {\mathit{E}}^{\left(0\right)}+\delta {\mathit{E}}^{\left(0\right)}:{\mathbb{C}}_{NM}:\Delta {\mathit{K}}^{\left(1\right)}\\ +\delta {\mathit{K}}^{\left(1\right)}:{\mathbb{C}}_{MN}:\Delta {\mathit{E}}^{\left(0\right)}+\delta {\mathit{K}}^{\left(1\right)}:{\mathbb{C}}_{MM}:\Delta {\mathit{K}}^{\left(1\right)}\end{array}\right)\mathrm{d}A}}\hfill \\ \hfill & +\underset{\mathrm{geometric}}{\underbrace{{\int }_{\mathcal{M}}\left(N^{ab}\Delta \delta E_{ab}^{\left(0\right)}+M^{ab}\Delta \delta K_{ab}^{\left(1\right)}\right)\mathrm{d}A}},\hfill \end{array}$$

(55)

with ${\mathbb{C}}_{NN}=\mathsf{A}H$, ${\mathbb{C}}_{NM}={\mathbb{C}}_{MN}=\mathsf{B}H$, ${\mathbb{C}}_{MM}=\mathsf{D}H$ and

$$\Delta \delta E_{ab}^{\left(0\right)}=sym\left({\partial }_a\delta \mathit{\phi }·{\partial }_b\Delta \mathit{\phi }\right),\Delta \delta K_{ab}^{\left(1\right)}=-\Delta \mathit{n}·{\stackrel{\circ }{\nabla }}_a{\stackrel{\circ }{\nabla }}_b\delta \mathit{\phi }-\mathit{n}·\Delta \left({\stackrel{\circ }{\nabla }}_a{\stackrel{\circ }{\nabla }}_b\right)\delta \mathit{\phi }.$$

(56)

The geometric term is what carries the destabilising effect of membrane compression; dropping it (as in a “modified Newton” scheme) does not change the equilibrium path but destroys quadratic convergence and invalidates eigenvalue-based stability assessment. Both terms in (55) are symmetric, here because the loading is dead and the material hyperelastic, so that (55) is a second variation of $\mathsf{\Pi }$. It is worth recording that symmetry survives in settings where that argument is unavailable: in the one-director Cosserat model the restriction that balance of angular momentum places on the admissible constitutive equations is by itself enough to make the consistent tangent symmetric even away from equilibrium ([41], §1.2(vii), Equation (6.34)). In Section 10 we obtain (55) to machine accuracy by complex-step differentiation of the discrete energy, which sidesteps the algebra entirely while remaining exact.

Remark 13 

(Continuum versus discrete linearisation). Equation (55) is the linearisation of the continuum weak form, and it is not in general the tangent of the discrete problem. If the unknowns are interpolated nonlinearly—as they unavoidably are whenever a director or a rotation is interpolated on a manifold — linearisation and interpolation do not commute, so substituting the interpolation into (55) yields an operator that is only consistent to leading order ([43], Rem. 4.4); the correct object is the derivative of the discrete residual [44]. For the present discretisation the distinction is vacuous, since the unknown is $\mathit{\phi }$ itself and enters its nodal parameters linearly, so the two operations commute. It is nevertheless the reason the results below are obtained by differentiating the discrete energy rather than by discretising (55). Nor is the point cosmetic: an inexact tangent destroys the asymptotic quadratic rate (Section 10.11) and, more seriously, corrupts the detection of limit and bifurcation points, which are defined by singularity of the tangent itself ([43], Rem. 5.2).

A related trap is worth naming, because it has the same source. In the second variation it matters which object is held fixed. Here $\delta \mathit{n}$ is not independent—it is determined by $\delta \mathit{\phi }$ through the Kirchhoff–Love constraint (Proposition 4)—so $\Delta \delta {\mathit{K}}^{\left(1\right)}$ acquires the second term in (55). The Cosserat counterpart is that the primitive variation is the rotation $\delta \mathit{\theta }$ and not the director increment $\delta \mathit{t}=\delta \mathit{\theta }\times \mathit{t}$, so that $\Delta \left(\delta \mathit{\theta }\right)=\mathit{0}$ while $\Delta \left(\delta \mathit{t}\right)=-(\delta \mathit{t}·\Delta \mathit{t})\mathit{t}\ne \mathit{0}$ ([41], Rem. 6.3, [43], §2.2.1). In both cases treating the dependent quantity as free loses a term of the geometric tangent and with it the quadratic rate.

10. Numerical Case Study

This section does three things: it verifies the claims of Section 3 directly, it verifies a discretisation of the theory against classical closed-form results, and it applies the whole apparatus to a genuinely nonlinear problem—the snap-through of clamped spherical caps—including the dynamic signature of the instability.

Throughout, the material is aluminium: $E=70$GPa, $\nu =0.3$, ${\rho }_0=2700$kg m⁻³.

10.1. Verification of Exact Objectivity

Theorem 2 is tested on a generic doubly curved Monge patch

$${\mathit{\phi }}_0(u,v)=\left(u,v,f(u,v)\right),f=0.60u^2-0.35v^2+0.25uv+0.20u^3-0.15u^{2}v+0.10u{v}^2+0.08v^3,$$

(57)

which has both elliptic and hyperbolic regions, non-orthogonal coordinate lines and non-constant curvature. All surface derivatives are computed analytically, so the only error present is floating-point round-off. The patch is subjected to $\mathit{\phi }=\mathit{Q}(\mathit{k},\theta ){\mathit{\phi }}_0+\mathit{c}$ with a generic axis $\mathit{k}\propto (0.4,-0.7,0.6)$ and translation $\mathit{c}=(1.2,-0.8,2.5)$, and ${\mathit{E}}^{\left(0\right)}$, ${\mathit{K}}^{\left(1\right)}$ are evaluated at $(u,v)=(0.30,-0.20)$.

Figure 2 shows the outcome. The exact measures sit at $6\times {10}^{-18}$–$7\times {10}^{-17}$ relative, i.e., at the round-off floor, uniformly in $\theta$—including at $\theta ={90}^{\circ }$ and ${180}^{\circ }$, where no expansion-based theory retains any accuracy. The linearised measures reproduce (23) to graphical precision. Selected values are given in Table 1.

Figure 2. Spurious strain generated by a pure rigid rotation of the doubly curved patch (57). (a) The Green–Lagrange measures ${\mathit{E}}^{\left(0\right)},{\mathit{K}}^{\left(1\right)}$ of Definition 2 remain at round-off ($\sim {10}^{-17}$) over the full range $0.{05}^{\circ }$–${180}^{\circ }$, while the corresponding linearised measures follow $1-cos\theta$ exactly, as predicted by (23). (b) The same data expressed as fictitious stress $E\parallel \mathit{\epsilon }\parallel$; the yield stress of Al 7075-T6 is exceeded by the numerical artefact alone at about $4.7^{\circ }$ of rigid rotation.

Figure 2. Spurious strain generated by a pure rigid rotation of the doubly curved patch (57). (a) The Green–Lagrange measures ${\mathit{E}}^{\left(0\right)},{\mathit{K}}^{\left(1\right)}$ of Definition 2 remain at round-off ($\sim {10}^{-17}$) over the full range $0.{05}^{\circ }$–${180}^{\circ }$, while the corresponding linearised measures follow $1-cos\theta$ exactly, as predicted by (23). (b) The same data expressed as fictitious stress $E\parallel \mathit{\epsilon }\parallel$; the yield stress of Al 7075-T6 is exceeded by the numerical artefact alone at about $4.7^{\circ }$ of rigid rotation.

Table 1. Relative spurious strain under rigid rotation of the patch (57). “Present” refers to ${\mathit{E}}^{\left(0\right)}$ of Definition 2; “linear” to ${\epsilon }^{\mathrm{lin}}$ of Remark 5. The last column is the equivalent fictitious stress $E\parallel {\epsilon }^{\mathrm{lin}}\parallel$.

$\mathit{\theta }$ [deg]	present $\parallel {\mathit{E}}^{\left(0\right)}\parallel$	present $\parallel {\mathit{K}}^{\left(1\right)}\parallel$	linear $\parallel {\mathit{\epsilon }}^{\mathbf{lin}}\parallel$	fictitious stress [MPa]
0.5	$7.3\times {10}^{-17}$	$2.8\times {10}^{-16}$	$2.96\times {10}^{-5}$	2.1
1.0	$7.3\times {10}^{-17}$	$1.5\times {10}^{-16}$	$1.18\times {10}^{-4}$	8.3
2.0	$1.3\times {10}^{-17}$	$1.5\times {10}^{-16}$	$4.73\times {10}^{-4}$	33
5.0	$6.3\times {10}^{-18}$	$1.4\times {10}^{-16}$	$2.95\times {10}^{-3}$	207
10.0	$6.3\times {10}^{-18}$	$6.8\times {10}^{-17}$	$1.18\times {10}^{-2}$	825
20.0	$7.1\times {10}^{-17}$	$1.4\times {10}^{-16}$	$4.68\times {10}^{-2}$	3277
45.0	$7.3\times {10}^{-17}$	$1.4\times {10}^{-16}$	$2.27\times {10}^{-1}$	15918
90.0	$7.4\times {10}^{-17}$	$1.2\times {10}^{-17}$	$7.76\times {10}^{-1}$	54341

Table 1. Relative spurious strain under rigid rotation of the patch (57). “Present” refers to ${\mathit{E}}^{\left(0\right)}$ of Definition 2; “linear” to ${\epsilon }^{\mathrm{lin}}$ of Remark 5. The last column is the equivalent fictitious stress $E\parallel {\epsilon }^{\mathrm{lin}}\parallel$.

$\mathit{\theta }$ [deg]	present $\parallel {\mathit{E}}^{\left(0\right)}\parallel$	present $\parallel {\mathit{K}}^{\left(1\right)}\parallel$	linear $\parallel {\mathit{\epsilon }}^{\mathbf{lin}}\parallel$	fictitious stress [MPa]
0.5	$7.3\times {10}^{-17}$	$2.8\times {10}^{-16}$	$2.96\times {10}^{-5}$	2.1
1.0	$7.3\times {10}^{-17}$	$1.5\times {10}^{-16}$	$1.18\times {10}^{-4}$	8.3
2.0	$1.3\times {10}^{-17}$	$1.5\times {10}^{-16}$	$4.73\times {10}^{-4}$	33
5.0	$6.3\times {10}^{-18}$	$1.4\times {10}^{-16}$	$2.95\times {10}^{-3}$	207
10.0	$6.3\times {10}^{-18}$	$6.8\times {10}^{-17}$	$1.18\times {10}^{-2}$	825
20.0	$7.1\times {10}^{-17}$	$1.4\times {10}^{-16}$	$4.68\times {10}^{-2}$	3277
45.0	$7.3\times {10}^{-17}$	$1.4\times {10}^{-16}$	$2.27\times {10}^{-1}$	15918
90.0	$7.4\times {10}^{-17}$	$1.2\times {10}^{-17}$	$7.76\times {10}^{-1}$	54341

The practical reading of Table 1 is the fourth row. A shell element rotating rigidly through $5^{\circ }$—a routine occurrence in a buckling or crash calculation—acquires, in a linearised formulation, a spurious strain of $3\times {10}^{-3}$. Physical membrane strains in a buckling shell are of exactly this order. The artefact is therefore not a small correction to the answer; it is the same size as the answer.

10.2. Verification of the Truncation Estimate

Corollary 1 claims that discarding ${\mathit{c}}_t-\mathit{c}$ costs $O\left(z^2\right)$. This is tested on the same patch under a genuine (non-rigid) deformation $\mathit{\phi }=\mathit{Q}\mathit{F}{\mathit{\phi }}_0$ with $\mathit{F}$ a non-symmetric stretch of $15\%$ and $\mathit{Q}$ a ${63}^{\circ }$ rotation, by comparing the exact $\mathit{E}\left(z\right)$ from (20) against ${\mathit{E}}^{\left(0\right)}+z{\mathit{K}}^{\left(1\right)}$. Figure 3 confirms the quadratic rate over three decades.

Figure 3. Thin-shell truncation error $\parallel \mathit{E}\left(z\right)-({\mathit{E}}^{\left(0\right)}+z{\mathit{K}}^{\left(1\right)})\parallel$ against fibre coordinate z, confirming the $O\left(z^2\right)$ estimate (21).

Figure 3. Thin-shell truncation error $\parallel \mathit{E}\left(z\right)-({\mathit{E}}^{\left(0\right)}+z{\mathit{K}}^{\left(1\right)})\parallel$ against fibre coordinate z, confirming the $O\left(z^2\right)$ estimate (21).

10.3. Axisymmetric Reduction of the General Theory

For the case study we specialise to shells of revolution undergoing axisymmetric, torsionless deformation. This retains both principal curvatures, non-constant Lamé parameters and non-zero Gaussian curvature, so it exercises the full geometric content of Section 2–Section 5 while reducing the field equations to one spatial dimension.

Let s be reference meridian arclength, $\theta$ the circumferential angle. The reference meridian is $(r_0\left(s\right),z_0\left(s\right))$ with tangent angle ${\varphi }_0$,

$$r_0^{\prime }=cos{\varphi }_0,z_0^{\prime }=sin{\varphi }_0,A_1=1,A_2=r_0,{\displaystyle \frac{1}{R_1}}={\varphi }_0^{\prime },{\displaystyle \frac{1}{R_2}}={\displaystyle \frac{sin{\varphi }_0}{r_0}}.$$

(58)

One checks immediately that (10) is trivial and (11) reduces to ${(sin{\varphi }_0)}^{\prime }=r_0^{\prime }{\varphi }_0^{\prime }$, an identity—an instance of Theorem 1.

The deformed meridian is $\left(r\right(s),z(s\left)\right)$; the Kirchhoff–Love constraint means the tangent angle $\varphi$ is not independent:

$${\lambda }_s=\sqrt{r^{\prime 2}+z^{\prime 2}},cos\varphi ={\displaystyle \frac{r^{\prime }}{{\lambda }_s}},sin\varphi ={\displaystyle \frac{z^{\prime }}{{\lambda }_s}},{\lambda }_{\theta }={\displaystyle \frac{r}{r_0}}.$$

(59)

Evaluating Definition 2 in the orthonormal reference coframe $\{\mathrm{d}s,r_0\mathrm{d}\theta \}$ gives

Proposition 8 

(Axisymmetric strain measures).

$$\begin{array}{|c|}\hline \begin{array}{cccc}\hfill E_s& =\frac{1}{2}\left({\lambda }_s^2-1\right),\hfill & \hfill E_{\theta }& =\frac{1}{2}\left({\lambda }_{\theta }^2-1\right),\hfill \\ \hfill K_s& ={\lambda }_s{\varphi }^{\prime }-{\varphi }_0^{\prime }={\displaystyle \frac{r^{\prime }{z}^{\prime \prime }-z^{\prime }{r}^{\prime \prime }}{{\lambda }_s}}-{\varphi }_0^{\prime },\hfill & \hfill K_{\theta }& ={\displaystyle \frac{1}{r_0}}\left({\displaystyle \frac{r}{r_0}}{\displaystyle \frac{z^{\prime }}{{\lambda }_s}}-sin{\varphi }_0\right).\hfill \end{array}\\ \hline\end{array}$$

(60)

Remark 14 

(Green–Lagrange versus naive curvature change). Note that $K_s={\lambda }_s^2{\kappa }_s-{\kappa }_s^0$ and $K_{\theta }={\lambda }_{\theta }^2{\kappa }_{\theta }-{\kappa }_{\theta }^0$, where ${\kappa }_s={\varphi }^{\prime }/{\lambda }_s$ and ${\kappa }_{\theta }=sin\varphi /r$ are the true current principal curvatures. The Green–Lagrange bending measure is therefore not $\kappa -{\kappa }^0$ but a stretch-weighted version of it; the two differ at relative order $\parallel {\mathit{E}}^{\left(0\right)}\parallel$. In the small-strain/large-rotation regime typical of thin shells the difference is immaterial, but it is not zero, and using $\kappa -{\kappa }^0$ with the energy (51) is variationally inconsistent. The correction is of the same type as the one distinguishing the Budiansky–Sanders “best” first-order theory [6,38] from Love’s, namely a term bilinear in the membrane strain and the second fundamental form; see ([42], Remarks 2.2), who report that the choice makes no significant numerical difference. As a concrete illustration, a sphere of radius $R_0$ inflated uniformly by a stretch λ has $K_s=K_{\theta }=(\lambda -1)/R_0\ne 0$ even though it remains a sphere.

The stored energy (51) in the thin-shell limit (47) becomes

$${\mathsf{\Pi }}_{\mathrm{int}}={\int }_0^L\left[\frac{C}{2}\left(E_s^2+2\nu E_s{E}_{\theta }+E_{\theta }^2\right)+\frac{B_f}{2}\left(K_s^2+2\nu K_s{K}_{\theta }+K_{\theta }^2\right)\right]2\pi r_0\mathrm{d}s.$$

(61)

Loading is a uniform live (follower) pressure, treated exactly through its potential

$${\mathsf{\Pi }}_{\mathrm{ext}}=pV\left(\mathit{\phi }\right),V=\pi {\int }_0^L{r}^2{z}^{\prime }\mathrm{d}s,$$

(62)

V being the enclosed volume. Because (62) is an exact potential the problem remains variational at finite deformation, and the tangent operator stays symmetric—a property that pressure loads treated as configuration-independent surface tractions do not possess.

10.4. Discretisation

Since (61) involves second derivatives of $(r,z)$, the trial space must be $C^1$ (Remark 4). We use cubic Hermite elements with nodal degrees of freedom $(r,r^{\prime },z,z^{\prime })$ and 5-point Gauss–Legendre quadrature. Three implementation choices are worth recording.

(1)

Discretely strain-free reference. The reference meridian is interpolated by the same Hermite basis and (60) is evaluated symmetrically in the current and reference fields. The discrete reference state then has exactly zero energy and zero residual (verified: $0.0$ and $0.0$ in double precision), rather than the $O\left(h_e^4\right)$ spurious initial stress produced by comparing an interpolated current state against an exact reference.

(2)

Complex-step gradients. Every expression in (60) is complex-analytic—in particular $\mathrm{atan}2$ is avoided by carrying $sin\varphi ,cos\varphi$ directly —— the complex-step derivative [32,46] applies, and the residual is obtained as $\partial \mathsf{\Pi }/\partial u_i=\Im \mathsf{\Pi }(u+\mathrm{i}\epsilon {\mathit{e}}_i)/\epsilon$ with $\epsilon ={10}^{-30}$, exact to machine precision and free of subtractive cancellation. The tangent (55) follows by central differencing of these gradients; it agrees with a finite-difference tangent of the residual to $7.6\times {10}^{-8}$ relative and is exactly symmetric, confirming the variational consistency of Theorem 3.

(3)

Path following. Equilibria are traced under apex-deflection control, with the pressure p carried as an extra unknown and the equilibrium equation at the controlled degree of freedom appended, so that limit points in p are passed without difficulty. This is a displacement-controlled variant of the arclength schemes of Wempner [54], Riks [37] and Crisfield [12].

10.5. Verification Against Classical Plate Solutions

The flat-plate specialisation (${\varphi }_0\equiv 0$, $r_0=s$) of (60) is tested against two classical results for a clamped circular plate with immovable edge, radius $a=100$mm, thickness $h=1$mm.

In the small-deflection limit the computed centre deflection matches $w_0=p{a}^4/\left(64B_f\right)$ to $1\times {10}^{-4}$ relative at $\overline{w}=w_0/h=1.2\times {10}^{-3}$, drifting to $0.06\%$ at $\overline{w}=0.034$ and $0.8\%$ at $\overline{w}=0.12$—the departure being physical membrane stiffening, not error. In the large-deflection regime, Table 2 compares against the classical one-term Galerkin result of Timoshenko and Woinowsky-Krieger ([52], Art. 97) $P^{*}\equiv p{a}^4/\left(E{h}^4\right)=5.86\overline{w}+2.86{\overline{w}}^3$.

Figure 4. Clamped circular plate under uniform pressure. (a) Large-deflection response against Timoshenko’s $P^{*}=5.86\overline{w}+2.86{\overline{w}}^3$ and the linear solution. (b) Departure of the computed centre deflection from the small-deflection result $w_0=p{a}^4/64B_f$, showing the onset of membrane stiffening.

Figure 4. Clamped circular plate under uniform pressure. (a) Large-deflection response against Timoshenko’s $P^{*}=5.86\overline{w}+2.86{\overline{w}}^3$ and the linear solution. (b) Departure of the computed centre deflection from the small-deflection result $w_0=p{a}^4/64B_f$, showing the onset of membrane stiffening.

Table 2. Clamped circular plate, large deflection. $P^{*}=p{a}^4/\left(E{h}^4\right)$. The growing deviation is expected: the reference expression is a one-term Galerkin approximation whose cubic coefficient is known to be low at large $\overline{w}$, whereas the present computation is a converged solution of the geometrically exact problem.

p [Pa]	$\overline{\mathit{w}}={\mathit{w}}_0/\mathit{h}$	${\mathit{P}}^{*}$ (present)	${\mathit{P}}^{*}$ (Timoshenko)	deviation
5.0	0.0012	0.007	0.007	$0.01\%$
140.0	0.0341	0.200	0.200	$0.02\%$
500.0	0.1209	0.714	0.714	$0.10\%$
2619.1	0.5477	3.742	3.679	$1.69\%$
7909.5	1.1240	11.299	10.648	$6.12\%$
19867.9	1.7428	28.383	25.352	$11.95\%$
60000.0	2.6949	85.714	71.769	$19.43\%$

Table 2. Clamped circular plate, large deflection. $P^{*}=p{a}^4/\left(E{h}^4\right)$. The growing deviation is expected: the reference expression is a one-term Galerkin approximation whose cubic coefficient is known to be low at large $\overline{w}$, whereas the present computation is a converged solution of the geometrically exact problem.

p [Pa]	$\overline{\mathit{w}}={\mathit{w}}_0/\mathit{h}$	${\mathit{P}}^{*}$ (present)	${\mathit{P}}^{*}$ (Timoshenko)	deviation
5.0	0.0012	0.007	0.007	$0.01\%$
140.0	0.0341	0.200	0.200	$0.02\%$
500.0	0.1209	0.714	0.714	$0.10\%$
2619.1	0.5477	3.742	3.679	$1.69\%$
7909.5	1.1240	11.299	10.648	$6.12\%$
19867.9	1.7428	28.383	25.352	$11.95\%$
60000.0	2.6949	85.714	71.769	$19.43\%$

10.6. Comparison with Exact and Published Solutions

The verifications so far are internal. This subsection compares the formulation against three independent reference solutions that it played no part in producing (Figure 5).

Figure 5. Comparison against exact reference solutions. (a) Axisymmetric free-vibration eigenvalues of a clamped circular plate against the exact Bessel-function values (dashed). (b) Linear deflection profile against Timoshenko and Woinowsky-Krieger’s closed form $w=p{(a^2-r^2)}^2/64B_f$. (c) Meridional resultant of a clamped spherical cap normalised by the membrane (Young–Laplace) value $-pR/2$; the three pressures collapse onto a single curve, confirming that the $7.4\%$ interior offset is a linear edge-restraint effect and not a nonlinear artefact.

Figure 5. Comparison against exact reference solutions. (a) Axisymmetric free-vibration eigenvalues of a clamped circular plate against the exact Bessel-function values (dashed). (b) Linear deflection profile against Timoshenko and Woinowsky-Krieger’s closed form $w=p{(a^2-r^2)}^2/64B_f$. (c) Meridional resultant of a clamped spherical cap normalised by the membrane (Young–Laplace) value $-pR/2$; the three pressures collapse onto a single curve, confirming that the $7.4\%$ interior offset is a linear edge-restraint effect and not a nonlinear artefact.

(a) Free vibration of a clamped circular plate.

The axisymmetric natural frequencies satisfy the classical characteristic equation [27,52]

$$J_0\left(\lambda \right)I_1\left(\lambda \right)+I_0\left(\lambda \right)J_1\left(\lambda \right)=0,{\lambda }^2=\omega a^2\sqrt{{\rho }_{0}h/B_f},$$

(63)

whose first three roots are ${\lambda }^2=10.2158,39.7711,89.1041$. Table 3 compares these with the eigenvalues of the discrete pencil $({\mathsf{K}}_T,\mathsf{M})$ at zero load. Four cubic Hermite elements already give the fundamental frequency to $4\times {10}^{-3}\%$. This is the only external validation of the mass matrix and of the tangent operator used for the stability analysis of Section 10.9, and it is therefore the comparison that most directly underwrites Figure 10.

Table 3. Axisymmetric vibration eigenvalues ${\lambda }^2=\omega a^2\sqrt{{\rho }_{0}h/B_f}$ of a clamped circular plate; relative error against the exact roots of (63) in parentheses. The error passes through a minimum near 16 elements and then grows: beyond that point the accuracy is limited not by the discretisation but by the finite-difference construction of the tangent operator, whose relative error ($\sim {10}^{-8}$) is amplified by the growing condition number of ${\mathsf{K}}_T$.

elements	${\mathit{\lambda }}_1^2$	${\mathit{\lambda }}_2^2$	${\mathit{\lambda }}_3^2$
4	$10.2163$ ($0.0042\%$)	$39.8431$ ($0.181\%$)	$90.1217$ ($1.142\%$)
8	$10.2160$ ($0.0022\%$)	$39.7760$ ($0.0123\%$)	$89.1765$ ($0.0812\%$)
16	$10.2166$ ($0.0077\%$)	$39.7724$ ($0.0032\%$)	$89.1098$ ($0.0063\%$)
32	$10.2188$ ($0.0289\%$)	$39.7748$ ($0.0092\%$)	$89.1083$ ($0.0047\%$)
64	$10.2246$ ($0.0857\%$)	$39.7820$ ($0.0273\%$)	$89.1157$ ($0.0130\%$)
exact	$10.2158$	$39.7711$	$89.1041$

Table 3. Axisymmetric vibration eigenvalues ${\lambda }^2=\omega a^2\sqrt{{\rho }_{0}h/B_f}$ of a clamped circular plate; relative error against the exact roots of (63) in parentheses. The error passes through a minimum near 16 elements and then grows: beyond that point the accuracy is limited not by the discretisation but by the finite-difference construction of the tangent operator, whose relative error ($\sim {10}^{-8}$) is amplified by the growing condition number of ${\mathsf{K}}_T$.

elements	${\mathit{\lambda }}_1^2$	${\mathit{\lambda }}_2^2$	${\mathit{\lambda }}_3^2$
4	$10.2163$ ($0.0042\%$)	$39.8431$ ($0.181\%$)	$90.1217$ ($1.142\%$)
8	$10.2160$ ($0.0022\%$)	$39.7760$ ($0.0123\%$)	$89.1765$ ($0.0812\%$)
16	$10.2166$ ($0.0077\%$)	$39.7724$ ($0.0032\%$)	$89.1098$ ($0.0063\%$)
32	$10.2188$ ($0.0289\%$)	$39.7748$ ($0.0092\%$)	$89.1083$ ($0.0047\%$)
64	$10.2246$ ($0.0857\%$)	$39.7820$ ($0.0273\%$)	$89.1157$ ($0.0130\%$)
exact	$10.2158$	$39.7711$	$89.1041$

(b) Linear plate deflection profile.

Beyond the centre deflection of Section 10.5, the whole computed profile is compared against $w\left(r\right)=p{(a^2-r^2)}^2/\left(64B_f\right)$ ([52], Art. 62). The two are graphically indistinguishable at $p=20$Pa; the maximum pointwise deviation is below $0.1\%$ of the centre deflection. This checks the bending operator and the clamped boundary conditions pointwise rather than at a single functional.

(c) Membrane limit of the spherical cap.

For a spherical shell under uniform pressure, membrane equilibrium (Young–Laplace) gives $N_s=N_{\theta }=-pR/2$. Figure 5(c) plots $N_s/(-pR/2)$ for three pressures spanning a factor of sixteen. The interior value is $1.0740$, $1.0750$ and $1.0792$ respectively—constant to better than $0.5\%$ across the load range, and dropping sharply only inside the edge layer of width $\sqrt{Rh}$. The formulation therefore reproduces the membrane state up to a fixed $7.4\%$ offset attributable to the immovable clamped edge, which for a shallow cap contributes a uniform additional compression that does not decay into the interior. That the offset is independent of p confirms it is a linear edge-restraint effect rather than an error in the pressure functional; the enclosed-volume functional (62) itself reproduces the exact cap volume $\pi (R{H}^2-H^3/3)$ to 8 significant figures.

Remark 15 

(What is not compared).  We deliberately do not overlay published critical pressures on Figure 7. The classical axisymmetric results of Budiansky [5] and the bifurcation results of Huang [18] and Weinitschke [53] are reported graphically as $p_{\mathrm{cr}}/p_{\mathrm{cl}}$ against a shallowness parameter whose normalisation differs between authors, and we have not had access to tabulated values. Digitising a published figure and presenting the result as a quantitative comparison would give a false impression of precision. The values in Table 4 should therefore be read as lying within the band those authors report, not as validated against them; the quantitative validations of this paper are those of Sections 10.1, 10.2, 10.5 and 10.6.

10.7. Mesh Convergence

Figure 6 shows h-convergence of the critical snap-through pressure for the $\mu =5$ cap of Section 10.8. The computed values are $p_{\mathrm{cr}}/p_{\mathrm{cl}}=0.624997,0.619907,0.618153,0.617324,0.616910$ for $5,10,20,40,80$ elements. Even five cubic Hermite elements are within $1.3\%$; the observed rate is consistent with the expected $O\left(h_e^2\right)$ for a limit-point load computed from a $C^1$ cubic discretisation. All subsequent results use 30 elements.

Figure 6. Mesh convergence of the critical pressure, $\mu =5$ cap; the 80-element value is taken as reference.

Figure 6. Mesh convergence of the critical pressure, $\mu =5$ cap; the 80-element value is taken as reference.

10.8. Snap-Through of Clamped Spherical Caps

We now apply the formulation to its intended regime. A clamped spherical cap of sphere radius $R=1$m and thickness $h=2$mm is loaded by uniform external pressure. Cap geometry is parameterised by the shallowness number of Budiansky [5] and Huang [18],

$$\mu =2{\left[3(1-{\nu }^2)\right]}^{1/4}\sqrt{H/h},$$

(64)

H being the rise; pressures are normalised by the classical buckling pressure of a complete sphere, $p_{\mathrm{cl}}=2E{h}^2/\left(R^2\sqrt{3(1-{\nu }^2)}\right)=338.9$kPa, the Zoelly result [55]. Three caps are considered, listed in Table 4.

Table 4. Cap geometries and computed critical loads. $p_{\mathrm{cr}}$ is the upper limit point, $p_{min}$ the following minimum (the lower limit point). The $\mu =3$ cap is too shallow to lose stability: its equilibrium path is monotone and its fundamental frequency never vanishes.

$\mathit{\mu }$	$\mathit{H}/\mathit{h}$	base radius a [mm]	${\mathit{p}}_{\mathbf{cr}}/{\mathit{p}}_{\mathbf{cl}}$	${\mathit{w}}_{\mathbf{cr}}/\mathit{h}$	${\mathit{p}}_{min}/{\mathit{p}}_{\mathbf{cl}}$	${\mathit{f}}_1$ at $\mathit{p}=0$ [Hz]
3	1.36	73.8	—(monotone)	—	—	1462.9
4	2.42	98.3	0.5648	1.01	0.4247	1237.5
5	3.78	122.8	0.6168	0.83	0.2846	1139.4

Table 4. Cap geometries and computed critical loads. $p_{\mathrm{cr}}$ is the upper limit point, $p_{min}$ the following minimum (the lower limit point). The $\mu =3$ cap is too shallow to lose stability: its equilibrium path is monotone and its fundamental frequency never vanishes.

$\mathit{\mu }$	$\mathit{H}/\mathit{h}$	base radius a [mm]	${\mathit{p}}_{\mathbf{cr}}/{\mathit{p}}_{\mathbf{cl}}$	${\mathit{w}}_{\mathbf{cr}}/\mathit{h}$	${\mathit{p}}_{min}/{\mathit{p}}_{\mathbf{cl}}$	${\mathit{f}}_1$ at $\mathit{p}=0$ [Hz]
3	1.36	73.8	—(monotone)	—	—	1462.9
4	2.42	98.3	0.5648	1.01	0.4247	1237.5
5	3.78	122.8	0.6168	0.83	0.2846	1139.4

Figure 7 shows the complete equilibrium paths. The transition between monotone and snapping behaviour occurs between $\mu =3$ and $\mu =4$, consistent with the classical picture [5,18,22,53]. The critical pressures, $0.56$ and $0.62$ times $p_{\mathrm{cl}}$, lie in the band reported for clamped shallow caps in the axisymmetric analyses of [5,7], though see Remark 15 on the limits of that statement and Section 10.13(ii) on axisymmetry.

Figure 8 makes the mechanism visible: the apex displacement at state F exceeds $2H$, so the shell has passed entirely through its own base plane. Rotations of the meridian normal exceed ${30}^{\circ }$ over much of the cap. This is precisely the regime in which the spurious strain of Figure 2 would dominate a linearised formulation.

Figure 9 isolates the term ${\widehat{N}}^{ab}{\tilde{b}}_{ab}$ identified in Remark 8. At the limit point the meridional resultant is compressive over the whole cap, and the geometric-stiffness density is large and negative near the clamped edge, where the bending boundary layer of width $\sim \sqrt{Rh}\approx 45$mm concentrates the curvature change. Omitting this term—or replacing $\widehat{N}$ by N — changes the computed critical pressure and destroys the correspondence between limit points and vanishing frequencies established next.

Figure 7. Equilibrium paths of clamped spherical caps under uniform pressure. Downward triangles mark upper limit points, upward triangles the following minima. The $\mu =3$ cap stiffens monotonically; $\mu =4$ and 5 exhibit classical snap-through, the unstable segment between the two limit points being traversable only under displacement control.

Figure 7. Equilibrium paths of clamped spherical caps under uniform pressure. Downward triangles mark upper limit points, upward triangles the following minima. The $\mu =3$ cap stiffens monotonically; $\mu =4$ and 5 exhibit classical snap-through, the unstable segment between the two limit points being traversable only under displacement control.

Figure 8. Deformation sequence for the $\mu =5$ cap. (a) Meridian profiles at the six marked states; the cap passes continuously from the initial convex shape (A) through the limit point (B) and the unstable branch (C, D) to a fully inverted, re-stiffening configuration (F). (b) Location of each state on the equilibrium path.

Figure 8. Deformation sequence for the $\mu =5$ cap. (a) Meridian profiles at the six marked states; the cap passes continuously from the initial convex shape (A) through the limit point (B) and the unstable branch (C, D) to a fully inverted, re-stiffening configuration (F). (b) Location of each state on the equilibrium path.

Figure 9. Stress resultants and geometric-stiffness density for the $\mu =5$ cap at the limit point (B), the minimum (D) and the fully inverted state (F). Solid lines: meridional; dashed: circumferential. The right panel shows ${\widehat{N}}^{ab}{\tilde{b}}_{ab}$ of (40), the term responsible for the instability.

Figure 9. Stress resultants and geometric-stiffness density for the $\mu =5$ cap at the limit point (B), the minimum (D) and the fully inverted state (F). Solid lines: meridional; dashed: circumferential. The right panel shows ${\widehat{N}}^{ab}{\tilde{b}}_{ab}$ of (40), the term responsible for the instability.

10.9. Dynamic Signature of the Instability

Because the formulation is variational, the second variation (55) evaluated at an equilibrium is the tangent stiffness of the linearised dynamics, and the generalised eigenproblem $\left({\mathsf{K}}_T-{\omega }^2\mathsf{M}\right)\mathit{\psi }=\mathbf{0}$—with $\mathsf{M}$ the consistent mass matrix from (35)—gives the natural frequencies of small vibration about that state. Stability and dynamics are therefore not separate calculations.

Figure 10(a) is the sharpest available check on the internal consistency of the whole formulation. The frequency computed from ${\mathsf{K}}_T$ crosses zero at $w_0/h=1.013$ and $2.825$ for $\mu =4$, and at $0.833$ and $5.164$ for $\mu =5$; the limit points located independently from the load–deflection curve are at $w_0/h=1.013$ and $2.878$, and $0.833$ and $5.164$. The two determinations agree to the resolution of the path discretisation, as they must if the tangent operator is genuinely the second variation of the same energy whose stationarity defines the path. For $\mu =3$ the frequency dips but never reaches zero, in agreement with the monotone path.

Figure 10. (a) Fundamental frequency along the equilibrium paths, plotted with the sign of the corresponding eigenvalue of ${\mathsf{K}}_T$: it passes through zero exactly at the limit points (triangles) and is imaginary (plotted negative) on the unstable branch. The $\mu =3$ cap never loses positive-definiteness. (b) Southwell-type frequency–load relation for $\mu =5$: ${(f_1/f_1^{\left(0\right)})}^2$ is close to linear in $p/p_{\mathrm{cr}}$ at low load and departs as the limit point is approached.

Figure 10. (a) Fundamental frequency along the equilibrium paths, plotted with the sign of the corresponding eigenvalue of ${\mathsf{K}}_T$: it passes through zero exactly at the limit points (triangles) and is imaginary (plotted negative) on the unstable branch. The $\mu =3$ cap never loses positive-definiteness. (b) Southwell-type frequency–load relation for $\mu =5$: ${(f_1/f_1^{\left(0\right)})}^2$ is close to linear in $p/p_{\mathrm{cr}}$ at low load and departs as the limit point is approached.

Panel (b) shows the frequency–load relation used in non-destructive buckling estimation. The near-linearity of ${(f_1/f_1^{\left(0\right)})}^2$ against $p/p_{\mathrm{cr}}$ at low load is the Southwell/Lurie behaviour that underwrites vibration-based load estimation; the pronounced departure above $p/p_{\mathrm{cr}}\approx 0.5$ quantifies the error incurred by extrapolating it, which for this cap would over-predict the critical pressure by about $35\%$.

10.10. Comparison with Published Finite-Element Results

The comparisons of Section 10.6 are against closed-form solutions. This subsection compares against published numerical results obtained with entirely different shell formulations. We use the benchmark suite of Simo, Fox and Rifai [43], which reports tabulated values—not only figures—and which was explicitly documented by its authors to promote comparison with other formulations.

Roll-up of a clamped strip ([43], Section 6.1.1).

An initially flat strip ($E=12\times {10}^6$, $\nu =0$, $L=10$, $w=1$, $t=0.1$) under an end moment $M=2\pi EI/L$ rolls into a complete circle of radius $\rho =EI/M=L/2\pi =1.59155$. Evaluating Definition 2 on the exact rolled-up configuration gives $\parallel {\mathit{E}}^{\left(0\right)}\parallel \le 1.1\times {10}^{-11}$ — the deformation is recognised as an exact isometry—and $K_{11}^{\left(1\right)}=0.62831858$ against the exact $1/\rho =0.62831853$. The residual is the finite-difference error of the surface derivatives used in this particular check, not of the theory. This is the shell-theoretic content of the benchmark: a ${360}^{\circ }$ rotation with identically zero membrane strain. In [43] the complete $2\pi$ roll-up is driven in a single load step on a 25-element mesh, which is the practical demonstration that the director parametrisation is singularity-free; the check performed here is of the strain measures alone and involves no incremental solution.

Between them the companion papers report a much larger suite. Part II ([42], Sections 7.1–7.8) covers the linear tests—constant tension, bending and in-plane twist patch tests on regular and skewed meshes, Cook’s membrane, the singularity-dominated rhombic plate, the ${90}^{\circ }$ twisted ribbon, the pinched hemisphere both with an ${18}^{\circ }$ hole and closed, the Scordellis–Lo barrel vault, the pinched cylinder with end diaphragms, and two mesh-distortion studies (a simply supported plate distorted to the point of a negative Jacobian, and a cantilever of ten arbitrarily skewed elements)—while Part III adds the nonlinear ones: in-plane bending on a deliberately distorted mesh, torsion of a flat strip through ${180}^{\circ }$, the finite-strain pinched hemisphere, buckling of a simply-supported plate, an L-shaped frame, and snap-through of a hinged cylindrical panel ([43], Sections 6.1–6.3). None of these is axisymmetric—in the Part II suite the hemisphere is the only axisymmetric geometry, and its loading, four alternating point forces at ${90}^{\circ }$, is not — so none is accessible to the reduction of Section 10.3, and we make no comparison against them. Part II does tabulate its results ([42], Tables 9–17) rather than reporting them graphically, so the obstacle is the kinematics of the benchmarks and not the availability of numbers.

One benchmark in Part III is an exception worth recording. The collapse of a clamped rubber sphere under a point load ([43], Section 6.3.3) is axisymmetric, and [43] compute it with transverse shear suppressed by a penalty precisely in order to match the shear-rigid theory of Taber [50], against which experimental data are reported. It is therefore a rare instance of a published, experimentally anchored, shear-rigid, axisymmetric, large-deflection benchmark—a natural target for the present formulation. We have not attempted it, because with $h/R=0.17$ and $\nu =0.5$ it lies well outside the thin, small-strain regime in which the Saint-Venant–Kirchhoff closure of Section 7 is defensible (Section 10.13(iii)). We flag it as the obvious next validation step. Table 2 of [42] is used separately in Section 7.2.

Axially compressed clamped cylinder ([43], Table 6.2.3).

A cylinder with $E=567$, $\nu =0.3$, $R=100$, $L=35.95$, $t=0.247$ ($R/t\approx 405$, $L/R=0.36$) is compressed axially between clamped ends. Simo et al. report that for clamped ends the critical mode is axisymmetric, which is what makes this benchmark accessible to the present axisymmetric reduction. Loads are normalised by $P_{\mathrm{cl}}=E{t}^2/\left(R\sqrt{3(1-{\nu }^2)}\right)$ ([51], p. 465).

Following [43] we perform a linearised buckling analysis: the material stiffness ${\mathsf{K}}_m={\mathsf{K}}_T{|}_{\delta =0}$ and the geometric stiffness ${\mathsf{K}}_g=\partial {\mathsf{K}}_T/\partial \delta$ are extracted from the tangent operator at imposed end shortening $\delta$, and the eigenproblem ${\mathsf{K}}_m\mathit{\psi }=-\Lambda {\mathsf{K}}_g\mathit{\psi }$ is solved. Table 5 and Figure 11 give the result.

Figure 11. Comparison with ([43], Table 6.2.3) for the clamped axially compressed cylinder. (a) The present value lies within the spread of the three published formulations. (b) Mesh convergence: the present value is converged to four figures from 40 elements onwards. (c) The computed critical mode is the axisymmetric ring-wave pattern reported in ([43], Figure 6.2.4(b)).

Figure 11. Comparison with ([43], Table 6.2.3) for the clamped axially compressed cylinder. (a) The present value lies within the spread of the three published formulations. (b) Mesh convergence: the present value is converged to four figures from 40 elements onwards. (c) The computed critical mode is the axisymmetric ring-wave pattern reported in ([43], Figure 6.2.4(b)).

Table 5. Normalised critical load $P_{\mathrm{cr}}/P_{\mathrm{cl}}$ for the clamped axially compressed cylinder of ([43], Table 6.2.3). Case (c): imposed end shortening; case (d): imposed compressive load. Reference values are quoted directly from ([43], Table 6.2.3); 4-RSDS is the resultant-stress degenerated element of Liu et al. [29] and 4-SRI the selectively-integrated degenerated element of Hughes and Liu [21].

Formulation	kinematics	case (c)	case (d)
4-RSDS [29]	degenerated solid, shear-flexible	$1.0369$	$1.0379$
4-SRI [21]	degenerated solid, shear-flexible	$1.1378$	$1.1361$
Simo, Fox & Rifai [43]	one-director Cosserat, shear-flexible	$1.0837$	$1.0833$
present	Kirchhoff–Love, shear-rigid	$\mathbf{1}.\mathbf{0981}$	—

Table 5. Normalised critical load $P_{\mathrm{cr}}/P_{\mathrm{cl}}$ for the clamped axially compressed cylinder of ([43], Table 6.2.3). Case (c): imposed end shortening; case (d): imposed compressive load. Reference values are quoted directly from ([43], Table 6.2.3); 4-RSDS is the resultant-stress degenerated element of Liu et al. [29] and 4-SRI the selectively-integrated degenerated element of Hughes and Liu [21].

Formulation	kinematics	case (c)	case (d)
4-RSDS [29]	degenerated solid, shear-flexible	$1.0369$	$1.0379$
4-SRI [21]	degenerated solid, shear-flexible	$1.1378$	$1.1361$
Simo, Fox & Rifai [43]	one-director Cosserat, shear-flexible	$1.0837$	$1.0833$
present	Kirchhoff–Love, shear-rigid	$\mathbf{1}.\mathbf{0981}$	—

The present value $1.0981$ differs from Simo et al. by $+1.33\%$, lies $3.5\%$ below 4-SRI and $5.9\%$ above 4-RSDS—that is, comfortably inside the $9.7\%$ spread of the three published results, and closest to the geometrically exact formulation of the three. Two remarks are in order.

(i)

The residual $+1.33\%$ against [43] has an identifiable origin: the present theory is Kirchhoff–Love and therefore shear-rigid, whereas all three reference formulations admit transverse shear. A shear-rigid model is stiffer and must over-predict the buckling load; the sign of the discrepancy is therefore the expected one, and its magnitude is consistent with $O(t^2/{\lambda }_b^2)$ where ${\lambda }_b\approx \pi \sqrt{Rt}/{\left[3(1-{\nu }^2)\right]}^{1/4}\approx 12$ is the axial buckle half-wavelength.

(ii)

We compare only case (c). Cases (a) and (b) of ([43], Table 6.2.3) have free ends, for which the reported critical loads fall to $0.5402$ and $0.5055$—about half the clamped value, in agreement with the analytical result of Hoff and Soong [17]—in a mode carrying two circumferential sine-waves and none in the axial direction. That mode is non-axisymmetric and therefore outside the reach of the present reduction (Section 10.13(ii)). We do not report a value for it.

(iii)

The axisymmetry of the mode is a property of the linearised analysis, not of the structure. When the eigenvalue problem is re-solved at the computed bifurcation point instead of at the reference configuration, the clamped cylinder’s critical mode is no longer axisymmetric but reticulated, at $P/P_{\mathrm{cl}}=1.04$ ([43], Section 6.2.4). What Table 5 compares is thus the linearised buckling load as computed by four formulations under one and the same restrictive assumption, not the true critical load of the cylinder. Simo et al. draw the explicit moral that linearised buckling results “should be interpreted with caution”; we adopt it. The same caveat attaches to the caps of Section 10.8—see Section 10.13(ii).

10.11. Asymptotic Convergence of the Newton Iteration

Remark 13 claims that the tangent used here is the exact derivative of the discrete residual. That claim has an observable consequence, and it is the sharpest single test of a nonlinear implementation available: an exact tangent forces the Newton error to contract quadratically, and nothing else does. Simo, Fox and Rifai report this rate for their formulation in ([43], Tables 6.4.1–6.4.3), which makes it directly comparable.

We take one large load step from the undeformed state and iterate with undamped Newton—no line search, no step limiting—recording the Euclidean norm of the out-of-balance force $\parallel R_k\parallel$ and the energy norm $e_k=|\Delta U_k·R_k|$. Table 6 gives the two sequences and Figure 12 plots them.

The energy norm falls by eighteen decades in eight iterations for the plate and by twenty-one in five for the cap. Fitting $log{e}_{k+1}=rlog{e}_k$ over the asymptotic range gives $r=1.9$–$2.0$ in both cases, and Figure 12(b) shows the points lying on the slope-2 reference line rather than the slope-1 line that an approximate tangent would produce. The behaviour matches that of ([43], Tables 6.4.1–6.4.3) closely enough that the two sets of sequences are hard to separate in Figure 12(a), which is the expected outcome: quadratic convergence is a property of exact linearisation, not of a particular shell theory, and its appearance in both is evidence that both linearisations are in fact exact.

Table 6. Newton convergence in a single load step with the consistent tangent. Left: clamped circular plate (30 elements) at $p=3\times {10}^4$Pa, i.e., apex deflection $w_c/h=2.07$, deep in the nonlinear regime. Right: clamped spherical cap ($\mu =5$, 30 elements) at $p=0.50p_{\mathrm{cl}}$. Both bottom out at the round-off floor of the residual evaluation, $\parallel R\parallel \sim {10}^{-6}$ in the units of the problem.

	clamped plate			spherical cap
k	$\parallel R_k\parallel$	$e_k$	k	$\parallel R_k\parallel$	$e_k$
0	$1.94\times {10}^{+2}$	$2.30\times {10}^{+0}$	0	$1.66\times {10}^{+3}$	$2.55\times {10}^{+0}$
1	$8.14\times {10}^{+4}$	$3.03\times {10}^{+1}$	1	$2.81\times {10}^{+3}$	$8.13\times {10}^{-2}$
2	$8.58\times {10}^{+3}$	$4.66\times {10}^{+0}$	2	$1.78\times {10}^{+2}$	$6.51\times {10}^{-4}$
3	$6.32\times {10}^{+3}$	$7.82\times {10}^{-1}$	3	$3.80\times {10}^{+0}$	$3.17\times {10}^{-7}$
4	$1.66\times {10}^{+3}$	$9.65\times {10}^{-2}$	4	$3.01\times {10}^{-3}$	$1.90\times {10}^{-13}$
5	$4.90\times {10}^{+2}$	$4.06\times {10}^{-3}$	5	$1.11\times {10}^{-6}$	$8.75\times {10}^{-22}$
6	$2.83\times {10}^{+1}$	$1.27\times {10}^{-5}$
7	$1.05\times {10}^{-1}$	$1.54\times {10}^{-10}$
8	$1.57\times {10}^{-6}$	$9.10\times {10}^{-19}$

Table 6. Newton convergence in a single load step with the consistent tangent. Left: clamped circular plate (30 elements) at $p=3\times {10}^4$Pa, i.e., apex deflection $w_c/h=2.07$, deep in the nonlinear regime. Right: clamped spherical cap ($\mu =5$, 30 elements) at $p=0.50p_{\mathrm{cl}}$. Both bottom out at the round-off floor of the residual evaluation, $\parallel R\parallel \sim {10}^{-6}$ in the units of the problem.

	clamped plate			spherical cap
k	$\parallel R_k\parallel$	$e_k$	k	$\parallel R_k\parallel$	$e_k$
0	$1.94\times {10}^{+2}$	$2.30\times {10}^{+0}$	0	$1.66\times {10}^{+3}$	$2.55\times {10}^{+0}$
1	$8.14\times {10}^{+4}$	$3.03\times {10}^{+1}$	1	$2.81\times {10}^{+3}$	$8.13\times {10}^{-2}$
2	$8.58\times {10}^{+3}$	$4.66\times {10}^{+0}$	2	$1.78\times {10}^{+2}$	$6.51\times {10}^{-4}$
3	$6.32\times {10}^{+3}$	$7.82\times {10}^{-1}$	3	$3.80\times {10}^{+0}$	$3.17\times {10}^{-7}$
4	$1.66\times {10}^{+3}$	$9.65\times {10}^{-2}$	4	$3.01\times {10}^{-3}$	$1.90\times {10}^{-13}$
5	$4.90\times {10}^{+2}$	$4.06\times {10}^{-3}$	5	$1.11\times {10}^{-6}$	$8.75\times {10}^{-22}$
6	$2.83\times {10}^{+1}$	$1.27\times {10}^{-5}$
7	$1.05\times {10}^{-1}$	$1.54\times {10}^{-10}$
8	$1.57\times {10}^{-6}$	$9.10\times {10}^{-19}$

Figure 12. Asymptotic convergence of the Newton iteration with the consistent tangent. (a) Energy norm against iteration number, normalised by its value at $k=0$; the dashed curves are the corresponding sequences reported by Simo, Fox and Rifai ([43], Tables 6.4.1–6.4.3), normalised the same way. (b) The same data as $e_{k+1}$ against $e_k$: a quadratically convergent method must fall on a line of slope 2, a method with an inexact tangent on a line of slope 1. The present sequences track slope 2 until the round-off floor of the residual evaluation is reached.

Figure 12. Asymptotic convergence of the Newton iteration with the consistent tangent. (a) Energy norm against iteration number, normalised by its value at $k=0$; the dashed curves are the corresponding sequences reported by Simo, Fox and Rifai ([43], Tables 6.4.1–6.4.3), normalised the same way. (b) The same data as $e_{k+1}$ against $e_k$: a quadratically convergent method must fall on a line of slope 2, a method with an inexact tangent on a line of slope 1. The present sequences track slope 2 until the round-off floor of the residual evaluation is reached.

Two caveats. First, the first one or two iterations are not in the asymptotic regime and need not contract at all—the plate residual rises by two orders of magnitude at $k=1$, and the same non-monotonicity is visible in ([43], Table 6.4.3), where the residual rises at the first iteration of a one-step collapse computation. Quadratic convergence is an asymptotic statement and nothing more. Second, undamped Newton from an undeformed state has a finite basin: at $p\gtrsim 0.7p_{\mathrm{cl}}$ the cap iteration diverges, and the damped scheme used for the production runs of Section 10.8 is required. This is a property of Newton’s method, not of the tangent.

10.12. Summary of Verification

Table 7 collects every quantitative check performed, with its reference solution and source. The distinction between the three columns matters: exact references are closed-form or machine-computable to arbitrary precision; published references are approximate solutions from the literature, so a discrepancy does not by itself indicate an error in the present work; internal checks test consistency between parts of this paper and cannot detect an error common to both.

Table 7. Complete verification record. “E” = exact reference, “P” = published approximate reference, “I” = internal consistency check.

Quantity	Reference		Agreement
Strain under finite rigid rotation	${\mathit{E}}^{\left(0\right)}={\mathit{K}}^{\left(1\right)}=0$, Theorem 2	E	$6\times {10}^{-18}$–$7\times {10}^{-17}$, all $\theta \le {180}^{\circ }$
Thickness truncation	$O\left(z^2\right)$, Corollary 1	E	measured slope $2.000$
Discrete reference state	zero energy, zero residual	E	$0.0$ (exact, double precision)
Plate centre deflection, small p	$p{a}^4/64B_f$ ([52], Art. 62)	E	$1\times {10}^{-4}$ relative
Plate deflection profile	$p{(a^2-r^2)}^2/64B_f$ [52]	E	$<0.1\%$ pointwise
Plate vibration ${\lambda }_1^2$	$10.2158$, roots of (63) [27]	E	$0.004\%$ with 4 elements
Plate vibration ${\lambda }_2^2,{\lambda }_3^2$	$39.7711$, $89.1041$ [27]	E	$0.003\%$, $0.006\%$ (16 el.)
Cap enclosed volume	$\pi (R{H}^2-H^3/3)$	E	8 significant figures
Cap interior resultant	$-pR/2$ (Young–Laplace)	E	$1.0740$–$1.0792$ over $16\times$ load range
Plate large deflection, $\overline{w}\le 0.12$	$5.86\overline{w}+2.86{\overline{w}}^3$ ([52], Art. 97)	P	$0.01$–$0.10\%$
Plate large deflection, $\overline{w}\le 2.7$	as above	P	$1.7$–$19\%$ (reference is a one-term Galerkin approximation)
Cap $p_{\mathrm{cr}}/p_{\mathrm{cl}}$	[5,18,53]	P	not compared—see Remark 15
Strip roll-up, ${360}^{\circ }$	$\parallel {\mathit{E}}^{\left(0\right)}\parallel =0$, $K=M/EI$ ([43], Section 6.1.1)	E	$1.1\times {10}^{-11}$; $8\times {10}^{-6}$ rel.
Higher-order stiffness	([42], Table 2)	P	structure recovered (hyperelastic, same ${\mathit{d}}^0$ dependence); coupling differs by $-\frac{3}{2}$ after conversion to bare resultants —Remark 11
Cylinder $P_{\mathrm{cr}}/P_{\mathrm{cl}}$	$1.0837$ ([43], Table 6.2.3)	P	$1.0981$ ($+1.33\%$); within the $9.7\%$ spread of three published formulations
Tangent operator	finite difference of the residual	I	$7.6\times {10}^{-8}$; exactly symmetric
Newton asymptotic rate	2, as for any exact tangent; cf. ([43], Tables 6.4.1–6.4.3)	I	measured $1.9$–$2.0$ over $\sim 20$ decades
Limit points vs. $f_1=0$	must coincide if variational	I	agree to path resolution
Mesh convergence of $p_{\mathrm{cr}}$	80-element value	I	$1.3\%$ with 5 elements

Table 7. Complete verification record. “E” = exact reference, “P” = published approximate reference, “I” = internal consistency check.

Quantity	Reference		Agreement
Strain under finite rigid rotation	${\mathit{E}}^{\left(0\right)}={\mathit{K}}^{\left(1\right)}=0$, Theorem 2	E	$6\times {10}^{-18}$–$7\times {10}^{-17}$, all $\theta \le {180}^{\circ }$
Thickness truncation	$O\left(z^2\right)$, Corollary 1	E	measured slope $2.000$
Discrete reference state	zero energy, zero residual	E	$0.0$ (exact, double precision)
Plate centre deflection, small p	$p{a}^4/64B_f$ ([52], Art. 62)	E	$1\times {10}^{-4}$ relative
Plate deflection profile	$p{(a^2-r^2)}^2/64B_f$ [52]	E	$<0.1\%$ pointwise
Plate vibration ${\lambda }_1^2$	$10.2158$, roots of (63) [27]	E	$0.004\%$ with 4 elements
Plate vibration ${\lambda }_2^2,{\lambda }_3^2$	$39.7711$, $89.1041$ [27]	E	$0.003\%$, $0.006\%$ (16 el.)
Cap enclosed volume	$\pi (R{H}^2-H^3/3)$	E	8 significant figures
Cap interior resultant	$-pR/2$ (Young–Laplace)	E	$1.0740$–$1.0792$ over $16\times$ load range
Plate large deflection, $\overline{w}\le 0.12$	$5.86\overline{w}+2.86{\overline{w}}^3$ ([52], Art. 97)	P	$0.01$–$0.10\%$
Plate large deflection, $\overline{w}\le 2.7$	as above	P	$1.7$–$19\%$ (reference is a one-term Galerkin approximation)
Cap $p_{\mathrm{cr}}/p_{\mathrm{cl}}$	[5,18,53]	P	not compared—see Remark 15
Strip roll-up, ${360}^{\circ }$	$\parallel {\mathit{E}}^{\left(0\right)}\parallel =0$, $K=M/EI$ ([43], Section 6.1.1)	E	$1.1\times {10}^{-11}$; $8\times {10}^{-6}$ rel.
Higher-order stiffness	([42], Table 2)	P	structure recovered (hyperelastic, same ${\mathit{d}}^0$ dependence); coupling differs by $-\frac{3}{2}$ after conversion to bare resultants —Remark 11
Cylinder $P_{\mathrm{cr}}/P_{\mathrm{cl}}$	$1.0837$ ([43], Table 6.2.3)	P	$1.0981$ ($+1.33\%$); within the $9.7\%$ spread of three published formulations
Tangent operator	finite difference of the residual	I	$7.6\times {10}^{-8}$; exactly symmetric
Newton asymptotic rate	2, as for any exact tangent; cf. ([43], Tables 6.4.1–6.4.3)	I	measured $1.9$–$2.0$ over $\sim 20$ decades
Limit points vs. $f_1=0$	must coincide if variational	I	agree to path resolution
Mesh convergence of $p_{\mathrm{cr}}$	80-element value	I	$1.3\%$ with 5 elements

Remark 16 

(On what the reader should conclude). The exact checks establish that the kinematics, the bending operator, the mass matrix, the pressure functional and the tangent operator are all correct, in the strong sense that each reproduces an independently known answer to several significant figures. What they do not establish is that the snap-through pressures of Table 4 agree with the literature, because we have not been able to make that comparison quantitatively (Remark 15). A reader who needs those numbers for design purposes should treat them as the output of a verified code applied to a problem whose published solutions we could not access in tabulated form, and not as validated predictions.

10.13. Limitations

Three restrictions of the case study should be stated plainly.

(i)

Deeper caps. For $\mu \gtrsim 6$ the axisymmetric equilibrium branch develops a fold in the apex deflection and in the enclosed volume near the upper limit point, so that neither load, apex-displacement nor volume control is a globally valid parameterisation. Tracing those branches requires true arclength continuation with adaptive step control. We have therefore restricted the reported results to $\mu \le 5$, where the paths are traversed without a single failed increment.

(ii)

Axisymmetry. The computation constrains deformation to be axisymmetric. Real spherical caps frequently buckle first into a non-axisymmetric mode [18,20,53], so the loads in Table 4 are upper bounds on the true critical pressures. Detecting those modes requires a Fourier decomposition in $\theta$ and a bifurcation check on each harmonic; the theory of Section 2–Section 9 supports this without modification, but the implementation does not. The restriction is not peculiar to caps: the axially compressed cylinder of Section 10.10 behaves the same way, its linearised clamped mode being axisymmetric while the mode computed at the bifurcation point is not ([43], Section 6.2.4).

(iii)

Constitutive scope. A Saint-Venant–Kirchhoff law is used, which is appropriate for the small-strain/large-rotation regime of these examples ($max\parallel {\mathit{E}}^{\left(0\right)}\parallel \approx 3\times {10}^{-3}$) but is not a sound finite-strain model. Nothing in Section 2–Section 5 depends on this choice; only (45) and (51) would change; see [2,47] for finite-strain shell constitutive theory.

11. Conclusions

Treating the shell body as a fibre bundle and its geometry through a moving coframe is not a change of notation. It changes what has to be assumed and what follows.

• Compatibility becomes an identity. The Gauss and Codazzi–Mainardi equations are the components of the curvature 2-form ${{\mathsf{\Omega }}^I}_J$, which vanishes because the reference configuration is embedded in flat space (Theorem 1). They need not be imposed or verified.
• Objectivity becomes exact. Because the strain measures are differences of pulled-back metrics, invariance under $\mathrm{SO}\left(3\right)⋉{\mathbb{R}}^3$ is an algebraic identity (Theorem 2), confirmed numerically at the ${10}^{-17}$ level over the full range of rotation angles. Table 1 quantifies what is lost otherwise: $207$MPa of fictitious stress at $5^{\circ }$.
• Dimensional reduction becomes fibre integration. Carrying the exact shifter determinant ${\mu }_0=1+2Hz+K{z}^2$ through the fibre integral yields resultant stiffnesses (34) including a curvature-induced membrane–bending coupling $\mathsf{B}=H{h}^3/6$ that the usual derivation discards without comment.
• The balance laws acquire their missing terms. Deriving Theorem 3 from the variational statement rather than by analogy exposes the effective resultant ${\widehat{N}}^{ab}=N^{ab}+M^{ac}{{\tilde{b}}_c}^b$, the shear–curvature coupling in the tangential equation, and the Kirchhoff corner forces (43). Figure 9 shows that the geometric-stiffness density ${\widehat{N}}^{ab}{\tilde{b}}_{ab}$ is precisely the instability mechanism.
• Stability and dynamics coincide. Because the theory is variational, the tangent operator is the second variation, and the vanishing of the fundamental frequency at the limit points (Figure 10) is automatic rather than coincidental. This is the strongest single check on the consistency of the formulation, and it is passed to the resolution of the path discretisation.

A further caveat deserves emphasis: the loads computed here are those of the perfect shell. Spherical shells are notoriously imperfection-sensitive [7,20,22], and experimentally observed buckling pressures fall well below the perfect-shell values. The present formulation is well suited to imperfection studies, since an imperfect reference surface enters only through ${\mathit{\phi }}_0$ and requires no change to any equation.

The natural continuations are non-axisymmetric bifurcation analysis via harmonic decomposition, finite-strain constitutive laws (for which only Section 7 changes), and isogeometric discretisation [19,23], which supplies the $C^1$ regularity that Remark 4 shows to be intrinsic to the Kirchhoff–Love constraint rather than a numerical convenience.