A Geometrically Exact Formulation of Elastic Shell Dynamics Based on Fiber Bundles and Differential Forms

1. Introduction

1.1. Historical Evolution and Geometric Awakening of Shell Theory

The mathematical theory of elastic shells is one of the oldest, most complex, and most physically captivating branches of continuum mechanics. Fundamentally, shell mechanics deals with how to reduce three-dimensional elastic equations to a two-dimensional surface, while preserving the additional stiffness and complex geometric coupling effects brought about by the curvature of the surface. The history of this field is a history of constant struggle against coordinate complexity, culminating in a geometric awakening.

The foundation of modern shell theory can be traced back to the mid-19th century. In 1850, Gustav Kirchhoff proposed two fundamental assumptions when studying elastic plates: the hypothesis of normal preservation (normals to the mid-surface before deformation remain normal and unstretched after deformation) and the hypothesis of negligible transverse normal stress [1]. These assumptions effectively reduced the number of unknown fields, allowing the three-dimensional elastic equations to be integrated through the thickness. Subsequently, Augustus Edward Hough Love extended this framework to surfaces with initial curvature in 1888, forming the renowned Kirchhoff-Love shell theory [2]. For much of the 20th century, this theory dominated civil, aerospace, and ship engineering. Scholars such as Love, Flügge [3], Novozhilov [5], Goldenveizer [6], and Kraus [4] derived extremely complex equations of motion in general curvilinear coordinate systems.

However, classical derivations face serious mathematical obstacles. In a general curvilinear coordinate system $({\alpha }^1,{\alpha }^2,z)$, the metric tensor $\mathit{g}$ of the shell depends explicitly on the thickness coordinate z. To derive equilibrium equations, one must calculate the covariant derivative of the stress tensor ${\nabla }_i{\sigma }^{ij}$, which requires introducing up to 18 independent Christoffel symbols of the second kind ${\mathsf{\Gamma }}_{jk}^i$. The expressions for these symbols are extremely lengthy and filled with nonlinear terms. When reducing the three-dimensional equilibrium equations via integration, a large number of cumbersome partial integrations and index contractions are involved. As pointed out by Green and Zerna [7], this coordinate-dependent derivation not only obscures the physical and geometric significance behind the equations but is also prone to algebraic errors. To simplify calculations, engineers have had to rely on special coordinate systems (such as lines of curvature coordinates), which compromises the universality of the theory.

In the 1970s, continuum mechanics underwent a profound paradigm shift. Élie Cartan’s method of moving frames [8] was introduced into mechanics. Cartan pointed out that by choosing a local orthonormal basis, the matrix of connection 1-forms is antisymmetric (${\omega }_{ij}=-{\omega }_{ji}$), which immediately halves the number of independent connection components. More importantly, exterior differential forms (including the exterior derivative $\mathrm{d}$, wedge product ∧, and interior product ${\iota }_X$) provide a completely coordinate-free calculus tool. Misner, Thorne, and Wheeler [9] demonstrated the powerful utility of these tools in general relativity. In the field of mechanics, Marsden and Hughes’ pioneering work Mathematical Foundations of Elasticity [10] systematically reformulated nonlinear elasticity using manifolds, tangent bundles, and Lie derivatives. They proved that the strain tensor is essentially the perturbation of the displacement flow on the Riemannian metric, and that equilibrium laws can be rigorously established via variational principles in a coordinate-free framework. Epstein [11] further extended these ideas to material defects and non-Euclidean elasticity, using principal bundles and connections to characterize the geometric structure of continua.

In recent years, important breakthroughs have been made in geometric methods based on Cartan’s moving frame method combined with Finite Element Exterior Calculus (FEEC). Dhas, N, Roy, and Reddy (2022) [23] systematically introduced Cartan’s moving frames into the mixed variational principle of nonlinear elasticity for the first time, interpreting physical quantities such as stress and deformation gradient as differential forms, and achieving geometrically exact discretization using FEEC. Subsequently, Jamun Kumar et al. extended this framework to Kirchhoff rods (2022) [24] and Kirchhoff shells (2024) [25], successfully constructing locking-free $C^0$ continuous finite element schemes by introducing independent orthogonal frame fields and the Hu-Washizu functional. Meanwhile, Rashad et al. (2023) [26] provided a unified exterior calculus formulation for intrinsic nonlinear elasticity starting from bundle-valued differential forms, revealing the principal bundle and de Rham complex structures underlying the theory. These advances indicate that combining fiber bundles, differential forms, and FEEC is a powerful tool for developing the next generation of geometrically exact, coordinate-free, and high-precision shell elements. This paper systematically constructs a geometrically exact shell dynamic theory based on fiber bundles and differential forms against this background, laying a foundation for subsequent numerical implementation.

1.2. The Necessity of Geometrically Exact Theory

With the development of modern engineering and physics, classical small-deformation linear shell theory can no longer meet demands. In fields such as flexible electronics [18], soft robotics, cell membrane mechanics [19], and nanomaterials (e.g., graphene [20]), shell structures often undergo large displacements and rigid body rotations, even though their intrinsic strains may remain small. Under these "large rotation, small strain" conditions, if linear kinematics are adopted, ignoring higher-order rotation terms, the structure will produce non-physical "spurious strains" under pure rigid body rotation, leading to completely erroneous stress calculations and numerical divergence.

To address this issue, Simo and Fox [14,30] proposed the "Geometrically Exact Shell Theory" in 1989. This theory retains full nonlinearity in kinematics, accurately describing arbitrarily large displacements and rotations, ensuring the invariance of rigid body motion. Subsequently, extensive research has been dedicated to embedding geometrically exact theory into finite element computations [14,15,30]. However, traditional geometrically exact theories are still mostly built within the framework of matrix vectors and tensor indices. For shells with complex initial curvature and non-Euclidean effects, index notation remains cumbersome.

