The 3D version of this sim reads off a closed-form critical crack length from a single algebraic formula and then animates a fixed propagation speed. This 2D version is a genuinely different computation: it numerically solves the Euler-Bernoulli beam-on-elastic-foundation equation for the slab bending over the collapsed weak layer, and derives both the critical length and the propagation dynamics from the actual elastic energy stored in that solved deflection field — the same finite-difference technique (Winkler-foundation beam models) used in real snow-mechanics literature to study weak-layer crack-tip stress concentration.
Governing PDE (beam on Winkler foundation):
EI·w''''(x) + k(x)·w(x) = q EI = E'D³/12
k(x) = 0 inside the crack, k₀ elsewhere
q = ρgD·cosψ (load pressing the slab onto the weak layer)
Solved by direct banded finite-difference elimination for w(x),
then the elastic + foundation energy Π(a) = U_bend + U_found − W_ext
gives the energy release rate: G(a) = −dΠ/da
Crack grows while G(a) > w_f (fracture energy); the critical length
a_c is found where G(a_c) = w_f — computed numerically every time
a slider changes, not read off a formula.
Independent numeric check (see verification script): for D=0.5 m, ψ=36°, w_f=2 J/m² this beam solve gives a_c ≈ 0.4–0.5 m, matching the field-measured Propagation Saw Test range of roughly 0.2–0.6 m. A plain dimensional check of the simple closed-form scaling law a_c = √(2E′w_f/(πτ_d²)) shows its units reduce to √m rather than m — it is missing a slab-thickness factor to be dimensionally a length. This solver sidesteps that entirely because thickness enters correctly through the bending rigidity EI = E′D³/12.
- G(a) < w_f — the beam solve shows the collapse stays local: an "whumpf" with no self-sustaining crack.
- G(a) ≥ w_f — because G(a) keeps increasing as the crack grows (verified numerically: more collapsed length always sags more, storing more bending energy), once critical the crack cannot arrest on its own — it accelerates toward the elastic-wave speed limit 0.35·√(E′/ρ) until the whole slab is isolated and slides at a = g(sinψ − μcosψ).
- ρ≈250 kg/m³ and E′≈2 MPa (slab) are fixed, matching the 3D sim's assumptions; the weak layer's own foundation stiffness (E≈0.25 MPa over 2 cm) is a separate fixed assumption unique to this beam model.