Source-linked AI summary
Phase field modeling and computer implementation: A review
X. Zhuang, S. Zhou, G. D. Huynh, P. Areias, T. Rabczuk
TL;DR
Fracture modeling requires approaches that handle complex crack processes without burdensome crack representation and tracking. This paper reviews phase field fracture theories and their computer implementation, alongside applications and benchmark examples. It presents phase fields as a scalar-field approach that determines fracture evolution through an additional differential equation while avoiding explicit crack-surface tracking and extra fracture criteria.
Problem
Fracture methods must address complex processes such as initiation, propagation, coalescence, branching, and crack interactions, which burden discontinuous approaches with topology tracking and additional criteria.
Method
The paper synthesizes phase field fracture theories, discretization methods, finite-element implementation, applications, and representative numerical examples.
Results
Phase field models represent discrete cracks with an additional scalar field whose evolution determines fracture shape and propagation without directly introducing physical discontinuities or tracking the fracture surface.
Takeaways & Limitations
The reviewed phase field framework is presented as practicable for complex fracture problems, including brittle, ductile, multi-field, quasi-static, and dynamic settings.
Takeaways & Limitations
Although phase field models can couple with several discretization methods, most models use finite-element discretization.
Abstract
from arXiv · showhide
This paper presents an overview of the theories and computer implementation aspects of phase field models (PFM) of fracture. The advantage of PFM over discontinuous approaches to fracture is that PFM can elegantly simulate complicated fracture processes including fracture initiation, propagation, coalescence, and branching by using only a scalar field, the phase field. In addition, fracture is a natural outcome of the simulation and obtained through the solution of an additional differential equation related to the phase field. No extra fracture criteria are needed and an explicit representation of a crack surface as well as complex track crack procedures are avoided in PFM for fracture, which in turn dramatically facilitates the implementation. The PFM is thermodynamically consistent and can be easily extended to multi-physics problem by 'changing' the energy functional accordingly. Besides an overview of different PFMs, we also present comparative numerical benchmark examples to show the capability of PFMs.
1 Introduction
The introduction motivates phase field models as continuous alternatives for fracture and outlines this review's coverage of theory, discretization, implementation, applications, and numerical examples.
- Motivation: Discontinuous fracture methods require crack-topology representation and tracking, with complex branching and interactions often needing additional criteria.Examples include triangular facets, level sets, and associated crack-tracking algorithms.
- Motivation: Phase field models avoid displacement discontinuities by smearing fracture over a finite-width localization band with an intrinsic length scale.They are presented alongside gradient-damage and screened-Poisson models as continuous fracture approaches.
- Motivation: PFMs naturally simulate crack initiation, propagation, coalescence, and branching on a fixed mesh without laborious crack-surface tracking.Avoiding crack-surface tracking is especially useful for three-dimensional simulations.
- Scope of the review: The review covers advances in phase-field theories, discretization methods, finite-element implementation, applications, and representative numerical examples.The examples span two- and three-dimensional settings and demonstrate practicability and capability.
2 Theories of phase field models for fracture
Phase-field fracture models represent cracks with scalar fields governed by energy-based or dissipative evolution equations. The reviewed formulations span dynamic and quasi-static fracture, incorporate tension–compression treatments, and support numerical regularization and implementation.
- Physics-based phase-field models: The Aranson model describes dynamic brittle fracture through coupled elastodynamic and phase-field equations, with damping and material parameters controlling the response.Its displacement equation includes density and viscous damping, while the order-parameter equation includes diffusion, model parameters, and displacement-gradient coupling.
- Physics-based phase-field models: Physics-based models use an order parameter that distinguishes intact material from fully cracked regions and evolves through dissipative dynamics derived from a free-energy functional.In the Aranson model, s = 1 outside cracks and s = 0 inside them; the phase evolution follows the variation of the free energy.
- Physics-based phase-field models: The reviewed physics-based models reproduce crack initiation, propagation, branching, instability, sound emission, fragmentation, and experimentally observed oscillatory behaviors, though discrepancies remain for finite-width Mode-I propagation.The Aranson model captures multiple fracture behaviors but differs from experiments for Mode-I cracks in rectangular strips; the Henry–Levine model reproduces branching, oscillations, and supercritical Hopf bifurcation.
- Variational and mechanics-based models: Variational fracture models formulate crack initiation, propagation, and branching as energy minimization, then regularize the discontinuous crack set with a scalar phase field.The scalar field transitions continuously from s = 1 in intact material to s = 0 in fully damaged material; the regularized formulation approaches the original one as ϵ → 0 in the sense of Γ-convergence.
- Variational and mechanics-based models: Ginzburg–Landau evolution equations govern phase-field crack propagation under irreversibility constraints, with mobility controlling dissipation and the infinite-mobility limit corresponding to quasi-static propagation.The reviewed formulation uses a mobility parameter M ≥ 0; finite M gives viscous quasi-static behavior, while M → ∞ yields quasi-static crack propagation.
2.3 Phase field approximation of the sharp crack topology
Phase-field fracture replaces the sharp crack topology with a regularized crack surface density function built from a scalar phase field and its gradient. Different models use distinct geometric functions and potentials, with the double-well choice providing an energy barrier that supports fracture irreversibility.
- The regularized functional approximates the crack surface area through a crack surface density function composed of the phase field and its spatial gradient.The density function approximates the Dirac-delta distribution along the crack.
- Kuhn et al. proposed a generic phase-field fracture energy containing an energetic degradation function, a local fracture-energy function, and a normalization constant.The normalization ensures convergence of the regularized fracture contribution to the crack surface measure as the internal length tends to zero.
- The local fracture-energy function includes double-well and monotonous families, with the convex quadratic case β = −1 commonly used.The reported families are w(s) = 16s2(1 −s2) and w(s) = (1 + βs)(1 −s), β ∈[−1, 1].
- The double-well potential supplies an energy barrier between broken and undamaged states, whereas w(s) = (1 −s)2 requires additional irreversibility to prevent crack healing.The distinction is illustrated by the different degradation functions shown in Fig. 1.
- Wu’s generic crack surface density function uses a geometric function, an internal length scale, and a scaling parameter to regularize and recover the sharp crack surface in the zero-length limit.The geometric function characterizes homogeneous phase-field evolution, while the scaling parameter recovers the crack area as b →0.
2.4 Energetic degradation function
Energetic degradation functions control how elastic energy and material stiffness decrease as the crack phase field evolves. The review summarizes admissibility conditions and several classical and parametric function families intended to improve fracture predictions.
- The degradation function models elastic-energy degradation and stiffness loss as the phase field approaches the broken state.It must increase monotonically, satisfy g(1) = 1 and g(0) = 0, and have g′(0) = 0.
- Classical choices include the quadratic function g(s) = s2 and the quartic function g(s) = 4s3 −3s4.Kuhn et al. also adopted a cubic degradation function.
- Wu’s unified damage model uses a monotonically decreasing energetic function v(ϕ) with v(0) = 1, v(1) = 0, and v′(1) = 0.The function describes degradation of the initial energy during crack phase-field evolution.
- Wu parameterizes the degradation function using a positive exponent and a continuous positive function Q(ϕ), which can be represented by a polynomial expansion.The polynomial is expressed as Q(ϕ) = α1ϕP(ϕ).
- Sargado et al. introduced a three-parameter family intended to improve critical-load predictions while preserving the bulk material’s linear-elastic response before fracture.Their numerical examples indicated superiority over the classical quadratic degradation function.
2.5 On constitutive assumptions
Constitutive assumptions determine how phase-field models degrade stress and prevent nonphysical crack behavior. Isotropic formulations are simpler and cheaper, while anisotropic and hybrid formulations address compression and crack-face interpenetration with differing implementation costs.
- The Amor constitutive model produces transverse isotropy about the crack normal, with stiffness degradation governed by the Lamé parameters and phase field.Its formulation permits crack evolution from positive volume changes and shape distortion.
- Volumetric-deviatoric and spectral energy splits preserve the variational character of phase-field fracture models.They also support the energy contribution to phase-field evolution and the non-interpenetration condition on crack surfaces.
- Isotropic and anisotropic formulations: Isotropic formulations have linear stress-strain relations and lower implementation cost, but they allow fracture in compression and crack-face interpenetration.Anisotropic formulations avoid these drawbacks through energy splitting but require laborious nonlinear implementation.
- Hybrid formulation: Hybrid formulations retain the isotropic model’s linear momentum balance while using anisotropic phase-field evolution to avoid fractures in compression.This design retains favorable computational cost while addressing a principal isotropic-model drawback.
- Zhang et al.’s separate mode-I and mode-II driving-force parameters cannot practically distinguish tension fractures from shear fractures.The review identifies this as a limitation of the formulation’s constitutive interpretation.
3 PFMs coupled with different discretization methods
Finite elements are the main computational framework for phase-field fracture, but alternative discretizations address higher-order equations and implementation demands. Isogeometric, meshfree, neural-network, and FFT approaches provide different routes for solving these models.
- Finite element method: Finite element methods commonly solve phase-field fracture equations because the governing partial differential equations can be discretized over the spatial domain.Suitable decoupling technology can be used for fully coupled fracture problems.
- Finite element method: High-order phase-field equations remain difficult to solve with finite elements, motivating coupling with isogeometric analysis and meshfree technologies.This difficulty is explicitly identified as a limitation of FEM implementations.
- Isogeometric methods: Fourth-order phase-field models have been implemented with isogeometric analysis, while isogeometric collocation can speed computation by reducing point evaluations.A hybrid collocation-Galerkin formulation is recommended for consistently enforcing Neumann conditions.
- Meshfree methods: Meshfree local maximum entropy approximants can directly solve fourth-order phase-field equations without splitting them into two second-order equations.Their higher-order continuity enables this direct treatment.
- Machine-learning methods: Physics-informed neural networks require few lines of code and may provide computational savings over FEM after training.The reported savings are conditional on the network being trained.
- FFT methods: FFT solvers use staggered updates to solve fracture and mechanical problems separately while retaining simple mesh generation and parallel implementation.FFT-based methods have also been extended to higher-order and multiphase-field fracture.
4 Finite element implementation of phase field methods
The section explains finite-element discretization, coupled solution schemes, and implementations of phase-field fracture models in commercial and open-source software.
- Finite-element discretization: Finite-element discretization represents displacement and phase fields through nodal values, shape-function matrices, and their gradients.The strain and phase-field gradients are obtained using the Bu and Bϕ derivative matrices.
- Finite-element discretization: Dynamic governing equations yield discretized inertial, internal, and external force terms for displacement and phase-field equations.For quasi-static problems, the initial inertial term vanishes.
- Solution schemes: Monolithic schemes solve displacement and phase fields simultaneously, whereas staggered schemes solve displacement first at fixed phase field and then update the phase field.Both implicit and explicit time schemes can be used in staggered formulations; monolithic dynamic calculations can use HHT integration.
- Software implementation: Commercial and open-source implementations include Abaqus, FEniCS, and COMSOL, using subroutines or automated partial-differential-equation solution frameworks.Abaqus implementations use UEL, UMAT, or VUMAT subroutines, while FEniCS supports separate field solves.
- Software implementation: COMSOL implementations use storage and solid-mechanics modules to update history strain, solve the phase field, and modify the elasticity matrix iteratively.The reviewed COMSOL implementations use staggered schemes, with Anderson acceleration applied to improve convergence.
- Element shape functions: Special exponential shape functions can accurately predict surface energy with lower element refinement, but require prior fracture-direction information.The reported forms are available for one-dimensional two-node elements and extend to two dimensions through tensor products.
5 Extensions and applications of the PFMs
The section surveys extensions of phase-field fracture models to ductile, dynamic, finite-deformation, shell, thermal, porous, geological, and length-scale-insensitive settings.
- Ductile fracture: Ductile-fracture extensions couple phase-field models with elasto-plasticity, including formulations for large strains, finite deformation, and elasto-viscoplastic materials.Reported developments modify degradation functions, plasticity formulations, or fracture energies to represent ductile behavior.
- Dynamic fracture: Dynamic phase-field formulations incorporate kinetic energy and support monolithic or staggered time integration for brittle and cohesive fracture.Several approaches use explicit integration, sub-stepping, or adaptive residual-based updates.
- Finite deformation: Finite-deformation formulations modify elastic energies, constitutive relations, balance equations, and phase-field evolution while retaining thermodynamic consistency.Approaches include Neo-Hookean materials, multiplicative principal-stretch decompositions, and energy-momentum-consistent integration.
- Thin structures: Phase-field methods for thin structures address consistency with shell kinematics and distinguish isotropic, anisotropic, and hybrid formulations.Published approaches vary in shell models, phase-field models, and discretization methods.
- Thermal and porous fracture: Thermo-elasto-plastic phase-field frameworks can model temperature-dependent void growth followed by macroscopic crack initiation and propagation.One framework combines modified GTN-type plasticity with phase-field fracture under large deformations.
- Length-scale-insensitive fracture: A length-scale-insensitive cohesive phase-field model makes failure strength, traction-separation behavior, and global fracture responses independent of the internal length scale.The model uses a linear-softening cohesive-zone formulation with optimal characteristic functions.
6 Representative numerical examples
Representative benchmarks demonstrate that phase field models reproduce diverse fracture patterns under tension, shear, mixed loading, and notched geometries, while model choice affects crack evolution and load response.
- The review presents representative numerical examples to demonstrate the capability of phase field fracture modeling.
- 2D notched square plate subjected to tension: A notched square plate under tension develops a horizontal crack through the plate’s middle.The benchmark uses E = 210 GPa, ν = 0.3, and Gc = 2700 J/m2.
- 2D notched square plate subjected to tension: All four tested PFMs obtain the same final tensile fracture pattern, but differ in initiation, propagation, and load-displacement response.Amor et al. (2009) predicts a much higher peak load than the other methods.
- 2D notched square plate subjected to shear: Under shear, Miehe et al. (2010b) and Ambati et al. (2015b) produce inclined fractures, whereas the isotropic and Amor models produce horizontal pure mode II fractures.The former models also allow the plate to sustain a larger shear load than the latter two.
- 2D notched square plate subjected to tension and shear: In the mixed tension-shear benchmark, fracture patterns are evaluated across critical energy release rates Gc = 25, 50, 75 and 100 J/m2.The plate contains two horizontal notches and is 200 mm long, 200 mm high, and 50 mm thick.
- 2D notched semi-circular bend test: In the NSCB test, fracture initiates at the upper tip of the pre-existing notch and propagates vertically, agreeing with experimental observations.
5. Propagation of multiple echelon flaws
For a square plate containing nine echelon flaws under tension, flaw interaction produces a final fracture pattern dominated by fractures from the bottom flaws.
- Nine 45° echelon flaws with varying lengths and spacings interact under tension, and fractures from the bottom flaws dominate the final pattern.The passage attributes this dominance to stress shielding and amplification effects from flaw interaction.
6. Propagation and coalescence of twenty parallel flaws
An anisotropic phase field model predicts symmetric fracture development in a tensile square plate with twenty parallel flaws, with cracking limited to the upper and bottom flaws.
- The anisotropic phase field model produces symmetric fractures in a square plate containing twenty parallel flaws under tension.
- Fractures initiate only from the upper and bottom flaws, with no interior fractures reported.
- In Brazilian discs with three vertically arranged notches, outer cracks propagate toward the disc ends and inner-crack coalescence depends on notch spacing.At S = 1 cm, one inner crack occurs; at S = 3 cm, two inner cracks coalesce.
7. Propagation and coalescence of three parallel flaws in Brazilian discs
Brazilian-disc simulations with three horizontal notches show that notch spacing changes whether only outer cracks or additional inner cracks develop.
- For three horizontal notches, S = 2 cm and S = 3 cm produce similar patterns with only two outer cracks propagating toward the specimen ends.
- At S = 1 cm, additional inner cracks initiate from notch tips and evolve between adjoining notches.
- Phase field modeling is also applied to compressive-shear fracture in rocks, a fracture mode that can form during rock loading.
8. Compressive-shear fracture
The review demonstrates phase-field fracture across static, dynamic, hydraulic, and thin-structure examples, including crack branching, propagation, and interaction. Results show sensitivity to fracture parameters while reproducing crack paths and load responses reported in previous studies.
- Dynamic fracture: Smaller Gc produces more complex Kalthoff crack patterns, including branching at Gc = 5 × 10^3 and 1 × 10^4 J/m^2.For larger Gc, simulations show a single crack; increasing Gc also decreases maximum crack-tip velocity and delays crack initiation.
- Dynamic fracture: Multiple crack branching occurs under dynamic tension for different Gc, while larger Gc produces a larger post-branching angle and lower maximum crack-tip velocity.
- Hydraulic fracture: Injected fluid drives a single pre-existing crack horizontally, while the pressure-concentrated region extends beyond the cracked region because of radial fluid penetration.The fracture pattern and pressure field are reported at t = 8.8 s.
- Hydraulic fracture: For two parallel hydraulic cracks, propagation spacing increases with time and the pressure field follows the phase field with a larger transition band.
- Hydraulic fracture: Three-crack hydraulic simulations show that only the left and right cracks propagate, while the middle crack remains arrested.
- Hydraulic fracture: Two parallel penny-shaped cracks in 3D grow beyond their initial 0.5 m radius and become bowl-shaped during propagation.The radii exceed 0.5 m at t = 13.2 s, with bowl-shaped cracks observed at t = 13.3 s.
- Thin structures: Thin-structure benchmarks reproduce previously reported crack paths and comparable load-deflection curves for solid and Kirchhoff-Love shell elements.The examples include central crack initiation followed by branching toward corners and straight axial crack propagation.
7 Conclusions
The review summarizes phase-field fracture theory, discretization, implementation, and applications. It concludes that a scalar phase field enables complex crack evolution without direct crack-surface tracking, while finite elements remain the predominant discretization.
- Phase-field models represent discrete cracks with an additional scalar field whose evolution determines fracture shape and propagation.
- The models simulate crack initiation, propagation, branching, and merging in arbitrary 2D and 3D geometries without algorithmically tracking fracture surfaces.
- Although isogeometric and mesh-free discretizations are possible, most phase-field fracture models use finite element discretization.
- The review presents finite-element computer-implementation details and representative quasi-static, dynamic, and 2D and 3D examples demonstrating phase-field modeling capability.