In the development of the aforementioned geometrically exact shell theory, the seminal contributions of Simo and Fox (1989) and their subsequent work (particularly Fox’s doctoral dissertation [29,30]) laid the foundation. They refined shell kinematics into a combination of the mid-surface position vector and a director field vector by introducing the single-director field variable pair $(\mathit{\phi },\mathit{d})$, thereby completely avoiding explicit operations of covariant derivatives and Christoffel symbols at the kinematic level; all field quantities are handled directly in vector space. This framework, due to its great simplicity and natural compatibility with finite element discretization, has become an important cornerstone of modern nonlinear shell numerical analysis.

However, although Simo-Fox’s "no complex geometry" framework is computationally convenient, it largely "implicitizes" the curvature and torsion deformation information of the surface into the gradient of the director field $\mathit{d}$. While efficient, this treatment makes the coordinate-independence of the theory less explicit and difficult to directly generalize to cases with non-trivial topology or non-Euclidean reference metrics (e.g., growth, defect-induced residual stress).

The starting point of this paper is complementary to the aforementioned classic work: we do not shy away from geometry but actively make it explicit. Specifically, we reintroduce Cartan’s moving frames, connection 1-forms on fiber bundles, and exterior differential operations. Although this choice introduces more abstract mathematical structures, it brings three irreplaceable advantages: (1) All physical quantities (strains, stresses) are interpreted as differential forms, and their integral and differential essence (e.g., Stokes’ theorem) is precisely retained; (2) The antisymmetry of connection 1-forms minimizes independent variables and automatically degenerates to classical results in the thin-shell limit; (3) The theory is naturally compatible with Finite Element Exterior Calculus (FEEC), providing a direct path for constructing locking-free, structure-preserving $C^0$ continuous shell elements. Therefore, this paper can be viewed as a geometric explicitation and deep structural supplement to the Simo-Fox classic framework.

Combining fiber bundle theory, exterior differential forms, and geometrically exact kinematics is the natural trend of current geometric mechanics development. A shell is essentially a fiber bundle structure where a three-dimensional manifold reduces to a two-dimensional manifold, with the thickness direction being the fiber. Viewing the dimensional reduction process as integration along the fiber, realized via generalized Stokes’ theorem and exterior covariant derivatives, maximizes the maintenance of the theory’s coordinate-independence and automatically reveals the geometric origin of curvature coupling terms [13,16].

1.3. Motivation and Contributions of This Paper

Although the work of Marsden-Hughes and Simo-Fox has laid a solid foundation, literature that starts from fiber bundles and differential forms to derive a complete geometrically exact Kirchhoff-Love shell theory and elaborate on its nonlinear dynamic coupling mechanisms in detail remains scarce and scattered. This paper aims to systematically construct this theoretical framework. Our core objective is to extend kinematics to finite deformation (geometric nonlinearity) to accurately capture large rotation effects, while retaining linear constitutive relations (Saint-Venant-Kirchhoff model) at the material level. This combination of "geometric nonlinearity + material linearity" explicitly targets "large rotation, small strain" conditions and is extremely applicable to the engineering analysis of most hard materials and flexible electronic devices.

The specific contributions of this paper are as follows:

• Fiber Bundle Geometry of the Reference Configuration: The reference shell is modeled as a fiber bundle with the mid-surface as the base. Through Cartan’s first structure equation, antisymmetric connection 1-forms are derived, encoding all geometric information in the lines-of-curvature coordinate system with minimal algebraic cost, and the Gauss-Codazzi equations are naturally recovered through the vanishing of the curvature 2-form.
• Nonlinear Kinematics Based on Pull-Back Metric: Abandoning the linearized Lie derivative, we introduce a nonlinear deformation mapping from the reference manifold to the current spatial configuration. Using the pull-back metric ${\mathsf{\Phi }}^{*}\mathit{\delta }$, we define the Green-Lagrange strain tensor and precisely separate nonlinear membrane and bending strains. Since the pull-back metric is used, the theory strictly produces zero strain under pure large rigid body rotation, eliminating the "spurious strains" in linear theory.
• Virtual Work Principle and Exterior Covariant Dynamics: In the Lagrangian description, we employ the second Piola-Kirchhoff (PK2) stress conjugate to the Green-Lagrange strain energy and formulate it as a vector-valued stress 2-form. The PK2 stress is designed for finite deformation and remains directionally invariant under rigid body rotation, ensuring the objectivity of the theory. Through the principle of virtual work, local equilibrium is transformed into a global variational statement of internal virtual work equaling external virtual work.
• Integration Along Fibers and Nonlinear Dimensional Reduction: We demonstrate in detail how to reduce the three-dimensional virtual work principle to the two-dimensional mid-surface through integration along the normal fiber. Nonlinear geometric coupling terms (e.g., the coupling of normal rotation and membrane forces, i.e., the geometric stiffness effect) naturally emerge through the variation of virtual strains and the use of exterior differential properties, without cumbersome coordinate expansion.
• Linear Constitutive Closure: The linear constitutive relation between stress and strain is retained. Since it is defined on the reference configuration, this "linear" relation still satisfies frame-invariance, making the entire theory suitable for large rotation problems while maintaining geometric nonlinearity.

2. Geometric Foundations of the Reference Configuration

To establish a geometrically exact shell theory, we must first lay a solid geometric foundation in the reference configuration. We will not rely on specific coordinate system choices but view the shell as a differentiable manifold equipped with specific metrics and connections.

2.1. Fiber Bundle Structure of the Shell Manifold

Let $\mathcal{S}$ be the undeformed reference configuration of the elastic shell, a three-dimensional smooth manifold with boundary. We endow $\mathcal{S}$ with a locally trivial fiber bundle structure, denoted as $(\mathcal{S},\mathcal{M},\pi )$, where $\mathcal{M}$ is a two-dimensional smooth, orientable Riemannian manifold representing the mid-surface of the shell; $\pi :\mathcal{S}\to \mathcal{M}$ is the bundle projection. For any point $p\in \mathcal{M}$, its preimage ${\pi }^{-1}\left(p\right)$ is a one-dimensional fiber representing the normal thickness interval $I=[-h/2,h/2]$ at point p, where $h\ll L$ (L is the characteristic length of the mid-surface). The shell element is illustrated in Figure 1.

