A crack that lives inside the mesh.

Getting K₁ without touching the crack tip

Engineers don't care about the displacement field for its own sake — they care about the stress-intensity factor \(K_I\), the single number that says how hard the crack tip is being pried open and whether the part will fail. The textbook way to extract it from a simulation is a special contour integral around the tip, which needs additional "branch-function" enrichment there. We deliberately took the cleaner route: the energy release rate.

Griffith's insight is that a crack advances when the structure can release stored elastic energy by growing. The energy released per unit of new crack area is \(G = \mathrm{d}U/\mathrm{d}a\), and for a linear-elastic material it ties directly to the stress-intensity factor:

\[ G = \frac{\mathrm{d}U}{\mathrm{d}a}, \qquad K_I = \sqrt{E'\,G}, \qquad E' = \frac{E}{1-\nu^2}\ \text{(plane strain)} \]

This is a beautiful fit for an engine, because \(U = \tfrac12\,\mathbf{u}^\top \mathbf{K}\,\mathbf{u}\) is just the strain energy of a solved model — something physicsbase already computes. We solve the specimen at two crack lengths, one element shorter and one longer, and take a central difference:

def stress_intensity(W, L, nx, ny, material, a, sigma):
    da = W / nx
    for aa in (a - da, a + da):              # crack one element shorter / longer
        m = XFEMCrack(W, L, nx, ny, material, aa)
        _, U[aa] = m.solve(sigma)            # strain energy of the solved model
    G  = (U[a+da] - U[a-da]) / (2*da)      # energy release rate  dU/da
    return np.sqrt(Eprime(material) * G)     # K_I = sqrt(E' G)

No tip enrichment, no interaction integral, no fragile near-tip stress sampling — just two solves and a finite difference. Because energy is a global quantity, it's far more forgiving than reading the singular stress field right at the tip, which is exactly why this route stays accurate with plain Heaviside enrichment.

What came out

We ran a plate of width \(W=1\), height \(2W\), Young's modulus 210 GPa, \(\nu=0.3\), under unit remote tension. The crack cuts through the middle row of elements; the enriched nodes and the recovered opening look like this:

The number that matters is \(K_I\). The handbook value for this geometry is \(K_I = \sigma\sqrt{\pi a}\,F(a/W)\) with the standard edge-crack polynomial \(F(\alpha) = 1.12 - 0.231\alpha + 10.55\alpha^2 - 21.72\alpha^3 + 30.39\alpha^4\). We compared across three crack depths and two mesh densities:

Two things to notice. First, the answers are close to the handbook across the whole range of crack depths — a few percent, from a method that never once meshed the crack. Second, and more telling, the error halves when the mesh is refined: from ~3.5% at 40×41 to ~1.5% at 80×81. That's the signature of a convergent method — the discretisation error is shrinking toward zero, not sitting at some irreducible modelling bias. The remaining gap is the price of using plain Heaviside enrichment without the extra crack-tip branch functions, and it keeps closing as the mesh gets finer.

Why this one is a good stress test

Topology optimization was one physics in a loop; phase-field fracture was two physics in a loop. XFEM is different again: it's a case where the approximation space itself is modified. The unknowns are no longer just "displacement at nodes" — some of them are jump amplitudes with a sign function baked in, and the quadrature has to respect a discontinuity that ignores element boundaries. An engine that only exposed a fixed catalogue of elements couldn't do this. Because physicsbase exposes the pieces an element is made of — the constitutive matrix, the shape-function derivatives — the enrichment can be layered on top without the engine's cooperation, and the part that is "just an elastic quad" is still the validated code underneath. That's the whole bet of treating simulation as composable primitives rather than a closed list of features.

Honest footnotes

We implemented Heaviside (jump) enrichment only, and placed the crack tip on an element boundary so the tip sits in a single column of nodes — which is why we could skip the near-tip branch-function enrichment that a production XFEM code adds. We recover \(K_I\) from the global energy release rate by finite-differencing the strain energy in the crack length, not from an interaction (M-) integral; the energy route is robust and mesh-convergent but doesn't separate mixed-mode \(K_I\)/\(K_{II}\), and here the loading is pure mode I by symmetry. The formulation is standard and common to the references below; the code is our own and reuses none of the published implementations. The point, as with the other case studies, is narrow: a general-purpose FE engine, handed nothing but its own constitutive matrix and shape functions, reproduces a fracture-mechanics benchmark to within a percent or two — and converges as you refine.

References

1. T. Belytschko & T. Black, "Elastic crack growth in finite elements with minimal remeshing," Int. J. Numer. Methods Eng. 45 (1999), 601–620.
2. N. Moës, J. Dolbow & T. Belytschko, "A finite element method for crack growth without remeshing," Int. J. Numer. Methods Eng. 46 (1999), 131–150.
3. D. Datta, "Introduction to eXtended Finite Element (XFEM) Method" (2013). arXiv:1308.5208
4. A. R. Jimenez et al., "An eXtended Finite Element Method Implementation in COMSOL Multiphysics: Solid Mechanics" (2021). arXiv:2109.03153
5. H. Tada, P. C. Paris & G. R. Irwin, The Stress Analysis of Cracks Handbook, 3rd ed. (ASME Press, 2000) — edge-crack geometry factor.