Research / Fracture and damage

Phase-field fracture of a notched plate

After Nathan Shauer, “Less than 500 lines self-contained Python finite element implementation of the phase-field method for fracture mechanics”, Advances in Engineering Software (2025), doi:10.1016/j.advengsoft.2025.104013.

Against the author’s code
30 of 30 checks pass· see the checks
Against the paper’s figures
Not compared: the full text was not accessible to us

Research, not a design tool. Nothing here is an engineering result.

A square plate with a notch from its left edge to the centre is pulled apart at the top. The AT2 phase-field model represents the crack as a damage field d between 0 (intact) and 1 (broken), smeared over a length ℓ. The experiment pane runs the Stabileo Rust implementation in Stabileo’s viewer; the orange curve is what the author’s Python code computed.

What the paper shows

The article presents a self-contained Python finite element code, under 500 lines, for brittle fracture with the phase-field method, and uses it on three examples: a bar in tension compared with the one-dimensional AT2 solution, a single-edge notched plate in tension (this page), and a three-point bending test. We read the code in full (the article’s text was not accessible to us) and everything below about the method comes from it.

In the replay the damage grows smoothly ahead of the notch while rises linearly to ; at the next increment , after , and the force drops to almost zero. shows the open crack on the deformed plate.

The model

Displacement u and damage d are both bilinear on quadrilaterals, integrated with 2×2 Gauss points; plane strain, small strain, quasi-static loading.

  • Elasticity at fixed d: div(g(d) C : ε(u)) = 0 with g(d) = (1 − d)², applied to the whole stress.
  • Damage (AT2) at fixed u: (Gc/ℓ + 2H) d − Gc ℓ Δd = 2H.
  • Driving energy: H = ψ+ = ½ Σ σi εi over the positive principal strains, with σi the undamaged principal stresses: elasticity is degraded fully while damage is driven by tension only, a “hybrid” formulation.
  • Alternating solution per load step: solve u with the current d, then d with the new u, until the elasticity residual ‖K(d)u − f‖ falls below 10−5.
  • Irreversibility: none in the reference code; damage is recomputed from the current energy and is not bounded by 1. to keep it on unloading (not in the reference).

Units are kN and mm. The homogeneous AT2 peak stress σc = √(27 E Gc / 256 ℓ) ≈ 2.82 kN/mm² here, so ℓ behaves as a strength as much as a regularisation length: and the peak load falls. The mesh matters too: the crack band needs elements a few times smaller than ℓ; .

How we reproduced it

  1. The author’s code (commit 61037c6, MIT licence) was run unmodified for its three examples in Python 3.12.15 with numpy 2.1.3 and scipy 1.14.1; every step’s force, iteration count, residual and damage range, and nodal fields at chosen steps, were saved with a manifest of versions and checksums.
  2. The Rust implementation regenerates the same meshes (bit-identical), applies the same boundary conditions, loop and stopping rule, and solves with the Stabileo engine’s sparse Cholesky factorisation.
  3. 30 named checks compare the two, each with a stated tolerance, in continuous integration. The browser bundle and the replay on this page are checked separately against the native build.

Comparison with the reference