The tangent bundle $T\mathcal{S}$ naturally splits into a direct sum of horizontal and vertical sub-bundles: $T\mathcal{S}=H\mathcal{S}\oplus V\mathcal{S}$. The horizontal sub-bundle $H\mathcal{S}$ is isomorphic to the pull-back bundle ${\pi }^{*}\left(T\mathcal{M}\right)$, containing vectors tangent to surfaces of constant thickness; the vertical sub-bundle $V\mathcal{S}$ is spanned by the normal vector field $\partial /\partial z$.

To parameterize the manifold, we introduce a local coordinate chart $({\alpha }^1,{\alpha }^2)$ on $\mathcal{M}$. By the fundamental theorem of surfaces, as long as the principal curvatures are unequal, we can always choose lines of curvature coordinates [12]. In this case, the coordinate lines coincide with the principal curvature directions. Let $R_1({\alpha }^1,{\alpha }^2)$ and $R_2({\alpha }^1,{\alpha }^2)$ be the principal radii of curvature, and the third coordinate z be taken along the unit normal vector $\mathit{N}$ of the mid-surface. Thus, the local coordinates on $\mathcal{S}$ are $({\alpha }^1,{\alpha }^2,z)$.

2.2. Reference Metric and Orthogonal Moving Frame

The geometric properties of the reference manifold $\mathcal{S}$ are determined by its Riemannian metric $\mathit{G}$. Consider a point on the manifold at a distance z from the mid-surface along the normal. In the principal direction a, the distance from this point to the corresponding center of curvature is $R_a+z$. The infinitesimal arc length is $\mathrm{d}{s}^a=(R_a+z)\mathrm{d}{\phi }^a$, where $\mathrm{d}{\phi }^a=A_a\mathrm{d}{\alpha }^a/R_a$ is the infinitesimal angle subtended at the center of curvature, and $A_a$ is the Lamé parameter of the mid-surface. Therefore, the reference metric tensor $\mathit{G}$ is:

$$\mathit{G}={\displaystyle \sum _{a=1}^2}{A}_a^2{\left(1+{\displaystyle \frac{z}{R_a}}\right)}^2\mathrm{d}{\alpha }^a\otimes \mathrm{d}{\alpha }^a+\mathrm{d}z\otimes \mathrm{d}z.$$

(1)

This metric depends explicitly on z, reflecting the physical fact that parallel surfaces have different intrinsic metrics. Traditional methods would proceed to calculate the inverse metric and 18 Christoffel symbols ${\mathsf{\Gamma }}_{jk}^i$. To avoid this coordinate disaster, we adopt Cartan’s method of moving frames [8,17].

Instead of using the coordinate co-frame $\{\mathrm{d}{\alpha }^1,\mathrm{d}{\alpha }^2,\mathrm{d}z\}$, we construct an orthogonal co-frame $\{{\mathsf{\Theta }}^1,{\mathsf{\Theta }}^2,{\mathsf{\Theta }}^3\}$ by observing equation (1):

$${\mathsf{\Theta }}^a=A_a\left(1+{\displaystyle \frac{z}{R_a}}\right)\mathrm{d}{\alpha }^a,{\mathsf{\Theta }}^3=\mathrm{d}z.$$

(2)

In this basis, the metric takes the standard Euclidean form:

$$\mathit{G}={\delta }_{IJ}{\mathsf{\Theta }}^I\otimes {\mathsf{\Theta }}^J,I,J\in \{1,2,3\}.$$

(3)

The dual moving frame (tangent basis) is ${\mathit{E}}_a={\displaystyle \frac{1}{A_a(1+z/R_a)}}{\displaystyle \frac{\partial }{\partial {\alpha }^a}}$ and ${\mathit{E}}_3={\displaystyle \frac{\partial }{\partial z}}$. The reference volume form is:

$$\mathrm{d}{V}_0={\mathsf{\Theta }}^1\wedge {\mathsf{\Theta }}^2\wedge {\mathsf{\Theta }}^3=A_1{A}_2{\mu }_0\left(z\right)\mathrm{d}{\alpha }^1\wedge \mathrm{d}{\alpha }^2\wedge \mathrm{d}z,$$

(4)

where ${\mu }_0\left(z\right)=\left(1+{\displaystyle \frac{z}{R_1}}\right)\left(1+{\displaystyle \frac{z}{R_2}}\right)$ is the area change factor due to normal displacement. The Hodge star operator ${★}_0$ acts very simply in this orthogonal basis; for example, ${★}_0{\mathsf{\Theta }}_1={\mathsf{\Theta }}^2\wedge {\mathsf{\Theta }}^3$.

2.3. Connection 1-Forms and Cartan’s First Structure Equation

The connectivity of the manifold is encoded in the connection 1-forms ${\omega }_J^I\in {\mathsf{\Omega }}^1\left(\mathcal{S}\right)$. Since we have chosen an orthogonal frame, the connection is metric-compatible ($\nabla \mathit{G}=\mathbf{0}$), meaning the connection matrix is antisymmetric: ${\omega }_{IJ}=-{\omega }_{JI}$. This property directly reduces the number of independent connection components from 18 to 3. Assuming a torsion-free connection (the hypothesis of no dislocations in standard continuum mechanics), Cartan’s first structure equation uniquely determines the connection:

$$\mathrm{d}{\mathsf{\Theta }}^I+{\omega }_J^I\wedge {\mathsf{\Theta }}^J=0.$$

(5)

To solve for ${\omega }_J^I$, we apply the exterior derivative to ${\mathsf{\Theta }}^a$:

$$\begin{array}{cc}\hfill \mathrm{d}{\mathsf{\Theta }}^a& =\mathrm{d}\left(A_a\left(1+{\displaystyle \frac{z}{R_a}}\right)\mathrm{d}{\alpha }^a\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{\partial }{\partial {\alpha }^b}}\left(A_a\left(1+{\displaystyle \frac{z}{R_a}}\right)\right)\mathrm{d}{\alpha }^b\wedge \mathrm{d}{\alpha }^a+{\displaystyle \frac{\partial }{\partial z}}\left(A_a\left(1+{\displaystyle \frac{z}{R_a}}\right)\right)\mathrm{d}z\wedge \mathrm{d}{\alpha }^a.\hfill \end{array}$$

