A discontinuous Galerkin / cohesive zone model approach for the computational modeling of fracture in geometrically exact slender beams

TL;DR

The paper addresses fracture modeling in slender beams undergoing large deformations, where buckling interacts with fracture. It develops a discontinuous Galerkin discretization of a geometrically exact torsion-free Kirchhoff beam, with inter-element jumps in position and tangent DOFs and a cohesive-zone model that activates at fracture onset via a stress-resultant criterion $f_{eq}(\langle\boldsymbol{f}\rangle,\langle\boldsymbol{m_{\perp}}\rangle) \ge f_c$. The DG/CZM framework yields energy-dissipative fracture behavior for tensile and bending modes and is validated against buckling benchmarks and the Audoly–Neukirch spaghetti experiments, while also reproducing spall and transverse-load fracture scenarios. This approach enables accurate, physics-based simulation of buckling–fracture transitions in slender-beam lattices and metamaterials under large deformations.

Abstract

Slender beams are often employed as constituents in engineering materials and structures. Prior experiments on lattices of slender beams have highlighted their complex failure response, where the interplay between buckling and fracture plays a critical role. In this paper, we introduce a novel computational approach for modeling fracture in slender beams subjected to large deformations. We adopt a state-of-the-art geometrically exact Kirchhoff beam formulation to describe the finite deformations of beams in three-dimensions. We develop a discontinuous Galerkin finite element discretization of the beam governing equations, incorporating discontinuities in the position and tangent degrees of freedom at the inter-element boundaries of the finite elements. Before fracture initiation, we enforce compatibility of nodal positions and tangents weakly, via the exchange of variationally-consistent forces and moments at the interfaces between adjacent elements. At the onset of fracture, these forces and moments transition to cohesive laws modeling interface failure. We conduct a series of numerical tests to verify our computational framework against a set of benchmarks and we demonstrate its ability to capture the tensile and bending fracture modes in beams exhibiting large deformations. Finally, we present the validation of our framework against fracture experiments of dry spaghetti rods subjected to sudden relaxation of curvature.

Figures (14)

• Figure 1: Illustration of the kinematics of a geometrically exact Kirchhoff beam in the case of a beam with circular cross-section. The beam configuration is characterized by the line of centroids $\boldsymbol{r}(s)$ and by the orthonormal intrinsic frame $\{\boldsymbol{g_1}(s), \boldsymbol{g_2}(s), \boldsymbol{g_3}(s)\}$. The Kirchhoff constraint enforces that $\boldsymbol{g_1}(s)$ be tangent to the line of centroids, that is $\boldsymbol{g_1}(s) = \frac{\boldsymbol{r}^{\prime}(s)}{||\boldsymbol{r}^{\prime}(s)||}$.
• Figure 2: The fracture modes of beams: (a) tension, (b) shear, (c) torsion and (d) bending.
• Figure 3: The force-separation cohesive law prescribes a linear decay of the scalar effective cohesive force $f_{coh}$ with the scalar effective separation $\Delta$, from a critical value $f_c$ to zero. Irreversibility is modeled by introducing a history variable $\Delta_{max}$ representing the maximum effective separation achieved. The unloading path follows a trajectory back to the origin, while reloading occurs on the same unloading path. The total area under the force-separation curve is equal to the effective fracture energy of the material while the black and grey areas at any point $\Delta = \Delta_{max}$ represent the dissipated and maximum recoverable energies at the cohesive boundary. Complete fracture is achieved when $\Delta \geq \Delta_{c}$.
• Figure 4: The discontinuous Galerkin discretization of the straight undeformed beam. Internal nodes, e.g. nodes $s_1^e$ and $s_0^{e+1}$, are duplicated to allow the embedding of potential discontinuities at the element interfaces $s_n$ for $n = 1, ..., E-1$.
• Figure 5: Convergence analysis of the discontinuous Galerkin (DG) finite element discretization for the bending of a cantilever beam benchmark. The plot shows the relative $L^2$ error in the deformed beam centerline positions as a function of the mesh size for different penalty parameter values ($\beta_{p,n}$ and $\beta_{t,n}$). Reference lines representing the convergence orders two and four are shown with dotted lines. The convergence order of the DG discretization approaches four, which is the expected value for a continuous Galerkin (CG) discretization, for higher values of the penalty parameters. Convergence deteriorates for smaller values of the penalty parameters.