How the Phase-Field Fracture Model Works in PhaseFieldX: A Technical Deep Dive

PhaseFieldX implements a staggered finite-element formulation that alternately solves mechanical equilibrium and phase-field evolution until convergence, using spectral energy splits and modular degradation functions to simulate brittle fracture.

The phase-field fracture model in PhaseFieldX provides a robust framework for simulating crack initiation and propagation without explicit tracking of discontinuities. This open-source library, built on FEniCSx/DOLFINx, couples a mechanical displacement problem with a phase-field (damage) problem through an efficient staggered algorithm. The following sections dissect the core components, governing equations, and implementation details found in the castillonmiguel/phasefieldx repository.

Core Components of the Phase-Field Fracture Model

Linear Elastic Material Model

The foundation of the fracture model rests on isotropic linear elasticity defined in src/phasefieldx/Materials/elastic_isotropic.py. This module computes the strain-energy density ψ and the strain tensor ε based on Young's modulus E and Poisson's ratio ν. The elastic constitutive law feeds directly into the energy split routines that determine how damage evolves under tensile versus compressive loading.

Spectral Energy Decomposition

To ensure that cracks only propagate under tension while preserving compressive stiffness, PhaseFieldX employs a spectral decomposition of the elastic energy. The file src/phasefieldx/Element/Phase_Field_Fracture/split_energy_stress_tangent_functions.py defines the functions psi_a and psi_b, which separate the strain energy into a tensile part ψₐ and a compressive part ψ_b. Corresponding stress functions sigma_a and sigma_b compute the degraded and undegraded stress contributions, respectively.

Degradation Functions

The degradation function g(φ) reduces the tensile stiffness as the phase-field variable φ ∈ [0,1] evolves from 0 (intact) to 1 (fully broken). PhaseFieldX offers multiple functional forms—including quadratic, Borden, Alessi, and Sargado—implemented in src/phasefieldx/Element/Phase_Field_Fracture/g_degradation_functions.py. The selector function g(phi, degradation_type) and its derivative dg allow users to switch models by changing the degradation_function parameter in the Input class.

Fracture Surface Energy