(7)

Through algebraic simplification and replacing the coordinate basis back to the orthogonal basis ${\mathsf{\Theta }}^a$, after careful index matching, we can extract the non-zero connection 1-forms from (5). The components related to normal curvature are:

$${\omega }_3^a={\displaystyle \frac{{\mathsf{\Theta }}^a}{R_a+z}},a\in \{1,2\}.$$

(8)

This represents the bending of the normal vector towards the center of curvature as one moves along the surface. The connection 1-form for in-plane rotation is:

$${\omega }_2^1=-{\omega }_1^2={\displaystyle \frac{{\left(A_1(1+z/R_1)\right)}_{,2}{\mathsf{\Theta }}^1-{\left(A_2(1+z/R_2)\right)}_{,1}{\mathsf{\Theta }}^2}{A_1{A}_2(1+z/R_1)(1+z/R_2)}}.$$

(9)

These two simple 1-forms completely replace the lengthy Christoffel symbol expressions in classical tensor calculus. The covariant derivatives of vector fields and differential forms are now simplified to the exterior derivative plus the wedge product with these sparse 1-forms.

2.4. Curvature Constraints and Mid-Surface Geometry

Since the reference configuration $\mathcal{S}$ is isometrically embedded in flat Euclidean space ${\mathbb{R}}^3$, its intrinsic Riemannian curvature must be zero. According to Cartan’s second structure equation:

$${\mathsf{\Omega }}_J^I=\mathrm{d}{\omega }_J^I+{\omega }_K^I\wedge {\omega }_J^K=0.$$

(10)

Substituting equations (8) and (9) into this identity imposes strict differential constraints, namely the famous Gauss-Codazzi-Mainardi equations. For example, calculating ${\mathsf{\Omega }}_3^1=0$ and projecting to the mid-surface $z=0$ yields the Codazzi equation ${\displaystyle \frac{A_{1,2}}{R_2}}={\left({\displaystyle \frac{A_1}{R_1}}\right)}_{,2}$. Calculating ${\mathsf{\Omega }}_2^1=0$ yields the Gauss equation, where the Gaussian curvature $K=1/\left(R_1{R}_2\right)$ is completely determined by the intrinsic metric. In our geometric framework, these are not added assumptions but necessary consequences of the flatness of the ambient space.

The geometry of the mid-surface $\mathcal{M}$ is induced by the pull-back of the inclusion mapping $\iota :\mathcal{M}\hookrightarrow \mathcal{S}$ at $z=0$. The mid-surface co-frame is ${\varphi }^a={\iota }^{*}{\mathsf{\Theta }}^a=A_a\mathrm{d}{\alpha }^a$. The first fundamental form is $\mathit{a}={\delta }_{ab}{\varphi }^a\otimes {\varphi }^b$, and the second fundamental form is $\mathit{d}={\displaystyle \frac{1}{R_1}}{\varphi }^1\otimes {\varphi }^1+{\displaystyle \frac{1}{R_2}}{\varphi }^2\otimes {\varphi }^2$. The intrinsic connection on the mid-surface is ${\tilde{\omega }}_2^1={\iota }^{*}{\omega }_2^1$, which completely determines the Levi-Civita covariant derivative $\tilde{\nabla }$ on the mid-surface. Since the frame is orthogonal, the connection coefficients are antisymmetric in the first two indices, which is the key to avoiding massive algebraic cancellation in the moving frame method.

3. Finite Deformation Kinematics and Pull-Back Metric

Classical shell theory usually defines strain as the Lie derivative of the displacement vector field $\mathit{U}$ with respect to the reference metric under the small deformation assumption. However, when the structure undergoes large displacements or rotations, this linearization generates spurious strains. To establish a geometrically exact theory, we must abandon the linear expansion of the displacement field and instead directly operate on the global nonlinear deformation mapping. This kinematic nonlinearity is the core distinction between this theory and classical linear theory.

3.1. Nonlinear Deformation Mapping and Kirchhoff-Love Constraint

Let the deformation mapping from the reference configuration $\mathcal{S}$ to the current spatial configuration ${\mathbb{R}}^3$ be $\mathbf{\Phi }:\mathcal{S}\to {\mathbb{R}}^3$. The Kirchhoff-Love hypothesis requires that normals to the mid-surface before deformation remain normal and unstretched after deformation. This means the deformation of the shell is completely determined by the mapping of the mid-surface $\mathit{\phi }:\mathcal{M}\to {\mathbb{R}}^3$ and the unit normal vector field of the current mid-surface $\mathit{n}:\mathcal{M}\to S^2$:

$$\mathbf{\Phi }({\alpha }^1,{\alpha }^2,z)=\mathit{\phi }({\alpha }^1,{\alpha }^2)+z\mathit{n}({\alpha }^1,{\alpha }^2).$$

(11)

Here, the normal vector $\mathit{n}$ is no longer a small rotation approximation of the reference normal $\mathit{N}$, but an independent kinematic variable on the tangent bundle, taking values in the unit sphere $S^2$ fiber bundle. Since $\mathit{n}$ must be perpendicular to the deformed mid-surface, it must satisfy the orthogonality constraint:

$$\mathit{n}·{\partial }_a\mathit{\phi }=0,a\in \{1,2\}.$$

(12)

The set of kinematic variables $(\mathit{\phi },\mathit{n})$ can accurately describe arbitrarily large rigid body translations and rotations. Under pure rigid body rotation $\mathit{Q}\in SO\left(3\right)$, $\mathit{\phi }\to \mathit{Q}\mathit{\phi }$, $\mathit{n}\to \mathit{Q}\mathit{n}$, and the form of the mapping $\mathbf{\Phi }$ remains invariant, ensuring the objectivity of the theory.

