Geometrically Exact Elastic Shell Crack Propagation Dynamics Based on Fiber Bundles and Differential Forms

1. Introduction

1.1. Geometric Complexity of Shell Fracture and the Need for Dimensional Reduction

Shell structures play a vital role in modern engineering, ranging from aerospace thin-walled structures and pressure vessels to emerging flexible electronic devices. However, the fracture and crack propagation of flexible shells under extreme loading remain a fundamental challenge in continuum mechanics. The fracture mechanics of shell structures differs fundamentally from that of traditional flat plates or three-dimensional blocks; its core difficulty lies in the fact that the stress field at the crack tip presents highly complex mixed-mode characteristics due to the presence of initial curvature.

In flat plate fracture, the stress field at the crack tip mainly includes tension and in-plane shear. However, in shells, due to geometric constraints imposed by curvature, crack propagation inevitably couples out-of-plane bending and torsional loads. This complex mixed-mode behavior makes the prediction of crack propagation direction extremely difficult, and the traditional superposition principle of stress intensity factors often fails. Furthermore, modern engineering structures often experience large displacements and rigid body rotations during service. Under such conditions, traditional fracture mechanics based on linear kinematics and small deformation assumptions leads to severe coordinate distortion.

For a long time, the development of shell fracture mechanics has faced a theoretical debate between "three-dimensional solid reduction" and "direct two-dimensional surface modeling." Many researchers tend to use three-dimensional continuum fracture mechanics to handle shell cracks. However, this approach is geometrically unnatural: the essence of a shell is a two-dimensional manifold with a small thickness. Therefore, there is an urgent need for an intrinsic geometric and coordinate-independent fracture mechanics framework that must be established directly and strictly on the two-dimensional surface representing the shell.

1.2. Integration of Configuration Mechanics and Differential Forms

The most natural theoretical tool for dealing with crack propagation is configuration mechanics. The energy-momentum tensor (the Eshelby tensor $\mathbb{E}$) proposed by Eshelby in the 1950s [12,13] can accurately describe the generalized forces acting on internal material defects. In the development of modern geometric mechanics, the Eshelby tensor is no longer a mere algebraic expression but a geometric invariant that can be naturally derived through the pullback metric of a manifold and variational principles [18,19,22,23].

Significant contributions to configurational mechanics in the context of fracture have been made by Steinmann and colleagues, who developed computational frameworks for material forces and their application to fracture mechanics [35,36,37]. The variational approach to fracture, pioneered by Francfort and Marigo [14], reformulated Griffith’s criterion as a variational problem, providing a rigorous mathematical foundation for crack propagation. Bourdin, Francfort, and Marigo [4,5] further developed numerical implementations of this variational approach using regularization techniques.

However, applying configuration mechanics to curved shell fracture still faces the algebraic complexity of local coordinate transformation and curvature coupling. Traditional shell theories often involve complex differential geometry concepts. To simplify calculations, Sun et al. [38], based on the work of Simo and Fox [32,33,34], reviewed and developed a single-director shell finite deformation model without complex geometric concepts. The motivation of this paper is to thoroughly internalize the concepts of configuration mechanics into the modern geometric framework of fiber bundles and exterior differential forms, while strictly limiting all derivations to the two-dimensional shell mid-surface manifold $\mathcal{M}$. The crack front is defined as a one-dimensional submanifold on $\mathcal{M}$, and its normal propagation direction and velocity are determined entirely by the configuration forces acting on it.

Differential forms and exterior calculus have brought a revolutionary way of expression to continuum mechanics [15]. By converting the two-dimensional stress tensor into a stress 1-form, the mechanical equilibrium equations can be expressed concisely as closed form conditions, and boundary integrals naturally satisfy Stokes’ theorem. For fracture mechanics, this means that the J-integral can be expressed as the flux integral of the configuration stress 1-form along any loop on the mid-surface manifold.

The Discrete Exterior Calculus (DEC) framework developed by Desbrun et al. [8,21] and the Finite Element Exterior Calculus (FEEC) of Arnold, Falk, and Winther [1,2] provide natural discretization schemes for structure-preserving numerical methods. These frameworks have been applied to various problems in continuum mechanics, including elasticity [39] and fluid dynamics [24].

This paper aims to extend the previously established geometrically exact Kirchhoff-Love shell theory based on fiber bundles and differential forms to a two-dimensional manifold topology containing evolving cracks. By introducing the evolving mid-surface, moving frames, and configuration variations, we construct a geometrically exact Eshelby stress 1-form and derive the emergence equations for crack propagation direction and velocity on the tangent space of the manifold.

Before proceeding to the detailed geometric formulation, it is worth clarifying a fundamental point regarding the dimensional reduction strategy adopted herein. Our theory is established intrinsically on the two-dimensional midsurface manifold, with all variational and differential-geometric operations performed on the tangent bundle of the midsurface. Through the Kirchhoff–Love kinematic hypothesis, the through-thickness effects are condensed into the membrane and bending rigidities.

2. Manifold Topology and Geometry of Cracked Shells

2.1. Manifold Structure of Evolving Mid-Surface and Crack Front

In the modern geometric formulation of continuum mechanics, the theory of elastic shells is strictly established on a two-dimensional mid-surface manifold $\mathcal{M}$. In the initial defect-free state, the manifold $\mathcal{M}$ is usually assumed to be a compact and smooth two-dimensional Riemannian manifold. However, when crack propagation occurs inside the shell, the topological structure of the manifold undergoes fundamental changes.

To describe this evolutionary process rigorously, we define the reference mid-surface manifold as a topological space that evolves with time. Let ${\mathcal{M}}_0$ be the initial, uncracked reference midsurface. The cracked reference manifold at time t is defined as:

$${\mathcal{M}}_t={\mathcal{M}}_0\setminus {\mathcal{C}}_t,$$

(1)

where ${\mathcal{C}}_t\subset {\mathcal{M}}_0$ is the open set representing the crack surface at time t. The crack surface is assumed to grow monotonically: ${\mathcal{C}}_s\subset {\mathcal{C}}_t$ for $s<t$, reflecting the irreversibility of fracture. The crack front (crack tip) is defined as the one-dimensional submanifold:

$$\Lambda \left(t\right)=\partial {\mathcal{C}}_t\subset {\mathcal{M}}_t.$$

(2)

The boundary $\partial {\mathcal{M}}_t$ of the cracked shell is divided into two disjoint parts: the external boundary $\partial {\mathcal{M}}_{\mathrm{ext}}$ and the internal crack boundary (i.e., crack surface) $\partial {\mathcal{M}}_{\mathrm{crack}}\left(t\right)$:

$$\partial {\mathcal{M}}_t=\partial {\mathcal{M}}_{\mathrm{ext}}\cup \partial {\mathcal{M}}_{\mathrm{crack}}\left(t\right),\partial {\mathcal{M}}_{\mathrm{crack}}\left(t\right)=\partial {\mathcal{C}}_t=\Lambda \left(t\right).$$

(3)

The external boundary is usually fixed or subjected to known loads, while the internal crack boundary evolves freely over time. From the perspective of differential topology, the existence of $\Lambda \left(t\right)$ divides the locally smooth two-dimensional manifold into cracked and uncracked regions.

2.2. Crack Surface Boundary Conditions and Moving Frames

On the reference configuration, the crack surface $\partial {\mathcal{M}}_{\mathrm{crack}}\left(t\right)$ is treated as a completely free surface that does not bear any external surface tractions. This free boundary condition corresponds to the homogeneous Neumann boundary condition in mechanics. In the context of geometrically exact shell theory, this means that the membrane force tensor $\mathit{N}$ and the bending moment tensor $\mathit{M}$ at the crack edge must simultaneously satisfy zero projection conditions:

$${\mathit{N}}^{ab}{\nu }_a=0,{\mathit{M}}^{ab}{\nu }_a=0\mathrm{on}\partial {\mathcal{M}}_{\mathrm{crack}}\left(t\right),$$

(4)

where $\mathit{\nu }$ is the unit outward normal vector of the crack surface on ${\mathcal{M}}_t$, and $a,b$ are the mid-surface coordinate indices; ${\mathit{N}}^{ab}$ and ${\mathit{M}}^{ab}$ are the components of the membrane stress tensor and bending moment tensor defined on the two-dimensional manifold [32,38].

To maintain the coordinate-independence of the theory, we employ the orthonormal moving frame $\left\{{\Theta }^I\right\}$ and its corresponding connection 1-forms ${\omega }_J^I$ on ${\mathcal{M}}_t\setminus \Lambda \left(t\right)$ [15]. The moving frame method attaches a set of orthonormal bases $\left\{{\mathit{E}}_I\right\}$ to the tangent space at each point of the manifold, expressing tensor components as invariants under local orthonormal bases.

2.3. Fiber Bundle Characterization of Shell Geometry and Pullback Metrics

To reveal more deeply the coupling between large shell deformations and crack geometry, we establish the kinematics of the shell on the structure of a fiber bundle. Let the two-dimensional reference mid-surface manifold be $\mathcal{M}$, endowed with a reference metric tensor $\mathit{G}$. The current configuration is embedded in Euclidean space ${\mathbb{R}}^3$, with the metric tensor denoted as $\mathit{g}$. The deformation mapping $\mathsf{\Phi }:\mathcal{M}\to {\mathbb{R}}^3$ induces a pullback of the metric:

$$\tilde{\mathit{g}}={\mathsf{\Phi }}^{*}\mathit{g}.$$

(5)

The two-dimensional Green-Lagrange strain tensor $\mathit{E}$ is essentially half the difference between the reference metric and the pullback metric:

$$\mathit{E}={\displaystyle \frac{1}{2}}({\mathsf{\Phi }}^{*}\mathit{g}-\mathit{G}).$$

(6)

In geometrically exact shell theory, the deformation mapping $\mathsf{\Phi }$ is decomposed into a combination of mid-surface displacement $\mathit{\phi }$ and normal rotation $\mathit{n}$. Based on the Kirchhoff-Love hypothesis, the normal rotation $\mathit{n}$ is automatically determined by the differentiation of the mid-surface displacement $\mathit{\phi }$, ensuring that the shell does not undergo transverse shear during deformation.

The presence of cracks directly affects the domain and regularity of the mapping $\mathsf{\Phi }$. In the region ${\mathcal{M}}_t\setminus \partial {\mathcal{M}}_{\mathrm{crack}}\left(t\right)$, $\mathsf{\Phi }$ is smooth; across the crack surface, $\mathsf{\Phi }$ experiences a jump or its gradient exhibits singularity.

2.4. Curvature Coupling Effects and Intrinsic Representation of Mixed Modes

The complexity of shell fracture stems largely from the strong coupling between the local geometric characteristics of the crack front and the initial curvature of the shell. In shells, the non-zero Gaussian curvature and mean curvature of the mid-surface cause the front normal $\mathit{\nu }$ and tangent $\mathsf{\tau }$ to rotate not only within the local tangent plane but also to be constrained by the bending of the surface.

In the moving frame, curvature coupling effects are implicitly contained within the connection 1-forms ${\omega }_J^I$. When the exterior differential operator acts on the moving frame, it automatically couples curvature information into the equilibrium equations via the Cartan structure equations:

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

(7)

Furthermore, the existence of curvature makes the Mode I, II, and III fracture modes at the crack tip no longer orthogonal decompositions in geometric terms. In the geometric framework of this paper, mixed modes manifest as the flux components of the configuration stress 1-form $\mathcal{C}$ in various basis vector directions of the tangent space $T\mathcal{M}$. The classical works of Sih, Paris, and Erdogan [10,28,29] on stress intensity factors for cracked shells provide the foundation for understanding these mixed-mode effects, while our geometric approach provides a coordinate-free reformulation.

3. Kinematics of Crack Propagation

Crack propagation is essentially the dynamic alteration of the topological structure of the reference manifold ${\mathcal{M}}_t$, not merely the displacement of material points in space. To describe this topological evolution, we need to introduce a rigorous kinematic mapping for crack propagation on the tangent bundle of the two-dimensional manifold.

