A crack inside the mesh.

Getting K₁ without touching the crack tip

Engineers do not want the displacement field for its own sake. They want the stress-intensity factor \(K_I\). This one number says how hard the crack tip is opened and whether the part will fail. The textbook method to get it from a simulation is a special contour integral around the tip. That method needs extra "branch-function" enrichment there. We took the simpler route on purpose: the energy release rate.

Griffith showed that a crack advances when the structure can release stored elastic energy as it grows. The energy released per unit of new crack area is \(G = \mathrm{d}U/\mathrm{d}a\). For a linear-elastic material it links 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 fits an engine well, because \(U = \tfrac12\,\mathbf{u}^\top \mathbf{K}\,\mathbf{u}\) is the strain energy of a solved model. physicsbase already computes it. We solve the specimen at two crack lengths, one element shorter and one longer. Then we 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)

There is no tip enrichment, no interaction integral, and no fragile near-tip stress sampling. There are only two solves and a finite difference. Energy is a global quantity. Thus it is more forgiving than a reading of the singular stress field at the tip. This is why this route stays accurate with plain Heaviside enrichment.

What came out

We ran a plate of width \(W=1\) and height \(2W\), with Young's modulus 210 GPa and \(\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 important number 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:

Notice two things. First, the answers are close to the handbook across the whole range of crack depths, within a few percent, from a method that never meshed the crack. Second, and more important, the error halves when the mesh gets finer: from about 3.5% at 40×41 to about 1.5% at 80×81. This is the sign of a convergent method. The discretisation error goes toward zero. It does not stay at a fixed modelling bias. The remaining gap is the cost of plain Heaviside enrichment without the extra crack-tip branch functions. 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. Here the approximation space itself changes. The unknowns are no longer only "displacement at nodes." Some are jump amplitudes with a sign function. The quadrature must respect a discontinuity that ignores element boundaries. physicsbase gives the constitutive matrix and shape-function derivatives, so enrichment can be added on top. The part that is "just an elastic quad" remains numerically verified code; the XFEM-specific evidence must still be assessed separately.

Honest footnotes

We used Heaviside (jump) enrichment only. We placed the crack tip on an element boundary, so the tip is in a single column of nodes. Thus we could skip the near-tip branch-function enrichment that a production XFEM code adds. We get \(K_I\) from the global energy release rate. We finite-difference the strain energy in the crack length. We do not use an interaction (M-) integral. The energy route is robust and mesh-convergent. But it does not separate mixed-mode \(K_I\) and \(K_{II}\). Here the loading is pure mode I by symmetry. The formulation is standard and common to the references below. The code is our own. It reuses none of the published implementations. The point, as in the other case studies, is narrow. A general-purpose FE engine, given only its own constitutive matrix and shape functions, reproduces a fracture-mechanics benchmark to within a percent or two. It also converges as you refine the mesh.

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.