I have been building a deliberately small finite-element demonstrator of cyclic phase-field fracture in a 2D single-edge-notched-tension (SENT-style) specimen. It is a research/teaching demo, not a validated gear model, and it must not be used for engineering allowables or life predictions. The goal is clarity: get the fatigue math and the cycle bookkeeping honest, then compare a fully resolved baseline against an accelerated scheme.
The specimen is a plane 2D rectangular domain with a narrow enforced-damage strip on the left edge that approximates a notch (it is not a physical notch cut). A prescribed displacement is applied and released in closed non-negative load triangles, with eight loading stages per cycle by default. Around one shared dependency-free core, the demo ships several executables: a cycle-resolved matching-law reference (sent_fatigue_reference.py), a constant-load accumulation run (sent_fatigue_cla.py), the historical PhaseFieldX-law baseline kept unchanged (sent_fatigue.py), and a separate coupled u/d Newton/SNESVI robustness experiment (sent_fatigue_monolithic.py). The core fatigue/history math is pure Python and imports no FEniCSx, PETSc, or NumPy; only the optional examples need a real DOLFINx/PETSc stack.
Let the phase field be phi (the code calls the field d), with phi = 0 intact and phi = 1 broken. The quadratic intactness function is
The elastic energy uses a volumetric–deviatoric tension approximation (not a spectral split; its deviatoric contribution stays in the positive part):
The driving fracture history is irreversible:
Fatigue is tracked through a normalized toughness scale
where Gc is the fracture energy and ell the regularization length. An accumulated driving variable alpha_bar degrades the fracture term through a fatigue factor (written F in the code):
The free energy the solver minimizes is
F is a coefficient inside the phase weak form, so a converged update affects the next state. A staggered loop alternates the split elasticity solve and the scalar phase solve, and clamps phi monotonically (phi_n >= phi_{n-1}) in addition to the history construction.
The reference resolves each closed R -> 1 -> R physical cycle and, only after each state converges, accumulates the positive driving increment:
The constant-load approximation (CLA) assumes the mechanics state is converged at maximum load and then represents a fixed number of identical physical cycles without solving the unloading path. With load ratio R:
This is deliberately different from the reference: CLA does not multiply psi+ by g(phi). So at a fixed converged psi+, changing the damage does not change the CLA increment. CLA also cannot be compared pointwise with the per-stage history of the resolved reference; only matched total physical-cycle bounds are meaningful.
To compare fairly, the campaign runner sends identical --cycles, --R, mesh, and displacement bounds to the reference and to two CLA runs (N_jump = 1 and N_jump = 4).
On the remote mapc4_wsl2 container, the matched reference completed 24 cycles. After the DG0 correction described below, the CLA run with N_jump = 1 matched the reference to 0.001% or better on the recorded summaries; the N_jump = 4 run showed a small cycle-jump error. That ordering is the expected one, and it is the point of the comparison rather than a surprise:
N_jump = 1 makes CLA a one-cycle-at-a-time update at maximum load, so with matching bounds it should agree with the closely related reference.I want to be careful here: this is a consistency check between two related laws, not a validation. The N_jump = 4 deviation is a numerical jump artifact, not a physical conclusion.
In some runs the global mean damage flattens. That is worth reading cautiously, because a flat mean damage is not the same as no crack growth. Mean damage averages over the whole domain, so a localized crack or notch contributes only a small, slowly changing global signal — a near-constant mean can hide ongoing local damage. Under displacement control, the specimen may also partially unload as damage localizes, which can suppress the global driving term and flatten the mean. A plateau could therefore mean genuine arrest (the driving energy is no longer enough to advance the front) or simply that the chosen global diagnostic is insufficient. Distinguishing the two needs local fields plus crack-length or reaction-force outputs, which the demo does not yet produce. A plateau is a prompt to look closer, not a result on its own.
The first real DG0 run surfaced a concrete numerical problem: a CG1 L2 projection of the non-negative psi+ can undershoot at nodes, creating spurious negative values and false positive unloading increments. The fix had several parts: put the non-negative quantities (psi+, history, alpha_bar, increment, and F) in DG0, where an L2 projection is a cell average and cannot go negative at nodes; keep the phase field phi itself in continuous CG1 so its gradient is represented in the weak form; apply the non-negativity guard to the internal state before the history/fatigue updates (not as a CSV-only display clamp), and reject non-finite projections. In addition, CLA now performs a documented non-physical zero-load, H = 0 notch-profile equilibration before its physical jumps. CLA starts directly at maximum load, whereas the reference first equilibrates at minimum load; DG0 made that initialization difference significant and CLA initially failed its staggered iteration cap. The setup solve leaves alpha_bar and F untouched and produces no physical cycle/XDMF/CSV record. After those changes the reference and CLA use CG1 damage with DG0 H/F coefficients, and the matched remote run converged again.
Deliberately missing pieces, in rough priority order:
N_jump with an error-indicator-driven adaptive cycle jump.All numbers above are demo/run artifacts on a small structured mesh. They are consistency and convergence observations, not validated fracture results.
FEniCSx, fatigue, fracture, phase-field — Sep 21, 2026