3.1. Front Extension Vector Field and Geometric Flow

Let $\Lambda \left(t\right)$ be the crack front curve at the current time t; it is a one-dimensional submanifold embedded within the mid-surface manifold ${\mathcal{M}}_t$. In the local neighborhood of $\Lambda \left(t\right)$, we can define two orthogonal unit vectors using the mid-surface metric $\mathit{G}$: the tangential unit vector $\mathsf{\tau }$ and the normal unit vector $\mathit{\nu }$, which satisfy the orthonormal relation $\mathit{G}(\mathsf{\tau },\mathit{\nu })=0$.

The evolution of the crack front is driven by its normal velocity field $\mathit{V}\in T\mathcal{M}$, which can be formulated geometrically as:

$${\displaystyle \frac{\mathrm{d}\Lambda }{\mathrm{d}t}}=\mathit{V}=V\mathit{\nu },V>0,$$

(8)

where V is a scalar field representing the magnitude of the local crack propagation velocity, and $\mathit{\nu }$ is a vector field indicating the direction of crack propagation [23]. The velocity field V here is not determined directly by the geometric curvature of the front itself but is controlled by the configuration mechanics equations to be derived later.

3.2. Configuration Variation and Shape Derivative

To obtain the generalized forces driving the evolution of the crack front, we must vary the two-dimensional reference configuration itself. This concept is at the core of configuration mechanics [12,18,23].

Let $\mathit{X}\in {\mathcal{M}}_t$ denote a material point on the reference manifold. Consider a one-parameter family of domain perturbations generated by a smooth vector field $\mathit{V}:{\mathcal{M}}_t\to {\mathbb{R}}^3$:

$${\mathcal{M}}_{t,\epsilon }=\{\mathit{X}+\epsilon \mathit{V}\left(\mathit{X}\right):\mathit{X}\in {\mathcal{M}}_t\}.$$

(9)

The shape derivative of a domain functional $J\left({\mathcal{M}}_t\right)={\int }_{{\mathcal{M}}_t}j\left(\mathit{X}\right)dA$ is given by the Hadamard formula [7,20]:

$${\left.{\displaystyle \frac{d}{d\epsilon }}J\left({\mathcal{M}}_{t,\epsilon }\right)\right|}_{\epsilon =0}={\int }_{{\mathcal{M}}_t}{\displaystyle \frac{\partial j}{\partial \mathit{X}}}·\mathit{V}dA+{\int }_{\partial {\mathcal{M}}_t}j\mathit{V}·\mathit{n}ds,$$

(10)

where $\mathit{n}$ is the outward unit normal to $\partial {\mathcal{M}}_t$.

Unlike the deformation variation (which changes the mapping $\mathsf{\Phi }$), the configuration variation changes the two-dimensional reference metric $\mathit{G}$ and the boundary topology of the manifold. The strain variation induced by the configuration variation can be given precisely by the Lie derivative:

$${\delta }_{\mathit{V}}{\mathit{E}}^{\left(0\right)}=-{\displaystyle \frac{1}{2}}{\mathcal{L}}_{\mathit{V}}\mathit{G}=-{\displaystyle \frac{1}{2}}({\nabla }_a{V}_b+{\nabla }_b{V}_a){\mathit{G}}^a\otimes {\mathit{G}}^b,$$

(11)

where ${\mathcal{L}}_{\mathit{V}}$ is the Lie derivative along the vector field $\mathit{V}$ [15]. The variation of the bending strain tensor is similarly determined by the Lie derivative of the normal vector field.

3.3. Decoupling of Physical Variation and Configuration Variation

In traditional nonlinear elastic variational principles, we usually perform a physical variation $\delta \mathsf{\Phi }$ on the deformation mapping $\mathsf{\Phi }$ to obtain the stress equilibrium equations. However, in fracture mechanics, we need to consider simultaneously the deformation in physical space and the evolution in material space.

The physical variation $\delta \mathsf{\Phi }$ keeps the reference manifold $\mathcal{M}$ and the metric $\mathit{G}$ constant, changing only the spatial mapping; its result is the classical deformation mechanics equations:

$${\int }_{{\mathcal{M}}_t}\left({\mathit{N}}^{ab}\delta {\mathit{E}}_{ab}^{\left(0\right)}+{\mathit{M}}^{ab}\delta {\mathit{K}}_{ab}^{\left(1\right)}\right)dA={\int }_{{\mathcal{M}}_t}{\rho }_{0}h\mathit{b}·\delta \mathit{\phi }dA-{\int }_{{\mathcal{M}}_t}{\rho }_{0}h\ddot{\mathit{\phi }}·\delta \mathit{\phi }dA.$$

(12)

The configuration variation ${\delta }_{\mathit{V}}$, however, keeps the spatial mapping $\mathsf{\Phi }$ constant but changes the reference metric $\mathit{G}$ and the manifold boundary. At the crack front, the tangential component of the configuration variation (along the $\mathsf{\tau }$ direction) represents only a reparameterization of the crack front and produces no actual topological change. In contrast, its normal component (along the $\mathit{\nu }$ direction) represents the real virtual propagation of the crack.

4. Configuration Mechanics and the 2D Eshelby Stress 1-Form

4.1. Variational Derivation of the Geometrically Exact 2D Eshelby Tensor

Configuration mechanics is the natural framework for describing the driving forces of internal material defect evolution [12,13]. To rigorously derive the Eshelby tensor in shell theory, we start directly from the total potential energy functional on the two-dimensional mid-surface manifold $\mathcal{M}$:

$$\Pi (\mathsf{\Phi },{\mathcal{M}}_t)={\int }_{{\mathcal{M}}_t}{\mathsf{\Psi }}_{\mathrm{shell}}({\mathit{E}}^{\left(0\right)},{\mathit{K}}^{\left(1\right)})d{A}_0-W_{\mathrm{ext}},$$

(13)

