- The paper introduces the Proximal Galerkin method, rigorously enforcing irreversibility and boundedness in phase field fracture models.
- It reformulates the fracture problem as a sequence of saddle-point PDEs using latent variable mapping, seamlessly integrating with Galerkin FE discretization.
- Numerical validations highlight superior energy control and sharp damage localization across tests including dynamic crack branching and impact experiments.
Proximal Galerkin for Phase Field Fracture: A Rigorous Solution Framework
Introduction
The phase-field (PF) method is an established variational framework for modeling brittle fracture, offering flexibility and robust simulation of crack nucleation, branching, and merging through a continuous field. Despite its advantages over discrete crack-tracking methods (e.g., XFEM, adaptive remeshing), PF models face significant numerical challenges, particularly regarding strict enforcement of physical inequality constraints — notably, irreversibility (crack growth cannot reverse) and boundedness of the phase-field variable. Standard approaches to enforcing these constraints, such as penalty and history-functional (HF) methods, introduce ad hoc parameters or only approximate constraint satisfaction. This paper proposes the Proximal Galerkin (PG) methodology as a mathematically rigorous, constraint-preserving framework that reformulates the PF fracture problem as a sequence of saddle-point PDEs using latent variables. PG strictly maintains both irreversibility and bounds, is agnostic to the choice of crack surface regularization (AT1 vs AT2), and integrates seamlessly with Galerkin FE discretization.
The PF model minimizes an energy functional dependent on displacement and phase-field variables, employing:
- Degradation Function: Quadratic g(ϕ)=(1−ϕ)2 reduces the active strain energy.
- Strain Energy Split: Spectral decomposition separates tensile and compressive effects, ensuring fracture is only driven by tensile stresses.
- Crack Surface Density Functional: Regularizes sharp cracks. The AT2 (quadratic) and AT1 (linear) models differ in their geometric function α(ϕ), influencing support and threshold behavior.
- Irreversibility Condition: Enforced through inequality constraints on the evolution of ϕ.
The PF variable is constrained by 0≤ϕ≤1 and monotonicity in time. Conventional penalty and HF approaches enforce these constraints only approximately, with convergence depending on parameter tuning and potential violation at the discrete level.
Proximal Galerkin Framework
PG methodology employs the Bregman proximal point algorithm to recast the constrained minimization as a recursive sequence of unconstrained (yet regularized) minimization problems. The key innovations are:
- Latent Variable Mapping: The constraint is encoded via an invertible transformation between ϕ and latent ξ, with ϕ=1+exp(ξ)ϕprev+exp(ξ). This dynamic update ensures strict enforcement of both bounds and irreversibility.
- Iterative Saddle-Point PDEs: Each PG step solves a coupled nonlinear system for (u,ϕ,ξ), which are readily discretized with standard C0 Galerkin FE spaces. PG requires no singular penalty parameters; proximity parameter βk affects only algorithmic convergence, not solution quality.

Figure 1: Analytical profiles of the phase-field α(ϕ)0 and latent variable α(ϕ)1 for different regularizations compared with the numerical solutions. AT2 is characterized by infinite support, while AT1 exhibits compact support, as seen by α(ϕ)2 at boundary.
PG delivers superior pointwise enforcement of constraints through nodal quadrature, and the saddle-point structure yields efficient Newton-based solution strategies.
Numerical Validation and Comparative Analysis
1D Bar and Crack Surface Regularization
Comparison of analytical and PG-computed solutions for both AT2 and AT1 regularizations demonstrates excellent agreement. AT1 produces sharply localized cracks (compact support) and requires explicit enforcement of lower bounds, while AT2 decays exponentially but lacks sharp cutoff.
Single Edge Notched Tension and Shear Tests
The PG approach strictly enforces bounding constraints:
- Reaction force and crack area evolution under tension closely match reference and penalty-based results with no negative α(ϕ)3 values appearing.
- Under shear loading, PG and HF formulations exhibit qualitatively identical crack propagation and force-displacement curves, though HF shows more diffuse cracks and only approximate irreversibility.



Figure 2: Single edge notched tension test; PF (top) and latent variable (bottom), demonstrating strict constraint enforcement and sharper damage localization for AT1.




Figure 3: Single edge notched shear test; crack progression at three applied displacements for both AT2 (top) and AT1 (bottom) models.
Figure 4: Single edge notched shear test; force-displacement curves for AT2 and AT1, confirming PG and HF performance parity, with PG providing stricter monotonicity.
Dynamic Crack Branching
Dynamic fracture simulations on progressively refined meshes show:







Figure 6: Dynamic crack branching; PF contours at α(ϕ)4 comparing PG, HF, and unconstrained approaches across mesh refinements.
Figure 7: Dynamic crack branching; elastic-strain (top) and dissipated (bottom) energy time history for AT2.

Figure 8: Dynamic crack branching; violation of irreversibility constraint visualized — PG (top) maintains constraints, HF (middle) and unconstrained (bottom) experience localized violations.




Figure 9: Dynamic crack branching; PF contours for AT1 via PG and history functional.
Figure 10: Dynamic crack branching; elastic-strain (top) and dissipated (bottom) energy for AT1.
Figure 11: Dynamic crack branching; latent variable contours for AT2 (top) and AT1 (bottom) — AT1 shows sharper, more localized constraint enforcement.
Kalthoff--Winkler Impact Experiment
The Kalthoff--Winkler impact experiment further validates the methodology:




Figure 13: Kalthoff--Winkler; PF contours at α(ϕ)6 for PG, HF, and unconstrained AT2 models over mesh refinement.
Figure 14: Kalthoff--Winkler; elastic-strain (left) and dissipated (right) energy time history for AT2.

Figure 15: Kalthoff--Winkler; irreversibility violation — PG maintains strict monotonicity, HF and unconstrained permit violations localized on cracks.

Figure 16: Kalthoff--Winkler; PF contours for AT1 showing sharp, localized crack growth in PG formulation.
Figure 17: Kalthoff--Winkler; elastic-strain (left) and dissipated (right) energy time histories for AT1.
Implications and Future Directions
The PG methodology provides a mathematically principled approach to enforcing physical constraints in PF fracture without introducing artificial parameters or suffering from ill-conditioning. The algorithm is robust across regularization choices (AT1/AT2), exhibits efficient convergence (one PG iteration per staggered iteration on average), and integrates seamlessly into existing staggered/non-monolithic PF solvers.
Practically, PG yields more accurate damage localization, better control of dissipated energy, and strict monotonicity, which is essential for high-fidelity fracture modeling in structural, geomechanical, and multiscale applications. Theoretical implications include the elimination of singular parameter regimes and the preservation of variational structure at both continuous and discrete levels. Future extensions may encompass adaptive FE discretizations, higher-order elements, and direct monolithic coupling of multiphysics phenomena.
Conclusion
The Proximal Galerkin methodology rigorously enforces both irreversibility and boundedness in PF fracture models, providing a robust, general-purpose solution framework agnostic to regularization choice and superior to penalty or HF methods in numerical accuracy and constraint enforcement. The PG algorithm's structure-preserving approach, mathematical consistency, and computational efficiency make it highly suitable for advanced fracture mechanics simulations and scalable to broader classes of constrained variational problems (2604.26210).