Cohesive phase-field fracture with an explicit strength surface: an eigenstrain-based return-mapping formulation

Technical Summary: Cohesive Phase-Field Fracture with an Explicit Strength Surface

Problem Statement
Standard phase-field fracture methods are fundamentally rooted in Griffith's brittle fracture theory. Consequently, they lack an inherent prescription for material strength, making crack nucleation dependent on ad-hoc criteria or the regularization parameter (phase-field length scale) rather than a physical yield surface. Furthermore, these methods struggle to capture cohesive fracture behavior and ductile responses, as the degradation functions are typically calibrated for uniaxial tension and do not naturally extend to multi-axial loading or pressure-dependent strengthening. While recent eigenstrain-based formulations (e.g., Vicentini et al., Bourdin et al.) successfully decouple material strength from fracture energy by introducing fracture eigenstrains, their implementation relies on global energy-minimization frameworks. This approach treats eigenstrains as additional field variables, introducing extra global degrees of freedom and requiring specialized software, which hinders integration into standard finite-element (FE) codes that evaluate constitutive models pointwise at integration points.

Methodology
This work reformulates the eigenstrain-based cohesive phase-field model as a local constitutive model, analogous to plasticity return-mapping schemes, to enable integration into standard FE codes without additional global degrees of freedom.

• Governing Equations: The framework utilizes a total energy functional comprising elastic energy (based on displacement and fracture eigenstrains), fracture energy (distributed via a phase-field parameter $\phi$), and a strength potential. The elastic strain is defined as $\varepsilon^{el} = \varepsilon - \eta$, where $\eta$ represents the fracture eigenstrain.
• Local Constitutive Formulation: A key innovation is the resolution of eigenstrain evolution entirely at the integration-point level. The eigenstrains are defined as $\eta = \lambda G(\varepsilon)$, where $\lambda$ is a multiplier and $G$ is a direction function dependent on the strain state. This allows the constitutive update to be solved locally using a Newton-Raphson scheme, similar to plasticity, rather than globally.
• Strength Criteria: Two distinct strength potentials are implemented:
  1. Non-smooth $r_1$ Criterion: Features independent tensile ($f_t$) and shear ($f_s$) strengths with a high compressive strength to enforce non-penetration.
  2. Smooth Drucker-Prager-like Criterion: Captures pressure-dependent strengthening under compression. It uses a reference strain parameter ($\varepsilon_{ref}$) to dictate the transition between tensile and compressive regimes, ensuring a smooth failure surface.
• Algorithmic Implementation: The authors derive consistent tangent operators for both smooth and non-smooth criteria to ensure robust convergence of the global Newton-Raphson solver. A staggered solution scheme is employed, solving for displacements and eigenstrains first, followed by the phase-field variable, iterating until convergence. The implementation uses quadratic Lagrange elements and is coded in FEniCSx.

Key Contributions

1. Local Return-Mapping Formulation: The primary contribution is the reformulation of the eigenstrain evolution into a local constitutive model. This eliminates the need for additional global degrees of freedom, allowing the cohesive phase-field model to be implemented as a "drop-in" replacement for standard phase-field models in existing FE codes.
2. Explicit Strength Surface: The method explicitly prescribes a material strength surface for crack nucleation, decoupling it from the fracture energy release rate ($G_c$). This allows for the independent control of nucleation stress and propagation energy.
3. Consistent Tangent Operators: The derivation of closed-form consistent tangent operators for both smooth and non-smooth strength surfaces ensures numerical stability and quadratic convergence rates in the global solver.
4. Pressure-Dependent Criterion: The introduction of a smooth Drucker-Prager-like criterion enables the modeling of pressure-dependent strengthening under compression, a feature often missing in standard phase-field formulations.

Results and Validation
The framework is validated against three benchmark problems:

• Plate with a Hole (Tension/Compression): The model demonstrates mesh-independent behavior when the element size is sufficiently small relative to the phase-field length scale ($\ell/dx \geq 2$). It correctly captures the transition from brittle to cohesive behavior as $G_c$ increases. Under compression, the model enforces non-penetration and correctly predicts shear band formation, though a competition between horizontal tensile cracking and diagonal shear cracking is observed at low $G_c$.
• Single-Edge Notched Plate (Shear): The model successfully captures mixed-mode fracture. Using the Drucker-Prager-like criterion with low reference strains (high compressive strength), the crack path transitions from shear-driven to tensile-driven, demonstrating the model's ability to handle complex stress states without ad-hoc splitting schemes.
• Dynamic Fracture: In dynamic loading scenarios, the model captures crack branching phenomena naturally. The results show that increasing $G_c$ suppresses crack branching and reduces propagation speed, while the strength surface governs the initial nucleation. The model produces mesh-independent and length-scale-independent results regarding peak load and crack paths.

Significance and Claims
The paper claims that this framework provides a powerful and flexible tool for simulating cohesive fracture processes. By decoupling the strength surface from the fracture energy and the phase-field length scale, the method allows for:

• The prediction of crack nucleation based on a prescribed material strength.
• The control of crack propagation and the brittle-to-cohesive transition solely through the fracture energy release rate ($G_c$).
• The simulation of complex phenomena, such as crack branching and mixed-mode propagation, without requiring ad-hoc criteria.
• Seamless integration into standard finite-element codes without modifying global solvers or element formulations.

The authors emphasize that the length scale $\ell$ acts solely as a numerical regularization parameter for the phase-field distribution width and does not influence the physical material strength or the energy dissipation, provided the mesh is sufficiently refined. The source code is made openly available to facilitate further application and verification.