where ${\mathsf{\Psi }}_{\mathrm{shell}}$ is the surface strain energy density of the shell, and $d{A}_0$ is the area form of the two-dimensional manifold. The variational approach to fracture, as formulated by Francfort and Marigo [14], provides the mathematical foundation for minimizing this energy functional with respect to both the deformation and the crack set.

Let $\mathit{V}$ be a smooth vector field representing a perturbation of the reference configuration. The first variation of the total potential energy with respect to the configuration (domain) is:

$${\delta }_{\mathit{V}}\Pi ={\delta }_{\mathit{V}}{\int }_{{\mathcal{M}}_t}{\mathsf{\Psi }}_{\mathrm{shell}}d{A}_0-{\delta }_{\mathit{V}}{W}_{\mathrm{ext}}.$$

(14)

The variation of the strain energy term consists of three parts: variation of the strain energy density due to metric changes, variation due to boundary movement, and variation of the area element. Using the Hadamard formula (10) and the Lie derivative expression (11), we obtain:

$$\begin{array}{cc}\hfill {\delta }_{\mathit{V}}{\int }_{{\mathcal{M}}_t}{\mathsf{\Psi }}_{\mathrm{shell}}d{A}_0& ={\int }_{{\mathcal{M}}_t}\left({\displaystyle \frac{\partial {\mathsf{\Psi }}_{\mathrm{shell}}}{\partial {\mathit{E}}^{\left(0\right)}}}:{\delta }_{\mathit{V}}{\mathit{E}}^{\left(0\right)}+{\displaystyle \frac{\partial {\mathsf{\Psi }}_{\mathrm{shell}}}{\partial {\mathit{K}}^{\left(1\right)}}}:{\delta }_{\mathit{V}}{\mathit{K}}^{\left(1\right)}\right)d{A}_0\hfill \\ \hfill & +{\int }_{{\mathcal{M}}_t}{\mathsf{\Psi }}_{\mathrm{shell}}\mathrm{div}\left(\mathit{V}\right)d{A}_0+{\int }_{\Lambda \left(t\right)}{\mathsf{\Psi }}_{\mathrm{shell}}\mathit{V}·\mathit{\nu }ds,\hfill \end{array}$$

(15)

where the last term arises from the moving boundary $\Lambda \left(t\right)$ (the crack front). Note that the crack front is a one-dimensional boundary of the domain ${\mathcal{M}}_t$, so the boundary integral is a line integral along $\Lambda \left(t\right)$.

Noting that ${\displaystyle \frac{\partial {\mathsf{\Psi }}_{\mathrm{shell}}}{\partial {\mathit{E}}^{\left(0\right)}}}=\mathit{N}$ (membrane stress) and ${\displaystyle \frac{\partial {\mathsf{\Psi }}_{\mathrm{shell}}}{\partial {\mathit{K}}^{\left(1\right)}}}=\mathit{M}$ (bending moment), and using integration by parts on the interior terms, we obtain:

$${\delta }_{\mathit{V}}\Pi =-{\int }_{{\mathcal{M}}_t}\mathrm{Div}\left(\mathbb{E}\right)·\mathit{V}d{A}_0+{\int }_{\Lambda \left(t\right)}\left[\mathbb{E}\mathit{\nu }\right]·\mathit{V}ds,$$

(16)

where the two-dimensional geometrically exact Eshelby configuration stress tensor is defined as:

$$\mathbb{E}={\mathsf{\Psi }}_{\mathrm{shell}}{\mathit{G}}^{♯}-\mathit{N}\otimes {\mathit{E}}^{\left(0\right)}-\mathit{M}\otimes {\mathit{K}}^{\left(1\right)}.$$

(17)

In component form, this is a mixed tensor:

$${\mathbb{E}}_b^a={\mathsf{\Psi }}_{\mathrm{shell}}{\delta }_b^a-{\mathit{N}}^{ac}{\mathit{E}}_{cb}^{\left(0\right)}-{\mathit{M}}^{ac}{\mathit{K}}_{cb}^{\left(1\right)},$$

(18)

where ${\delta }_b^a$ is the Kronecker delta. This expression accurately reflects the joint contribution of mid-surface stretching and bending deformation mechanisms in the shell to the crack driving force. Since $\mathit{N},\mathit{M},{\mathit{E}}^{\left(0\right)},{\mathit{K}}^{\left(1\right)}$ are all invariant under rigid rotations, $\mathbb{E}$ is an objective two-dimensional geometric invariant.

The configurational force formulation developed by Steinmann [35,36] provides the computational framework for evaluating such Eshelby tensors in finite element settings. Our geometric approach extends this framework to the shell manifold setting.

4.2. Configuration Stress 1-Form and the J-Integral

To perform integral calculations on the two-dimensional manifold, we convert $\mathbb{E}$ into a vector-valued configuration stress 1-form $\mathcal{C}\in {\Omega }^1(\mathcal{M},T^{*}\mathcal{M})$ on the reference manifold. This is achieved via the Hodge star operator ${★}_0$:

$${\mathcal{C}}^a={\mathbb{E}}_b^a{★}_0{\Theta }^b={\mathbb{E}}_b^a{ϵ}_c^b\mathrm{d}{\xi }^c,$$

(19)

where ${ϵ}_c^b$ is the two-dimensional antisymmetric tensor. According to Noether’s theorem, in defect-free regions $\mathrm{d}\mathcal{C}=0$.

We take a closed loop $\Gamma \subset \mathcal{M}$ around the crack front and integrate the configuration stress form to obtain the energy release rate for crack propagation. From the variational result (16), the configurational force on the crack front is:

$${\mathit{f}}_{\mathrm{config}}={\int }_{\Lambda \left(t\right)}\left[\mathbb{E}\mathit{\nu }\right]ds.$$

(20)

The energy release rate $\mathcal{G}$ (the J-integral) is the projection of this configurational force onto the crack propagation direction $\mathit{\nu }$:

$$\mathcal{G}={\mathit{f}}_{\mathrm{config}}·\mathit{\nu }={\int }_{\Lambda \left(t\right)}\mathit{\nu }·\left[\mathbb{E}\mathit{\nu }\right]ds.$$