30 validation checks run in CI, each with its tolerance
CheckObservedToleranceResult
ex1: mesh identical to the reference (node and element fingerprints)equal—✓ pass
ex1: same number of load stepsequal—✓ pass
ex1: same number of staggered iterations at every stepequal—✓ pass
ex1: maximum damage at every step8,3 × 10−121 × 10−9✓ pass
ex1: stress at the bar centre at every step1,1 × 10−121 × 10−9✓ pass
ex1: nodal damage at step 504,5 × 10−101 × 10−8✓ pass
ex1: nodal damage at step 733,8 × 10−101 × 10−8✓ pass
ex1: nodal displacement at step 501,6 × 10−101 × 10−8✓ pass
ex1: centre stress vs 1D AT2 homogeneous solution while max damage < 0.257,3 × 10−51 × 10−3✓ pass
ex2: mesh identical to the reference (node and element fingerprints)equal—✓ pass
ex2: same number of load stepsequal—✓ pass
ex2: same number of staggered iterations at every stepequal—✓ pass
ex2: peak force at the same load stepequal—✓ pass
ex2: peak force1,2 × 10−121 × 10−9✓ pass
ex2: force at every step up to the peak1,6 × 10−121 × 10−9✓ pass
ex2: maximum damage at every step up to the peak4,2 × 10−121 × 10−9✓ pass
ex2: force at every step, including fracture (comparison tolerance)2,2 × 10−71 × 10−5✓ pass
ex2: nodal damage at the peak step 565 × 10−111 × 10−8✓ pass
ex2: nodal displacement at the peak step 568,5 × 10−111 × 10−8✓ pass
ex2: nodal damage at step 57, reference or its mirror image2,9 × 10−51 × 10−4✓ pass
ex2: nodal damage at step 78, reference or its mirror image3,1 × 10−101 × 10−8✓ pass
ex3: mesh identical to the reference (node and element fingerprints)equal—✓ pass
ex3: same number of load stepsequal—✓ pass
ex3: same number of staggered iterations at every stepequal—✓ pass
ex3: peak force at the same load stepequal—✓ pass
ex3: force at every step up to the peak4 × 10−111 × 10−9✓ pass
ex3: force at every step, including fracture (comparison tolerance)1 × 10−61 × 10−4✓ pass
ex3: nodal damage at the peak step 566,9 × 10−111 × 10−8✓ pass
ex3: nodal displacement at the peak step 561,1 × 10−101 × 10−8✓ pass
ex3: nodal damage at the last step 98, reference or its mirror image4,1 × 10−101 × 10−8✓ pass

One finding: the plate is symmetric about the notch plane, and the crack runs along the row of elements just above it in one code and just below it in the other. Which row is chosen is decided by floating-point round-off; against the mirror image the final damage fields differ by 3.085970368132962 × 10−10. The thresholds in the table are comparison tolerances, not mathematical error bounds. Displacement fields are compared at example 1 step 50, example 2 step 56, and example 3 step 56.