3.2. Pull-Back Metric and Green-Lagrange Strain

In the geometrically nonlinear framework, the strain measure cannot be calculated directly through the spatial metric on the current configuration; it must be pulled back to the reference configuration to eliminate the influence of rigid body motion. Let $\mathit{\delta }$ be the Euclidean metric of ${\mathbb{R}}^3$. The pull-back metric $\mathit{C}={\mathbf{\Phi }}^{*}\mathit{\delta }$ is a second-order covariant tensor defined on the reference manifold $\mathcal{S}$, physically known as the right Cauchy-Green deformation tensor.

The Green-Lagrange strain tensor $\mathit{E}$ is defined as half the difference between the pull-back metric and the reference metric:

$$\mathit{E}={\displaystyle \frac{1}{2}}({\mathbf{\Phi }}^{*}\mathit{\delta }-\mathit{G})={\displaystyle \frac{1}{2}}(\mathit{C}-\mathit{G}).$$

(13)

Since both $\mathit{C}$ and $\mathit{G}$ are defined on the reference configuration, $\mathit{E}$ is also a tensor on the reference configuration. Under pure rigid body motion, $\mathit{C}=\mathit{G}$, and thus $\mathit{E}\equiv \mathbf{0}$, which excludes spurious strains by definition.

Remark (Note on Not Explicitly Using the Deformation Gradient Tensor). In classical continuum mechanics, the deformation gradient tensor $\mathit{F}$ is usually a core physical quantity. However, in the geometrically exact framework of this paper, we do not explicitly use $\mathit{F}$ for three reasons: First, the right Cauchy-Green deformation tensor $\mathit{C}$ is essentially the inner product of the deformation gradient ($\mathit{C}={\mathit{F}}^T·\mathit{F}$). By directly introducing the pull-back metric $\mathit{C}={\mathbf{\Phi }}^{*}\mathit{\delta }$ and differentiating the deformation mapping when calculating its components $C_{ab}={\partial }_a\mathbf{\Phi }·{\partial }_b\mathbf{\Phi }$, we have implicitly and completely retained all information of the deformation gradient. Second, the direct expansion based on the Kirchhoff-Love kinematic variables $(\mathit{\phi },\mathit{n})$ naturally derives membrane and bending strains, skipping the step of explicitly constructing the $3\times 3$ deformation gradient matrix, which greatly simplifies algebraic operations. Third, the component representation of the deformation gradient tensor depends highly on the choice of coordinate basis. This paper aims to build a coordinate-independent theory using moving frames and differential forms. Using the pull-back metric to replace the deformation gradient matrix maximizes the preservation of the geometric covariance and objectivity of the theory, reflecting the core philosophy of modern geometric mechanics (e.g., the Marsden-Hughes framework).

3.3. Thickness Expansion and Nonlinear Strain Separation

To reduce the three-dimensional strain, we need to calculate the components of the pull-back metric and expand them along the thickness z. From (11), the tangential components of the pull-back metric are:

$$\begin{array}{cc}\hfill C_{ab}={\displaystyle \frac{\partial \mathbf{\Phi }}{\partial {\alpha }^a}}·{\displaystyle \frac{\partial \mathbf{\Phi }}{\partial {\alpha }^b}}& =({\partial }_a\mathit{\phi }+z{\partial }_a\mathit{n})·({\partial }_b\mathit{\phi }+z{\partial }_b\mathit{n})\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & ={\partial }_a\mathit{\phi }·{\partial }_b\mathit{\phi }+z({\partial }_a\mathit{\phi }·{\partial }_b\mathit{n}+{\partial }_b\mathit{\phi }·{\partial }_a\mathit{n})+z^2{\partial }_a\mathit{n}·{\partial }_b\mathit{n}.\hfill \end{array}$$

(15)

We define the first fundamental form ${\mathit{a}}_t$ and the second fundamental form ${\mathit{d}}_t$ of the current mid-surface:

$$\begin{array}{cc}\hfill {\mathit{a}}_t& ={\delta }_{ab}({\partial }_c\mathit{\phi }·{\partial }_d\mathit{\phi }){\varphi }^c\otimes {\varphi }^d,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\mathit{d}}_t& ={\delta }_{ab}({\partial }_c\mathit{\phi }·{\partial }_d\mathit{n}){\varphi }^c\otimes {\varphi }^d=-{\delta }_{ab}({\partial }_{cd}\mathit{\phi }·\mathit{n}){\varphi }^c\otimes {\varphi }^d.\hfill \end{array}$$

(17)

Neglecting higher-order terms in $z^2$ (the standard practice in thin shell approximation, corresponding to ignoring normal stress effects in the thickness direction), the pull-back metric can be compactly represented as:

$${\mathbf{\Phi }}^{*}\mathit{\delta }={\mathit{a}}_t+2z{\mathit{d}}_t+\mathrm{d}z\otimes \mathrm{d}z.$$

(18)

Simultaneously, the expansion of the reference metric is $\mathit{G}=\mathit{a}+2z\mathit{d}+\mathrm{d}z\otimes \mathrm{d}z$. Substituting into the strain definition (13), we obtain:

$$\mathit{E}={\displaystyle \frac{1}{2}}({\mathit{a}}_t-\mathit{a})+z({\mathit{d}}_t-\mathit{d}).$$

(19)

This naturally separates the Green-Lagrange strain into a zero-order membrane strain ${\mathit{E}}^{\left(0\right)}$ and a first-order bending strain ${\mathit{K}}^{\left(1\right)}$:

$$\begin{array}{cc}\hfill {\mathit{E}}^{\left(0\right)}& ={\displaystyle \frac{1}{2}}({\mathit{a}}_t-\mathit{a}),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\mathit{K}}^{\left(1\right)}& ={\mathit{d}}_t-\mathit{d}.\hfill \end{array}$$

(21)