(21)

Equivalently, by Stokes’ theorem, this can be expressed as a contour integral along any loop $\Gamma$ surrounding the crack front:

$$\mathcal{G}={\int }_{\Gamma }{\mathcal{C}}^a{\nu }_a={\int }_{\Gamma }\left({\mathsf{\Psi }}_{\mathrm{shell}}{\nu }_a-{\mathit{N}}^{cb}{\mathit{E}}_{cb}^{\left(0\right)}{\nu }_a-{\mathit{M}}^{cb}{\mathit{K}}_{cb}^{\left(1\right)}{\nu }_a\right)ds,$$

(22)

where ${\nu }_a$ is the component of the unit normal to the contour $\Gamma$ on the manifold, and $ds$ is the loop arc length element. Since $\mathrm{d}\mathcal{C}=0$, the integral $\mathcal{G}$ is independent of the path $\Gamma$ by Stokes’ theorem.

This geometric formulation of the J-integral generalizes the classical path-independent integral of Rice [27] and Cherepanov [6] to curved shell manifolds. The shell J-integral has been extensively studied in the shell fracture literature [26,30], and our geometric formulation provides a coordinate-invariant generalization.

4.3. Physical Decomposition and Objectivity of Shell Configuration Forces

Shell fracture is often accompanied by complex mixed modes. In the geometrically exact framework, this physical decomposition can be naturally realized through the components of the Eshelby tensor. Since $\mathbb{E}$ is a mixed tensor, its symmetric and antisymmetric components correspond respectively to the driving forces of different fracture modes. The second Piola-Kirchhoff stress $\mathit{S}$ and Green-Lagrange strain $\mathit{E}$ are both invariant under rigid rotations, and the reference metric $\mathit{G}$ is also objective; therefore $\mathbb{E}$ is likewise an objective geometric invariant [32,38].

The work of Sih and Paris [29,30] on the fracture mechanics of shells established the importance of correctly accounting for mixed-mode effects in curved structures. Our geometric approach naturally captures these effects through the structure of the Eshelby tensor on the curved manifold.

4.4. Topological Singularities in Cracked Manifolds and Configuration Force Divergence

In a continuum, the divergence $\mathrm{d}\mathcal{C}$ of the configuration stress 1-form can be understood as a "source of configuration forces" on the manifold. In homogeneous, defect-free regions, Noether’s theorem guarantees $\mathrm{d}\mathcal{C}=0$ due to the translational symmetry of the reference configuration. However, the presence of a crack destroys this symmetry.

At the crack front $\Lambda \left(t\right)$, a topological mutation occurs in the manifold. This topological mutation mathematically leads to discontinuities or divergences in the configuration stress form $\mathcal{C}$ as it crosses $\Lambda \left(t\right)$. The jump $\left[\mathbb{E}\mathit{\nu }\right]$ across the crack front gives the configurational force density per unit length along the front:

$${\mathit{f}}_{\mathrm{config}}=\left[\mathbb{E}\mathit{\nu }\right]\mathrm{on}\Lambda \left(t\right).$$

(23)

The crack front $\Lambda \left(t\right)$ thus serves as a "line source" of configuration forces. The J-integral $\mathcal{G}$ extracts the normal component of this line source.

5. Crack Propagation Laws: Emergence of Direction and Velocity

5.1. Griffith Criterion and Energy Balance

The cornerstone of continuum fracture mechanics is the energy balance theory proposed by Griffith [17]. Let $G_{\mathrm{c}}$ be the critical energy release rate (fracture toughness) of the shell material. The Griffith criterion states:

$$\mathcal{G}\le G_{\mathrm{c}},\mathrm{and}\mathcal{G}=G_{\mathrm{c}}\mathrm{when}\mathrm{the}\mathrm{crack}\mathrm{propagates}.$$

(24)

From the perspective of irreversible process thermodynamics, crack propagation is a dissipative process. The dissipation power density can be expressed as:

$$\mathcal{D}=(\mathcal{G}-G_{\mathrm{c}})V\ge 0.$$

(25)

The variational formulation of Griffith’s criterion, as developed by Francfort and Marigo [14], provides a rigorous mathematical framework for understanding crack propagation as an energy minimization problem. Bourdin et al. [4,5] demonstrated the numerical implementation of this variational approach using phase-field regularization, which has become a standard method for computational fracture mechanics.

5.2. Determination of Propagation Direction: Maximum Energy Release Rate Principle

At any point on the crack front $\Lambda \left(t\right)$, the propagation direction vector $\mathit{\nu }$ is determined by the principle of maximization of the local energy release rate. The energy release rate as a function of the crack extension direction $\theta$ (measured from the tangent direction $\mathsf{\tau }$) is given by the projection of the configurational force:

$$\mathcal{G}\left(\theta \right)={\mathit{f}}_{\mathrm{config}}·(cos\theta \mathsf{\tau }+sin\theta \mathit{\nu }).$$

(26)

The maximum energy release rate principle requires:

$${\displaystyle \frac{\partial \mathcal{G}\left(\theta \right)}{\partial \theta }}=0,{\displaystyle \frac{{\partial }^2\mathcal{G}\left(\theta \right)}{\partial {\theta }^2}}<0.$$

(27)

Solving these conditions yields:

$$tan{\theta }^{*}={\displaystyle \frac{{\mathit{f}}_{\mathrm{config}}·\mathit{\nu }}{{\mathit{f}}_{\mathrm{config}}·\mathsf{\tau }}},$$

(28)

or equivalently, the optimal propagation direction is:

$${\mathit{\nu }}^{*}={\displaystyle \frac{{\mathit{f}}_{\mathrm{config}}-({\mathit{f}}_{\mathrm{config}}·\mathsf{\tau })\mathsf{\tau }}{\parallel {\mathit{f}}_{\mathrm{config}}-({\mathit{f}}_{\mathrm{config}}·\mathsf{\tau })\mathsf{\tau }\parallel }}.$$

(29)