The regularized crack surface energy follows the Ambrosio-Tortorelli (AT1/AT2) functionals. The module src/phasefieldx/Element/Phase_Field/energy.py provides calculate_crack_surface_energy, which computes the fracture energy term Gc * ( (1/(c₀ l) φ + (l/c₀) |∇φ|² ). The constant c₀ depends on the selected crack density functional, distinguishing between AT1 (linear) and AT2 (quadratic) models.

The Staggered Solution Algorithm

PhaseFieldX solves the coupled displacement-damage system using a staggered scheme rather than a fully monolithic approach. This algorithm, implemented in src/phasefieldx/Element/Phase_Field_Fracture/solver/solver.py, alternates between mechanical and phase-field sub-problems until convergence.

The staggered loop follows these steps:

1. Solve mechanical equilibrium using the current phase-field φ to degrade the tensile stiffness.
2. Update the history variable H to enforce irreversibility (damage cannot heal).
3. Solve the phase-field evolution using the updated history variable.
4. Check convergence by evaluating the L²-norm change of both displacement u and phase-field φ.

# Pseudocode extracted from solver.py (lines 29-49, 84-106)

while t < final_time:
    # Outer time loop

    while not converged and stagger_iter < max_iter:
        # 1) Solve mechanical problem (Newton → problem_u.solve())

        #    update displacement field u_new → u_old

        
        # 2) Update history variable H (irreversibility)

        project(psi_a(u_new, Data), V_c)          # line 55

        H.x.array[:] = np.maximum(V_c.x.array, V_n.x.array) \
                         if Data.irreversibility == "miehe" else V_c.x.array
        
        # 3) Solve phase-field problem (Newton → problem_phi.solve())

        #    update phase field Φ_new → Φ_old

        
        # 4) Optional fatigue update (lines 62-77)

        
        # 5) Check L²-norm errors (eval_error_L2_normalized)

        #    break when both < stagger_error_tol

The convergence check uses eval_error_L2_normalized to compute normalized L² norms, ensuring that both fields reach equilibrium before advancing to the next time step.

Governing Equations

The weak forms implemented in PhaseFieldX follow the standard phase-field fracture theory with spectral splits.

Mechanical equilibrium (weak form):

∫_Ω [(g(φ) + k) σ_a(u) : ε(δu) + σ_b(u) : ε(δu)] dx = ∫_Γ_N t̄ · δu ds

Here, g(φ) is the degradation function, k is a small stabilization parameter (Data.k), σ_a is the degraded tensile stress, and σ_b is the undegraded compressive stress.

Phase-field evolution (variational inequality):

∫_Ω [g'(φ) (σ_a(u) : ε(u)) δφ + (G_c / (c₀ l)) δφ + G_c l ∇φ · ∇δφ] dx = 0

When fatigue is enabled, an additional term Fatigue * Gc * (...) appears in the phase-field residual, as implemented in lines 17-21 of solver.py.

Optional Extensions

Fatigue Degradation

For cyclic loading scenarios, PhaseFieldX includes a fatigue extension in src/phasefieldx/Element/Phase_Field_Fracture/fatigue_degradation_functions.py. This introduces a history-dependent degradation factor that evolves with a cumulative fatigue variable α, allowing the simulation of subcritical crack growth under repeated loading cycles.

Energy Monitoring and Post-Processing

Accurate energy tracking is essential for verifying the thermodynamic consistency of fracture simulations. PhaseFieldX provides dedicated utilities in src/phasefieldx/Element/Phase_Field_Fracture/energy.py.

The function compute_total_energies calculates:

• Degraded elastic energy using the current phase-field
• Split components ψₐ (tensile) and ψ_b (compressive)
• Total elastic energy E = degraded_energy + ψ_b

# compute_total_energies (energy.py, lines 13-17)

psi_a_val, psi_b_val = compute_elastic_energy_components(u, Data, comm, dx=dx)
degraded_energy = compute_degraded_elastic_energy(u, phi, Data, comm, dx=dx)
E = degraded_energy + psi_b_val

Additionally, calculate_crack_surface_energy from src/phasefieldx/Element/Phase_Field/energy.py computes the regularized crack surface energy γ, its phase-field component γ_φ, and its gradient component γ_∇φ. The solver writes these values to total.energy at each time step, enabling post-processing of energy release rates and crack resistance curves.

Practical Implementation Example

The following minimal driver script demonstrates how to set up and run a phase-field fracture simulation using the Input class and staggered solver:

# example_driver.py

from phasefieldx.Element.Phase_Field_Fracture.Input import Input
from phasefieldx.Element.Phase_Field_Fracture.solver.solver import solve
from phasefieldx.Mesh.generator import create_mesh   # hypothetical helper

from phasefieldx.FunctionSpaces import create_function_spaces   # hypothetical helper

# ------------------------------------------------------------------

# 1. Define simulation parameters

# ------------------------------------------------------------------

data = Input(
    E=210.0,               # Young's modulus [GPa]

    nu=0.3,                # Poisson ratio

    Gc=2.7e-3,             # Critical fracture energy [N/mm]

    l=0.015,               # Length-scale parameter

    degradation="isotropic",
    split_energy="spectral",
    degradation_function="quadratic",
    irreversibility="miehe",
    fatigue=False,
    save_solution_vtu=True,
)

# ------------------------------------------------------------------

# 2. Build mesh and function spaces (replace with your own utilities)

# ------------------------------------------------------------------

msh = create_mesh()                     # returns dolfinx.mesh

V_u, V_phi = create_function_spaces(msh)

# ------------------------------------------------------------------

# 3. Boundary conditions (simple example)

# ------------------------------------------------------------------

bc_u = []      # list of dolfinx.DirichletBC for displacement

bc_phi = []    # list of DirichletBC for phase field

# ------------------------------------------------------------------

# 4. Run the staggered solver

# ------------------------------------------------------------------

solve(
    Data=data,
    msh=msh,
    final_time=1.0,
    V_u=V_u,
    V_Φ=V_phi,
    bc_list_u=bc_u,
    bc_list_phi=bc_phi,
    dt=0.01,
)

Key implementation details include the Input class definition in src/phasefieldx/Element/Phase_Field_Fracture/Input.py and the solver entry point in src/phasefieldx/Element/Phase_Field_Fracture/solver/solver.py (lines 29-49).

Summary

• PhaseFieldX implements a staggered finite-element formulation that alternates between mechanical equilibrium and phase-field evolution until convergence.
• The spectral energy split in split_energy_stress_tangent_functions.py ensures damage only affects tensile strains, preventing unphysical cracking under compression.
• Modular degradation functions in g_degradation_functions.py support quadratic, Borden, Alessi, and Sargado models via the degradation_function parameter.
• The history variable H enforces irreversibility through the Miehe approach or allows healing when disabled.
• Energy monitoring via compute_total_energies and calculate_crack_surface_energy provides thermodynamic consistency checks and post-processing capabilities.

Frequently Asked Questions

What is the staggered solution algorithm in PhaseFieldX?

The staggered algorithm decouples the displacement and phase-field problems to improve robustness and reduce computational cost per iteration. At each time step, the solver first computes mechanical equilibrium with a fixed damage field, then updates the history variable to enforce irreversibility, and finally solves the phase-field evolution using the new strain history. Convergence is assessed via normalized L²-norms of both fields using eval_error_L2_normalized in solver.py.

How does PhaseFieldX prevent damage under compressive loading?

The library employs a spectral decomposition of the elastic energy tensor, implemented in split_energy_stress_tangent_functions.py. This split separates the strain energy into positive (tensile) and negative (compressive) components. The degradation function g(φ) only multiplies the tensile part ψₐ, while the compressive part ψ_b remains undegraded, ensuring that cracks do not propagate under pure compression.

What degradation function options are available in PhaseFieldX?

PhaseFieldX provides four degradation models selectable via the degradation_function parameter in the Input class: quadratic (standard polynomial), Borden (modified for better crack surface approximation), Alessi, and Sargado. These functions and their derivatives dg are defined in src/phasefieldx/Element/Phase_Field_Fracture/g_degradation_functions.py, allowing users to switch models without modifying the solver logic.

How is the critical fracture energy Gc incorporated into the model?

The critical fracture energy Gc appears in the regularized crack surface energy term computed in src/phasefieldx/Element/Phase_Field/energy.py. The functional form depends on the selected AT1 or AT2 model, scaling the phase-field gradient and the field itself by Gc/(c₀l) and Gc*l/c₀, respectively, where l is the length-scale parameter and c₀ is a model-specific constant. This term drives the evolution of the damage field when the elastic energy release exceeds the fracture resistance.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}