Step table (79 steps)
stepu, mmF, kNiterationsresidualmax dreference F, kN
00,000100,0134528712,8 × 10−80,00010,01345287
10,000200,0269047511,7 × 10−70,00040,02690475
20,000300,0403544314,3 × 10−70,00090,04035443
30,000400,053800717,9 × 10−70,00160,0538007
40,000500,0672423611,3 × 10−60,00240,06724236
50,000600,0806781911,9 × 10−60,00350,08067819
60,000700,0941069812,6 × 10−60,00480,09410698
70,000800,107527513,4 × 10−60,00630,1075275
80,000900,120938614,4 × 10−60,00790,1209386
90,001000,13433915,5 × 10−60,00980,134339
100,001100,147727616,7 × 10−60,01190,1477276
110,001200,16110318 × 10−60,01420,161103
120,001300,174464219,4 × 10−60,01660,1744642
130,001400,187788723,5 × 10−70,01940,1877887
140,001500,201114324,5 × 10−70,02230,2011143
150,001600,214421825,9 × 10−70,02550,2144218
160,001700,227709827,6 × 10−70,02890,2277098
170,001800,240977129,6 × 10−70,03240,2409771
180,001900,254222621,2 × 10−60,03620,2542226
190,002000,26744521,5 × 10−60,04030,267445
200,002100,280642921,8 × 10−60,04460,2806429
210,002200,293815322,2 × 10−60,04910,2938153
220,002300,306960822,7 × 10−60,05390,3069608
230,002400,320078223,2 × 10−60,05890,3200782
240,002500,333166223,8 × 10−60,06410,3331662
250,002600,346223524,5 × 10−60,06970,3462235
260,002700,359248925,3 × 10−60,07550,3592489
270,002800,37224126,1 × 10−60,08160,372241
280,002900,385198627,1 × 10−60,08800,3851986
290,003000,398120328,3 × 10−60,09470,3981203
300,003100,411004829,6 × 10−60,10170,4110048
310,003200,423841132 × 10−60,10920,4238411
320,003300,436645632,4 × 10−60,11680,4366456
330,003400,449408632,9 × 10−60,12490,4494086
340,003500,462128633,5 × 10−60,13330,4621286
350,003600,474804134,3 × 10−60,14210,4748041
360,003700,487433735,1 × 10−60,15130,4874337
370,003800,500015636,2 × 10−60,16090,5000156
380,003900,512548337,4 × 10−60,17100,5125483
390,004000,5250338,9 × 10−60,18160,52503
400,004100,53745243,3 × 10−60,19290,537452
410,004200,549824944,1 × 10−60,20460,5498249
420,004300,56214145,1 × 10−60,21700,562141
430,004400,57439846,5 × 10−60,22990,574398
440,004500,586593548,2 × 10−60,24370,5865935
450,004600,598718754,2 × 10−60,25840,5987187
460,004700,610781455,5 × 10−60,27390,6107814
470,004800,622773357,4 × 10−60,29030,6227733
480,004900,634690459,9 × 10−60,30800,6346904
490,005000,646520566,8 × 10−60,32720,6465205
500,005100,6582769,7 × 10−60,34790,65827
510,005200,669918178 × 10−60,37080,6699181
520,005300,681459687,6 × 10−60,39620,6814596
530,005400,692876598,5 × 10−60,42500,6928765
540,005500,704134118,1 × 10−60,45880,704134
550,005600,7151725148,7 × 10−60,50060,7151725
560,005700,7257948229,3 × 10−60,56060,7257948
570,005800,00047543893088,3 × 10−61,00810,0004752194
580,005900,000454193924,6 × 10−71,00810,0004541974
590,006000,000433467822,5 × 10−71,00830,0004334679
600,006100,000413999222,3 × 10−71,00850,0004139992
610,006200,000395713222,1 × 10−71,00870,0003957132
620,006300,000378519822 × 10−71,00880,0003785198
630,006400,000362336521,9 × 10−71,00890,0003623365
640,006500,000367766319,6 × 10−61,00910,0003677663
650,006600,000352428319,2 × 10−61,00920,0003524283
660,006700,000337749318,9 × 10−61,00930,0003377493
670,006800,000323894718,5 × 10−61,00940,0003238947
680,006900,000310807618,2 × 10−61,00940,0003108076
690,007000,00029843417,9 × 10−61,00950,000298434
700,007100,000286724817,6 × 10−61,00950,0002867248
710,007200,000275635117,4 × 10−61,00950,0002756351
720,007300,000265123317,1 × 10−61,00960,0002651233
730,007400,000255151616,9 × 10−61,00960,0002551516
740,007500,000245684916,7 × 10−61,00960,0002456849
750,007600,00023669116,4 × 10−61,00960,000236691
760,007700,0002281416,2 × 10−61,00960,00022814
770,007800,000220004516 × 10−61,00960,0002200045
780,007900,00021225915,8 × 10−61,00950,000212259

Comparison with the paper

Not done: the article’s full text was not accessible to us, so no number printed in it has been compared. Figures here come from our own runs.

Limitations

  • Agreement with a reference code is not validation against experiments, and none of this is an engineering certification.
  • The reference mesh refines only a 0.04 × 0.04 mm box around the notch tip; beyond it the crack crosses elements wider than 3ℓ.
  • No bound constraints: damage can exceed 1 slightly (the reference reaches 1.0096 here and 1.0021897551918642 in example 1) and undershoot 0. The recorded Python reference outputs have minimum nodal damage 3.076336950017797e-08 (example 1), -0.0031224948220514573 (example 2), and -0.06024191772700118 (example 3). This also occurs in the reference code’s unconstrained AT2 discretization. Negative damage is unphysical: (1-d)² at the example 3 minimum is about 1.124, above intact stiffness. Assembly uses interpolated damage at quadrature points, so this is not a 12.4% assembled-stiffness error. The solver does not clamp damage; the replay’s quantized nodal field and colour scale clamp it for display.
  • The alternating scheme is not accelerated; at unstable crack growth it can take hundreds of iterations in one step.
  • Browser runs are limited in mesh size and parameter range.

Code and references

Cite this

Stabileo Research, “Phase-field fracture of a notched plate”, reproduction of Nathan Shauer (2025), doi:10.1016/j.advengsoft.2025.104013. Implementation: research/fracture/shauer-2025 at commit 0ce459e9eb75. Engine: the Stabileo engine at the same commit. Permanent links to this exact build are planned.