This derivation is performed entirely within the intrinsic space of the two-dimensional manifold, automatically handling the coupling problem of mixed-mode cracks through the projection of differential forms.

The maximum energy release rate criterion for mixed-mode fracture has been extensively studied in the literature [10,25,31]. Our geometric formulation provides a coordinate-invariant implementation of this criterion on curved shell manifolds.

5.3. Evolution Equations for Propagation Velocity

When $\mathcal{G}\left({\mathit{\nu }}^{*}\right)>G_{\mathrm{c}}$, the crack propagates dynamically. During dynamic propagation, the motion of the crack front excites inertial effects, and kinetic energy terms consume part of the driving energy. Following Freund [16], the dynamic energy release rate is:

$${\mathcal{G}}_{\mathrm{dyn}}=\mathcal{G}·f\left(V\right),$$

(30)

where $f\left(V\right)$ is a universal function of crack speed (typically $f\left(V\right)=1-V/c_R$ for mode I, with $c_R$ the Rayleigh wave speed).

To establish a closed dynamic equation, we adopt a power-law model:

$$V=C_m{\left({\displaystyle \frac{{\mathcal{G}}_{\mathrm{dyn}}\left({\mathit{\nu }}^{*}\right)-G_{\mathrm{c}}}{G_{\mathrm{c}}}}\right)}^m,$$

(31)

where $C_m$ is a material-dependent parameter for maximum crack propagation velocity, and m is the dynamic exponent. This evolution equation, combined with the preceding geometric kinematic equations, constitutes a closed and objective dynamic system for crack propagation.

The dynamic fracture mechanics framework developed by Freund [16] and the cohesive zone models of Barenblatt [3] and Dugdale [9] provide complementary approaches to modeling crack propagation velocities. Our geometric framework can accommodate these various constitutive laws for crack velocity.

5.4. Intrinsicality of the Geometric Fracture Criterion and Numerical Advantages

The crack propagation law proposed in this paper is entirely built on the intrinsic geometric structure of the manifold. In the geometrically exact framework, the configurational force vector ${\mathit{f}}_{\mathrm{config}}$ is obtained directly by integrating the Eshelby stress 1-form $\mathcal{C}$; its direction directly indicates the manifold path of fastest energy release. The propagation direction ${\mathit{\nu }}^{*}$ is obtained through simple tangential projection without complex eigenvalue analysis. This intrinsicality makes this fracture criterion particularly suitable for implementation in the Discrete Exterior Calculus (DEC) framework.

6. Coupled Dynamic Equations and Topology Updates

Synthesizing the preceding geometric and mechanical derivations, the complete geometrically exact shell crack propagation dynamic system consists of three mutually coupled subsystems: deformation dynamics, configuration calculation, and propagation evolution.

6.1. Geometric Discretization of the Deformation Dynamics Subsystem

On the fixed current topology ${\mathcal{M}}_t$, the deformation dynamics subsystem must first be solved, i.e., finding the mid-surface mapping $(\mathit{\phi },\mathit{n})$ that satisfies the principle of virtual work with cracked free boundary conditions. In geometrically exact shell theory, this weak form is:

$${\int }_{{\mathcal{M}}_t}\left({\mathit{N}}^{ab}\delta {\mathit{E}}_{ab}^{\left(0\right)}+{\mathit{M}}^{ab}\delta {\mathit{K}}_{ab}^{\left(1\right)}\right)d{A}_0={\int }_{{\mathcal{M}}_t}{\rho }_{0}h\mathit{b}·\delta \mathit{\phi }d{A}_0-{\int }_{{\mathcal{M}}_t}{\rho }_{0}h\ddot{\mathit{\phi }}·\delta \mathit{\phi }d{A}_0.$$

(32)

The crack surface boundary conditions ${\mathit{N}}^{ab}{\nu }_a=0$ and ${\mathit{M}}^{ab}{\nu }_a=0$ are homogeneous Neumann conditions, automatically included as natural boundary conditions.

6.2. Local Extraction of the Configuration Calculation Subsystem

Based on the deformation field solved at the current moment, we calculate the Eshelby stress 1-form $\mathcal{C}$ within a local tubular neighborhood of the crack front $\Lambda \left(t\right)$. By integrating along the closed loop $\Gamma$ surrounding the front, we obtain the direction-dependent generalized J-integral $\mathcal{G}\left(\mathit{\nu }\right)$:

$$\mathcal{G}\left(\mathit{\nu }\right)=\sum _{e\in \Gamma }{\mathcal{C}}_e^a{\nu }_a.$$

(33)

The configurational force principal vector is reconstructed as:

$${\mathit{f}}_{\mathrm{config}}=\sum _{e\in \Gamma }{\mathcal{C}}_e^a{\mathit{E}}_a,$$

(34)

where $\left\{{\mathit{E}}_a\right\}$ is the discrete orthonormal frame on the edge.

6.3. Topology Updates of the Propagation Evolution Subsystem

Using the maximum energy release rate principle to solve for the variational extremum, we first obtain the optimal propagation direction ${\mathit{\nu }}^{*}$; subsequently, we calculate the scalar velocity V according to the dynamic propagation law (31). The topological update of the reference manifold is governed by the front evolution equation:

$${\displaystyle \frac{\mathrm{d}\Lambda \left(t\right)}{\mathrm{d}t}}=V{\mathit{\nu }}^{*}.$$

(35)

In numerical implementation, this corresponds to adaptive reconstruction of the discrete simplicial complex.

6.4. Numerical Advantages of Structure-Preserving Algorithms

This decoupling-coupling process of topology-deformation has natural advantages in the Discrete Exterior Calculus (DEC) or Finite Element Exterior Calculus (FEEC) frameworks [1,8,21]. Since discrete differential operators maintain strict algebraic structures when the mesh topology changes, solving mechanical equilibrium does not require modifying the core matrices of the discrete format. This structure-preserving algorithm guarantees conservation of energy and momentum and improves robustness in handling large-scale dynamic crack propagation.

7. Discrete Exterior Calculus Implementation Framework

