Growing a crack from scratch.

Making a material break on a computer is famously hard — cracks are discontinuities, and finite elements hate discontinuities. The phase-field method sidesteps that by turning the crack into a smooth field. We reproduced the benchmark on physicsbase, and it took two engines talking to each other.

Ask a structural engineer where a part will fail and they can usually point to it. Ask a computer to show the crack forming and propagating, and you've walked into one of the hardest problems in computational mechanics. A crack is a moving discontinuity in the displacement field; classic finite elements assume the field is smooth, so for decades cracking meant re-meshing around the tip at every step, or enriching elements with special crack functions. Then came a beautifully sneaky idea: don't track the crack as a surface at all. Represent it as a smooth phase field d that goes from 0 (intact) to 1 (fully broken) over a thin band, and let a PDE decide where it grows. That's the phase-field method, and it's now the dominant way to simulate brittle fracture.

It's also a perfect showcase for physicsbase, because it isn't one physics — it's two, coupled. The displacement field needs a structural solve; the crack field needs a field solve; and they hand results back and forth until they agree. So this is a story about the boundary: what each solver sends the other.

The two fields, and how they talk

The benchmark is the single-edge-notched tension test (Bourdin et al. 2000; Miehe et al. 2010): a unit square with a horizontal pre-notch cut halfway in, clamped at the bottom, pulled straight up at the top. Two unknowns live on the same mesh:

• Displacement u — a structural problem. Cracked material carries no load, so each element's stiffness is scaled by \(g(d) = (1-d)^2\): where the crack field is high, the element goes soft.
• Crack d — a field problem. It solves a screened-Poisson equation driven by how much elastic energy the material is storing. Highly-strained material wants to crack.

Neither can be solved alone — the stiffness depends on d, and d depends on the strain from u. So the solver alternates (a "staggered" scheme), and the interesting thing is exactly what crosses between them each step.

Exchange 1 — the structural solve

The crack solver hands the structural solver a degraded structure: same mesh, but each element's effective stiffness reduced by its current damage. physicsbase solves for the displacement and reports the strain energy in every element — the fuel the crack feeds on.

# in:  the current design + how hard the top is pulled
elements: quad4 mesh, each stiffness scaled by  g = (1 - d)**2
supports: bottom edge clamped;  top edge  uy = Δ  (this load step)

# out: physicsbase returns the deformation and per-element stress
displacements: [ {"ux":…, "uy":…}, … ]
elements:      [ {"stress":[sxx,syy,txy], "strain":[…]}, … ]
reaction (top): 0.699          # the force the specimen is carrying

Exchange 2 — the crack solve

From those stresses the driver computes a crack-driving energy H for each element (kept monotonic, so cracks never heal), and hands that field to the field solver. physicsbase returns the updated crack field — the geometry of the damage.

# in:  a screened-Poisson PDE assembled from physicsbase field operators
#   ( Gc·l ∇²  +  (Gc/l + 2H) )  d  =  2H
conductivity: Gc·l        # crack surface energy × length scale
reaction:     Gc/l + 2H   # H = the driving energy from Exchange 1
source:       2H

# out: physicsbase returns the crack field, clamped irreversible
d: [ 0.02, 0.05, …, 1.00, 1.00, … ]   # 0 intact → 1 broken

Using physicsbase: two kernels, four operators

Written as coupled PDEs this looks like a research code. Written as physicsbase calls, the entire physics is four operator evaluations — two from the structural element, three from the field element (one shared idea). Here is every piece the engine provides; everything else is glue.

The elasticity side — from Quad4

Two calls. The first is the element stiffness, degraded per element by the crack; the second is the element's strain-displacement operator, which is how we read the elastic energy that drives cracking straight out of the displacement solution:

from femengine.elements import Quad4
from femengine.materials import Material

sq = [[0,0],[h,0],[h,h],[0,h]]                          # one h×h cell
KE = Quad4.stiffness(sq, Material(E, nu), {"kind":"plane_strain"})  # (8,8)
B, _ = Quad4._B_at(sq, 0.0, 0.0)                          # (3,8) strain-disp at centre
D  = Material(E, nu).plane_strain_D()                      # (3,3) constitutive

# per element, from the solved displacements u_e:
eps = B @ u_e                 # strain  [εxx, εyy, γxy]
psi = 0.5 * eps @ (D @ eps)     # elastic energy density -> crack driving force

