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

> Explore the phase-field fracture model in PhaseFieldX. Discover its staggered finite-element formulation, spectral energy splits, and modular degradation functions for brittle fracture simulation.

- Repository: [Miguel Castillón/phasefieldx](https://github.com/castillonmiguel/phasefieldx)
- Tags: deep-dive
- Published: 2026-02-26

---

**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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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 `φ`.

```python

# 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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`

```python

# 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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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:

```python

# 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`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field_Fracture/Input.py) and the solver entry point in [`src/phasefieldx/Element/Phase_Field_Fracture/solver/solver.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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.