This section presents a detailed discretization scheme based on Discrete Exterior Calculus (DEC) [8,21,24]. Unlike traditional finite element methods, DEC preserves the fundamental algebraic structures of exterior calculus—the generalized Stokes theorem ${\int }_{\Omega }\mathrm{d}\omega ={\int }_{\partial \Omega }\omega$ holds exactly on the discrete level.

7.1. Discrete Manifold Representation

The continuous shell midsurface $\mathcal{M}$ is approximated by a simplicial complex—a triangulated surface ${\mathcal{M}}_h=(V,E,F)$ consisting of vertices V, oriented edges E, and oriented triangular faces F. The crack is modeled by splitting the nodes along the crack path in the reference configuration, creating free internal boundaries $\partial {\mathcal{M}}_{\mathrm{crack}}$. The crack front $\Lambda \left(t\right)$ is represented as a chain of connected edges in the simplicial complex.

7.1.1. Discrete Differential Operators

The fundamental operators of exterior calculus admit exact discrete analogues on ${\mathcal{M}}_h$:

• Discrete exterior derivative ${\mathsf{d}}_k$: For a 0-form f, ${\left({\mathsf{d}}_{0}f\right)}_e=f_j-f_i$ on oriented edge $e=(i,j)$. For a 1-form $\alpha$, ${\left({\mathsf{d}}_1\alpha \right)}_f={\alpha }_{e_1}+{\alpha }_{e_2}+{\alpha }_{e_3}$ on face f. These satisfy the discrete Poincaré lemma: $D_1{D}_0=0$.
• Discrete Hodge star ${★}_k$: For a 0-form, ${★}_0$ is the diagonal matrix of dual vertex areas $V_i^{*}={\displaystyle \frac{1}{3}}{\sum }_{f\ni i}{A}_f$. For 1-forms, ${★}_1$ is diagonal with entries $|e^{*}|/|e|$, where $|e^{*}|={\displaystyle \frac{1}{2}}{\sum }_{f\supset e}{A}_f/h_f$ is the dual edge length and $h_f$ is the altitude to edge e in face f.

7.2. Discrete Kinematics

The deformation mapping $\mathsf{\Phi }:{\mathcal{M}}_h\to {\mathbb{R}}^3$ is discretized by vertex displacements ${\mathit{u}}_i\in {\mathbb{R}}^3$, so that the current position is ${\mathit{x}}_i^{\mathrm{cur}}={\mathit{x}}_i^{\mathrm{ref}}+{\mathit{u}}_i$.

7.2.1. Membrane Strain

On each triangular face $f=(i,j,k)$, for edge $e=(i,j)$, the Green-Lagrange strain along the edge is:

$$E_e^{\left(0\right)}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\ell }_{\mathrm{cur},e}^2-{\ell }_{\mathrm{ref},e}^2}{{\ell }_{\mathrm{ref},e}^2}}\right),$$

(36)

where ${\ell }_{\mathrm{ref},e}=\parallel {\mathit{x}}_j^{\mathrm{ref}}-{\mathit{x}}_i^{\mathrm{ref}}\parallel$ and ${\ell }_{\mathrm{cur},e}=\parallel ({\mathit{x}}_j^{\mathrm{ref}}+{\mathit{u}}_j)-({\mathit{x}}_i^{\mathrm{ref}}+{\mathit{u}}_i)\parallel$. A least-squares reconstruction on each face yields the full $2\times 2$ covariant membrane strain tensor ${\mathit{E}}_f^{\left(0\right)}$.

7.2.2. Bending Strain

The bending strain ${\mathit{K}}_f^{\left(1\right)}$ is computed from the variation of vertex normals. Let ${\mathit{n}}_i$ be the discrete unit normal at vertex i. For edge $e=(i,j)$, the bending strain component is:

$$K_e^{\left(1\right)}={\displaystyle \frac{({\mathit{n}}_j-{\mathit{n}}_i)·{\mathit{t}}_e}{{\ell }_{\mathrm{ref},e}}},$$

(37)

where ${\mathit{t}}_e$ is the unit tangent vector along edge e in the reference configuration.

7.3. Discrete Equilibrium Equations

For a fixed topology, the principle of virtual work is discretized as: find $\mathit{u}\in {\mathbb{R}}^{3N_V}$ such that for all virtual displacements $\delta \mathit{u}$,

$$\sum _{f\in F}{A}_f\left({\mathit{N}}_f:\delta {\mathit{E}}_f^{\left(0\right)}+{\mathit{M}}_f:\delta {\mathit{K}}_f^{\left(1\right)}\right)=\sum _{f\in F}{A}_f{\rho }_{0}h\mathit{b}·\delta {\mathit{u}}_f-\sum _{f\in F}{A}_f{\rho }_{0}h{\ddot{\mathit{u}}}_f·\delta {\mathit{u}}_f,$$

(38)

where ${\mathit{N}}_f=\partial {\mathsf{\Psi }}_{\mathrm{shell}}/\partial {\mathit{E}}_f^{\left(0\right)}$ and ${\mathit{M}}_f=\partial {\mathsf{\Psi }}_{\mathrm{shell}}/\partial {\mathit{K}}_f^{\left(1\right)}$. For linear elasticity:

$${\mathit{N}}_f={\mathbf{D}}_{\mathrm{membrane}}:{\mathit{E}}_f^{\left(0\right)},{\mathit{M}}_f={\mathbf{D}}_{\mathrm{bending}}:{\mathit{K}}_f^{\left(1\right)},$$

(39)

with ${\mathbf{D}}_{\mathrm{membrane}}={\displaystyle \frac{Eh}{1-{\nu }^2}}\mathcal{I}$ and ${\mathbf{D}}_{\mathrm{bending}}={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}\mathcal{I}$.

Linearization yields the discrete Newton system:

$${\mathit{K}}_{\tan \mathrm{gent}}\left({\mathit{u}}^{\left(k\right)}\right)\Delta \mathit{u}={\mathbf{F}}_{\mathrm{ext}}-{\mathbf{F}}_{\mathrm{int}}\left({\mathit{u}}^{\left(k\right)}\right),$$