Nothing here is bespoke: KE, B, and D are the same objects the engine uses for an ordinary stress analysis. We just scale KE by \(g(d) = (1-d)^2\) before assembling, and read the energy from the returned strain.

The crack side — from FieldQuad4

The phase-field equation \[ \big(G_c\, l\, \nabla^2 - (G_c/l + 2H)\big)\, d = -2H \] is a screened-Poisson problem, and physicsbase's field element hands back the three canonical operators it's built from — the same ones that assemble a heat-conduction problem:

from femengine.field.elements import FieldQuad4

Kc = FieldQuad4.conductivity(sq, 1.0, {})   # (4,4)  ∫∇N·∇N   — the diffusion term
Mc = FieldQuad4.capacitance(sq, 1.0, {})    # (4,4)  ∫N·N     — the reaction term
Fc = FieldQuad4.source(sq, 1.0, {})         # (4,)   ∫N       — the forcing term

# assemble the crack system element by element:
#   A = Gc·l · Kc   +   (Gc/l + 2H) · Mc            (H = driving energy per element)
#   b = 2H · Fc
# solve A d = b, then clamp d to [d_prev, 1]  (cracks never heal)

That's the whole "field" side. conductivity gives the crack its surface-energy smoothing, capacitance supplies the reaction coefficient that couples in the driving energy H, and source is the forcing. Three matrix templates from the identical engine that solves temperature fields, reused verbatim to grow a crack.

The loop that couples them

for step in load_steps:                # ratchet the top displacement up
    for _ in range(2):                   # staggered iterations
        u = solve_elasticity(KE, g(d))   # ← Quad4, degraded by the crack
        H = np.maximum(H, psi(u))         # energy from the solve, kept monotonic
        d = solve_phase(Kc, Mc, Fc, H)   # ← FieldQuad4, driven by the energy

What came out

We pulled the top edge up in 40 steps on a 56×56 mesh. For most of them, nothing dramatic: the specimen stretches, the reaction climbs almost linearly, the pre-notch sits quietly. Then, right around a top displacement of 0.008, the notch tip reaches its breaking energy — and the crack runs.

That path was not prescribed. We seeded a notch on the left and pulled from the top; the horizontal, mode-I crack is what the energy competition produced. And the force tells the other half of the story — the signature of brittle fracture, a near-linear climb to a peak followed by a sudden collapse as the crack severs the specimen:

Why this is the hard case for an engine

Topology optimization was one physics in a loop. This is two physics in a loop, and it's a genuine stress test of an engine's composability: the same mesh has to serve a vector structural problem and a scalar field problem; element stiffness has to be scalable per-element; and the field solver has to accept an arbitrary spatially-varying reaction term. physicsbase's Quad4 and FieldQuad4 — validated separately against closed-form solutions — snap together here with nothing between them but two arrays of numbers passing back and forth. That's the payoff of treating simulation as an API: the coupling logic lives in fifty lines of glue, and each half of it stands on a solve the engine has already proven correct.

Honest footnotes

We implemented the standard AT2 phase-field model with a staggered solve — the formulation is common to the references below and to any modern fracture-mechanics course; the code is our own and reuses none of the published implementations. We used the simplest isotropic energy (no tension/compression split), a modest mesh, and two staggered iterations per step, which is enough to reproduce the benchmark's behaviour cleanly but not to chase the last few percent on the peak load or handle crack branching and compression-induced closure. Those refinements — spectral energy splits, finer meshes, monolithic or interior-point solvers — are exactly what the cited papers contribute. Our point was narrower and, we think, more interesting: a general-purpose FE engine, asked to couple two of its own solvers through nothing but fields, grows the textbook crack.

References

1. B. Bourdin, G. A. Francfort & J.-J. Marigo, "Numerical experiments in revisited brittle fracture," Journal of the Mechanics and Physics of Solids 48 (2000), 797–826.
2. C. Miehe, M. Hofacker & F. Welschinger, "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits," Comput. Methods Appl. Mech. Engrg. 199 (2010), 2765–2778.
3. Y. Bai et al., "An AsFem implementation for quasi-static brittle fracture phase field model" (2023). arXiv:2302.01510
4. Interior-point methods for the phase-field approach to brittle and ductile fracture (2020). arXiv:2011.10125