This derivation is purely geometric and nonlinear. ${\mathit{E}}^{\left(0\right)}$ precisely measures the stretching and shearing of the mid-surface due to deformation, while ${\mathit{K}}^{\left(1\right)}$ precisely measures the true change in curvature. Even if the mid-surface is not stretched (${\mathit{a}}_t=\mathit{a}$), non-zero bending strain will arise as long as the curvature changes (${\mathit{d}}_t\ne \mathit{d}$).

4. Dynamics Based on the Principle of Virtual Work

Since the kinematics are completely nonlinear, traditional direct force balance derivations become extremely complex due to the expansion of nonlinear terms. In the Lagrangian description, establishing dynamic equations through the principle of virtual work is the most rigorous approach that naturally preserves objectivity.

4.1. Stress 2-Forms on the Reference Configuration

To be energy-conjugate to the Green-Lagrange strain $\mathit{E}$, we adopt the second Piola-Kirchhoff (PK2) stress tensor $\mathit{S}$. Unlike the Cauchy stress, PK2 stress is defined on the reference configuration and does not change direction with rigid body rotation. In finite deformation mechanics, PK2 stress is a stress measure specifically designed for the Lagrangian description; it ensures that the stress-strain relationship on the reference configuration remains objective regardless of how large the rotation the structure undergoes. This makes it extremely convenient for handling large rotation problems.

In the language of modern differential geometry, we do not view stress as a second-order tensor but as a vector-valued stress 2-form $\Xi \in {\mathsf{\Omega }}^2(\mathcal{S},T^{*}\mathcal{S})$ defined on the reference manifold $\mathcal{S}$:

$${\Xi }^I=S^{IJ}{\iota }_{{\mathit{E}}_J}\mathrm{d}{V}_0=S^{IJ}{★}_0{\mathsf{\Theta }}_J.$$

(22)

Here, ${★}_0{\mathsf{\Theta }}_J$ is a 2-form representing the directed area element orthogonal to the basis vector ${\mathit{E}}_J$. The stress 2-form $\Xi$ intrinsically represents the flux of force through that area element in the ${\mathit{E}}_I$ direction. The Hodge star operator automatically handles the orientation of the surface, avoiding the complex normal vector cross-products in tensor index methods.

4.2. Principle of Virtual Work

The equilibrium of the system is governed by the principle of virtual work: for any kinematically admissible virtual displacement $\delta \mathbf{\Phi }$, the internal virtual work equals the external virtual work plus the inertial virtual work. On the reference configuration, the three-dimensional principle of virtual work is expressed as:

$${\int }_{\mathcal{S}}{S}^{IJ}\delta E_{IJ}\mathrm{d}{V}_0={\int }_{\mathcal{S}}{\rho }_0\mathit{B}·\delta \mathbf{\Phi }\mathrm{d}{V}_0+{\int }_{\partial {\mathcal{S}}_t}\overline{\mathit{T}}·\delta \mathbf{\Phi }\mathrm{d}{A}_0-{\int }_{\mathcal{S}}{\rho }_0\ddot{\mathbf{\Phi }}·\delta \mathbf{\Phi }\mathrm{d}{V}_0,$$

(23)

where $\delta E_{IJ}$ is the variation of the Green-Lagrange strain, $\mathit{B}$ is the body force per unit mass, $\overline{\mathit{T}}$ is the prescribed surface traction on the Neumann boundary $\partial {\mathcal{S}}_t$, ${\rho }_0$ is the reference density, and $\ddot{\mathbf{\Phi }}$ is the acceleration. Equation (23) is completely coordinate-independent and implicitly implies conservation of linear and angular momentum.

4.3. Variation of Virtual Strain and Nonlinear Geometric Coupling

To reduce the equation dimension, we need to calculate the virtual strain $\delta \mathit{E}$. Since $\mathit{E}={\displaystyle \frac{1}{2}}({\mathbf{\Phi }}^{*}\mathit{\delta }-\mathit{G})$ and the reference metric $\mathit{G}$ is invariant, the virtual strain is the variation of the pull-back metric:

$$\delta \mathit{E}={\displaystyle \frac{1}{2}}\delta \left({\mathbf{\Phi }}^{*}\mathit{\delta }\right).$$

(24)

The variation of the deformation mapping is $\delta \mathbf{\Phi }=\delta \mathit{\phi }+z\delta \mathit{n}$. Through algebraic operations, it can be shown that the virtual strain can also be expanded in thickness:

$$\delta \mathit{E}=\delta {\mathit{E}}^{\left(0\right)}+z\delta {\mathit{K}}^{\left(1\right)}.$$

(25)

The explicit expressions for membrane virtual strain and bending virtual strain are:

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

(26)

$$\begin{array}{cc}\hfill \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}).\hfill \end{array}$$

(27)

The most critical geometric feature here is that the virtual strain depends not only on the mid-surface variation $\delta \mathit{\phi }$ but also highly on the variation of the normal vector field $\delta \mathit{n}$. Due to the orthogonality constraint (12), $\delta \mathit{n}$ must be perpendicular to $\mathit{n}$, i.e., $\delta \mathit{n}·\mathit{n}=0$. Thus, $\delta \mathit{n}$ is actually an infinitesimal rotation vector. This intertwining of $\delta \mathit{\phi }$ and $\delta \mathit{n}$ is the core distinction between geometrically nonlinear theory and linear theory. In the virtual work integral, $\delta {\mathit{K}}^{\left(1\right)}$ will generate coupling terms between membrane forces and normal rotation, which cannot be captured by classical linear theory.

5. Integration Along Fibers and Dimensional Reduction Dynamics

The core of shell theory lies in dimensional reduction. We integrate the three-dimensional virtual work equation (23) along the normal fiber z.

5.1. Reduction of Internal Virtual Work

Substituting $\delta \mathit{E}=\delta {\mathit{E}}^{\left(0\right)}+z\delta {\mathit{K}}^{\left(1\right)}$ into the internal virtual work integral and using the reference volume form $\mathrm{d}{V}_0={\mu }_0\left(z\right)\mathrm{d}{A}_0\mathrm{d}z$, we have:

$${\int }_{\mathcal{S}}{S}^{ab}\delta E_{ab}\mathrm{d}{V}_0={\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}_0.$$

(28)

We define the membrane force tensor $\mathit{N}$ and moment tensor $\mathit{M}$ based on the PK2 stress $\mathit{S}$:

$$\begin{array}{cc}\hfill N^{ab}& ={\int }_{-h/2}^{h/2}{S}^{ab}{\mu }_0\left(z\right)\mathrm{d}z,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill M^{ab}& ={\int }_{-h/2}^{h/2}z{S}^{ab}{\mu }_0\left(z\right)\mathrm{d}z.\hfill \end{array}$$

(30)

Since $\delta {\mathit{E}}^{\left(0\right)}$ and $\delta {\mathit{K}}^{\left(1\right)}$ depend only on mid-surface coordinates (independent of z), the integral can be separated directly. The internal virtual work simplifies to:

$$\delta W_{int}={\int }_{\mathcal{M}}\left(N^{ab}\delta E_{ab}^{\left(0\right)}+M^{ab}\delta K_{ab}^{\left(1\right)}\right)\mathrm{d}{A}_0.$$

(31)

This is the internal part of the nonlinear principle of virtual work on the two-dimensional mid-surface.

5.2. Reduction of External and Inertial Forces

Similarly, integrating the body force $\mathit{B}$ along the thickness yields the equivalent area force $\mathit{p}$ and equivalent area moment $\mathit{m}$ on the mid-surface:

$$\begin{array}{cc}\hfill \mathit{p}& ={\int }_{-h/2}^{h/2}{\rho }_0\mathit{B}{\mu }_0\left(z\right)\mathrm{d}z+\mathrm{boundary}\mathrm{term}\mathrm{integration},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \mathit{m}& ={\int }_{-h/2}^{h/2}z{\rho }_0\mathit{B}{\mu }_0\left(z\right)\mathrm{d}z.\hfill \end{array}$$

(33)

The reduction of the acceleration term requires utilizing the kinematics $\mathbf{\Phi }=\mathit{\phi }+z\mathit{n}$, thus $\ddot{\mathbf{\Phi }}=\ddot{\mathit{\phi }}+z\ddot{\mathit{n}}$. Neglecting higher-order small quantities of rotational inertia (order $h^3$), the inertial virtual work reduces to:

$$\delta W_{iner}={\int }_{\mathcal{M}}{\rho }_{0}h\ddot{\mathit{\phi }}·\delta \mathit{\phi }\mathrm{d}{A}_0.$$

(34)

5.3. Nonlinear Principle of Virtual Work on the Mid-Surface

Synthesizing the above results, the geometrically exact nonlinear principle of virtual work on the reference mid-surface $\mathcal{M}$ is stated as: find $(\mathit{\phi },\mathit{n})$ such that for all kinematically admissible variations $(\delta \mathit{\phi },\delta \mathit{n})$, the following holds:

$${\int }_{\mathcal{M}}\left(N^{ab}\delta E_{ab}^{\left(0\right)}+M^{ab}\delta K_{ab}^{\left(1\right)}\right)\mathrm{d}{A}_0={\int }_{\mathcal{M}}\mathit{p}·\delta \mathit{\phi }\mathrm{d}{A}_0+{\int }_{\partial {\mathcal{M}}_t}\overline{\mathit{n}}·\delta \mathit{\phi }\mathrm{d}s-{\int }_{\mathcal{M}}{\rho }_{0}h\ddot{\mathit{\phi }}·\delta \mathit{\phi }\mathrm{d}{A}_0.$$

(35)

5.4. Emergence of Strong Form and Nonlinear Equilibrium Equations

By expanding $\delta {\mathit{E}}^{\left(0\right)}$ and $\delta {\mathit{K}}^{\left(1\right)}$ in equation (35), and utilizing the covariant derivative $\tilde{\nabla }$ on the mid-surface and partial integration (i.e., using Stokes’ theorem $\int \mathrm{d}(\dots )=0$ in exterior calculus), we can extract the local strong form equilibrium equations. Since the variations $\delta \mathit{\phi }$ and $\delta \mathit{n}$ are arbitrary within the domain (except for the constraint $\delta \mathit{n}·\mathit{n}=0$), we obtain: Linear Momentum Balance:

$${\tilde{\nabla }}_b{N}^{ab}+{\overline{n}}^a={\rho }_{0}h{\ddot{\phi }}^a,$$

(36)

Angular Momentum Balance:

$${\tilde{\nabla }}_b{M}^{ab}+Q^a=0,$$

(37)

Normal Balance (Core of Nonlinear Coupling):

$${\tilde{\nabla }}_a{Q}^a+N^{ab}(\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi })+p^n={\rho }_{0}h{\ddot{\phi }}^n.$$

(38)

Here, $Q^a$ is the transverse shear force. The key lies in the nonlinear term $N^{ab}(\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi })$ in the third equation. Since $\mathit{n}·{\partial }_a{\partial }_b\mathit{\phi }$ is precisely the second fundamental form $d_{t,ab}$ of the current mid-surface, this term represents the coupling of membrane forces with the current curvature (i.e., geometric stiffness). It is this nonlinear term that allows the structure to exhibit real nonlinear behaviors such as "stress stiffening" or "stress softening" under large rotations. In our geometric framework, this complex nonlinear term naturally emerges from the $\delta \mathit{n}$ term contained in the virtual strain variation $\delta {\mathit{K}}^{\left(1\right)}$, without relying on specific coordinate expansion, completely preserving tensor covariance.

6. Linear Constitutive Relations and System Closure

To close the system, we need to establish constitutive relations between stress and strain. As mentioned earlier, this paper focuses on engineering problems of large rotation and small strain, so we retain the linear elastic hypothesis of the material.

6.1. Objectivity of Saint-Venant-Kirchhoff Material