(40)

where ${\mathit{K}}_{\tan \mathrm{gent}}={\sum }_{f\in F}\left({\mathit{K}}_{\mathrm{membrane}}^{\left(f\right)}+{\mathit{K}}_{\mathrm{bending}}^{\left(f\right)}+{\mathit{K}}_{\mathrm{geo}}^{\left(f\right)}\right)$.

7.4. Discrete Eshelby Stress and the J-Integral

On each triangular face f, the discrete Eshelby tensor is computed as:

$${\mathbb{E}}_{b,f}^a={\mathsf{\Psi }}_f{\delta }_b^a-N_f^{ac}{E}_{cb,f}^{\left(0\right)}-M_f^{ac}{K}_{cb,f}^{\left(1\right)},$$

(41)

where ${\mathsf{\Psi }}_f={\displaystyle \frac{1}{2}}({\mathit{N}}_f:{\mathit{E}}_f^{\left(0\right)}+{\mathit{M}}_f:{\mathit{K}}_f^{\left(1\right)})$.

The configuration stress 1-form $\mathcal{C}$ is discretized on each edge $e\in \Gamma$ as:

$${\mathcal{C}}_e^a=\sum _{f\supset e}{\displaystyle \frac{1}{2}}{\mathbb{E}}_{b,f}^a{\left({★}_0{\Theta }^b\right)}_e.$$

(42)

7.5. Algorithmic Summary

The complete numerical procedure is summarized in Algorithm 1.

Algorithm 1 DEC solution procedure for geometrically exact shell fracture

1: procedurePreprocessing 2:     Generate triangulated shell midsurface ${\mathcal{M}}_h$ with refinement near crack tips 3:     Split nodes along crack paths to enforce free boundary conditions 4:     Construct discrete operators: $D_0$, $D_1$, ${★}_0$, ${★}_1$, ${★}_2$ 5:     Identify crack front chain $\Lambda \left(t\right)$ and integration loop $\Gamma$ 6: end procedure 7: procedureMechanicalEquilibrium(${\mathit{u}}^{\left(0\right)}$) 8:     for $k=0,1,2,\dots$ until convergence do 9:         Compute face strains ${\mathit{E}}_f^{\left(0\right)},{\mathit{K}}_f^{\left(1\right)}$ from ${\mathit{u}}^{\left(k\right)}$ 10:         Compute stresses ${\mathit{N}}_f,{\mathit{M}}_f$ and assemble ${\mathbf{F}}_{\mathrm{int}}$ 11:         Assemble tangent stiffness ${\mathit{K}}_{\tan \mathrm{gent}}$ 12:         Solve ${\mathit{K}}_{\tan \mathrm{gent}}\Delta \mathit{u}={\mathbf{F}}_{\mathrm{ext}}-{\mathbf{F}}_{\mathrm{int}}$ 13:         Update ${\mathit{u}}^{(k+1)}={\mathit{u}}^{\left(k\right)}+\Delta \mathit{u}$ 14:     end for 15:     return $\mathit{u}$ 16: end procedure 17: procedureFracturePostProcessing($\mathit{u}$) 18:     for each face f do 19:         Compute ${\mathbb{E}}_{b,f}^a$ via (41) 20:     end for 21:     for each edge $e\in \Gamma$ do 22:         Compute ${\mathcal{C}}_e^a$ and accumulate $\mathcal{G}\left(\mathit{\nu }\right)$ via (33) 23:     end for 24:     Compute ${\mathit{f}}_{\mathrm{config}}$ and optimal direction ${\mathit{\nu }}^{*}$ via (29) 25: end procedure 26: Preprocessing 27: $\mathit{u}\leftarrow$MechanicalEquilibrium($\mathbf{0}$) 28: FracturePostProcessing($\mathit{u}$)

8. Conclusion

This paper constructs a rigorous, intrinsic, and coordinate-independent geometric mechanics framework for fracture and crack propagation in flexible shells. The main contributions are:

1.

Strict characterization of the manifold topology of cracked shells: The cracked shell is modeled as an evolving manifold ${\mathcal{M}}_t={\mathcal{M}}_0\setminus {\mathcal{C}}_t$ with crack front $\Lambda \left(t\right)=\partial {\mathcal{C}}_t$.

2.

Complete variational derivation of the energy release rate: Using the Hadamard shape derivative formula, the first variation of the total potential energy is derived rigorously, yielding the configurational force ${\mathit{f}}_{\mathrm{config}}=\left[\mathbb{E}\mathit{\nu }\right]$ on the crack front.

3.

Geometrically exact 2D Eshelby stress tensor: The tensor $\mathbb{E}={\mathsf{\Psi }}_{\mathrm{shell}}{\mathit{G}}^{♯}-\mathit{N}\otimes {\mathit{E}}^{\left(0\right)}-\mathit{M}\otimes {\mathit{K}}^{\left(1\right)}$ is derived and expressed as a vector-valued 1-form.

4.

Intrinsic crack propagation laws: The propagation direction is determined by the maximum energy release rate principle, yielding ${\mathit{\nu }}^{*}={\displaystyle \frac{{\mathit{f}}_{\mathrm{config}}-({\mathit{f}}_{\mathrm{config}}·\mathsf{\tau })\mathsf{\tau }}{\parallel \dots \parallel }}$. The velocity follows a power-law model $V=C_m{(({\mathcal{G}}_{\mathrm{dyn}}-G_{\mathrm{c}})/G_{\mathrm{c}})}^m$.

5.

DEC discretization framework: A complete discrete framework preserving the algebraic structure of the continuous theory is developed.

The profound significance of this theoretical framework lies in implicitly containing the complex curvature-fracture coupling within the structures of exterior differentiation and pullback metrics, ensuring that the objectivity axiom of continuum mechanics holds strictly on manifolds with evolving defects.

Looking to the future, the geometric framework proposed in this paper paves the way for structure-preserving adaptive dynamic fracture algorithms. Further research may explore generalization to multiphysics coupling and non-local effects in geometric fracture mechanics.