The Saint-Venant-Kirchhoff constitutive model assumes that the PK2 stress $\mathit{S}$ is proportional to the Green-Lagrange strain $\mathit{E}$. Since both quantities are defined on the reference configuration, they are invariant under pure rigid body rotation. Therefore, this "linear" relationship is actually objective and will not produce non-physical stresses due to large rotations. This allows us to safely use linear constitutive relations without undermining the overall geometric nonlinearity of the theory.

For isotropic materials, on the tangential plane of the reference mid-surface, the constitutive relation is:

$$S^{ab}=H^{abcd}{E}_{cd},$$

(39)

where $H^{abcd}$ is the fourth-order isotropic elasticity tensor on the reference mid-surface:

$$H^{abcd}={\displaystyle \frac{E}{1-{\nu }^2}}\left[{\displaystyle \frac{1-\nu }{2}}({\delta }^{ac}{\delta }^{bd}+{\delta }^{ad}{\delta }^{bc})+\nu {\delta }^{ab}{\delta }^{cd}\right],$$

(40)

where E is Young’s modulus and $\nu$ is Poisson’s ratio.

6.2. Reduced Constitutive Equations

Substituting the strain expansion $\mathit{E}={\mathit{E}}^{\left(0\right)}+z{\mathit{K}}^{\left(1\right)}$ into the constitutive equation (39) and combining the definitions of membrane force and moment, integration along the thickness yields:

$$\begin{array}{cc}\hfill N^{ab}& ={\int }_{-h/2}^{h/2}{H}^{abcd}(E_{cd}^{\left(0\right)}+z{K}_{cd}^{\left(1\right)}){\mu }_0\left(z\right)\mathrm{d}z,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill M^{ab}& ={\int }_{-h/2}^{h/2}z{H}^{abcd}(E_{cd}^{\left(0\right)}+z{K}_{cd}^{\left(1\right)}){\mu }_0\left(z\right)\mathrm{d}z.\hfill \end{array}$$

(42)

For thin shells, we approximate ${\mu }_0\left(z\right)\approx 1$. Since ${\mathit{E}}^{\left(0\right)}$ and ${\mathit{K}}^{\left(1\right)}$ are independent of z, and the integral $\int z\mathrm{d}z=0$, membrane and bending terms decouple. Defining the extensional stiffness $C={\displaystyle \frac{Eh}{1-{\nu }^2}}$ and bending stiffness $B={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}$, we obtain the two-dimensional constitutive equations:

$$\begin{array}{cc}\hfill N^{ab}& =C{H}^{abcd}{E}_{cd}^{\left(0\right)},\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill M^{ab}& =B{H}^{abcd}{K}_{cd}^{\left(1\right)}.\hfill \end{array}$$

(44)

Although the constitutive relation is linear in form, since ${\mathit{E}}^{\left(0\right)}$ and ${\mathit{K}}^{\left(1\right)}$ contain the nonlinear deformation mapping $\mathit{\phi }$ and normal $\mathit{n}$, the final relationship between force and deformation is highly nonlinear. These equations, together with the preceding nonlinear equilibrium equations and the kinematic constraint (12), constitute a closed geometrically exact Kirchhoff-Love shell theory.

7. Conclusions and Outlook

This paper systematically constructs a geometrically exact elastic shell dynamic theory based on fiber bundles and differential forms. By viewing the shell as a fiber bundle on the reference manifold, we utilize Cartan’s method of moving frames to greatly simplify the geometric description, replacing lengthy Christoffel symbols with antisymmetric connection 1-forms. In kinematics, we introduce a nonlinear deformation mapping and define the Green-Lagrange strain through the pull-back metric, precisely separating nonlinear membrane and bending strains, completely eliminating spurious strains under large rotations. In dynamics, by adopting a vector-valued stress 2-form based on PK2 stress, we take the principle of virtual work as the core and naturally reduce it to the mid-surface through integration along the normal fiber. Complex geometric nonlinear coupling terms (such as the interaction of membrane forces with current curvature) automatically emerge during the variational process and fully preserve coordinate-independence. Finally, linear Saint-Venant-Kirchhoff constitutive relations are retained, making the theory explicitly applicable to engineering structures with large rotations and small strains.

It is worth emphasizing that the geometric framework constructed in this paper is completely equivalent in physical substance to the single-director model of Simo-Fox (1989). When we degenerate the Cartan co-frame in this paper to a coordinate basis and explicitly substitute the connection 1-forms ${\omega }_J^I$ into the strain expressions, the membrane, bending, and shear strain components in classical literature can be accurately recovered. The difference between the two lies only in the language of formulation: the classical framework "hides" geometry within vector derivatives, while this paper makes it "explicit" as algebraic operations of differential forms. Although this explicitness adds to the initial abstraction, it provides a unified mathematical platform for subsequent structure-preserving discretization, multi-field coupling (e.g., piezoelectric, thermo-mechanical), and non-Euclidean geometry extensions. Therefore, this work is not a negation of the classics but a deep excavation and formal enhancement of their geometric content.

The theoretical framework of this paper is entirely based on exterior differentiation and differential forms, which gives it a natural affinity with Discrete Exterior Calculus (DEC) or Finite Element Exterior Calculus (FEEC) [27,28]. By directly mapping the continuous theory to piecewise polynomial differential form spaces, efficient numerical algorithms that preserve geometric structures and avoid locking (especially shear and membrane locking) can be constructed. Existing work [23,24,25] has successfully implemented discretization based on mixed variational principles and FEEC, demonstrating excellent numerical performance. Therefore, the theory in this paper lays a solid continuous foundation for developing the next generation of geometrically exact, structure-preserving shell elements.

This geometric framework is not only an elegant rewrite of classical nonlinear shell theory but also opens broad avenues for future research. First, since the entire derivation is based on exterior differentiation, the theory can be directly mapped into the Discrete Exterior Calculus (DEC) framework [22], thereby developing structure-preserving finite element algorithms that naturally avoid shear locking and membrane locking. Second, by modifying the reference metric $\mathit{G}$ to have non-zero Riemann curvature, the theory can be seamlessly extended to non-Euclidean elasticity [21] to describe residual stresses caused by growth, thermal expansion, or defects. These will be the focus of our future work.