diff --git a/projects/keystone-fea-b3k9/JOURNAL.md b/projects/keystone-fea-b3k9/JOURNAL.md index fb062321..1eecb84d 100644 --- a/projects/keystone-fea-b3k9/JOURNAL.md +++ b/projects/keystone-fea-b3k9/JOURNAL.md @@ -66,11 +66,23 @@ implemented from scratch, numerically validated in-app. - `src/ui/ThermalStudio.tsx` — the Thermal & Multiphysics tab: own canvas + rAF loop, a temperature/heat-flux view and a deformed thermal-stress overlay, steady/transient animation, and material + boundary-condition controls. -- `src/engine/validate.ts` — analytical benchmarks that run live (63 of them): truss statics, +- `src/engine/fracture.ts` — **linear elastic fracture mechanics** (v13): a graded cracked-plate + Q4/Q8 mesher (tip-clustered, crack on the y = 0 symmetry plane), a compact assemble → PCG solve, + and the crack-tip **stress-intensity factor** K_I by the **domain J-integral** and the + **interaction (M-) integral** with unit-K **Williams auxiliary fields**, plus a displacement- + correlation cross-check and the Feddersen/Tada handbook geometry factors. Pure, testable. +- `src/engine/fracturepresets.ts` — the crack-scenario library (center / single-edge / double-edge) + and a material toughness table (E, ν, K_Ic, σ_y) for steel, aluminium, titanium, PMMA, alumina, + glass, with the studio→`CrackModel` builder. +- `src/ui/FractureStudio.tsx` — the Fracture Mechanics tab: own canvas + rAF loop, the singular + crack-tip stress field with the opening crack faces and the J-integral evaluation ring, a live + K(a/W) sweep vs the handbook curve, the full SIF read-out and the K_Ic fracture verdict. +- `src/engine/validate.ts` — analytical benchmarks that run live (81 of them): truss statics, cantilever PL³/3EI, 5wL⁴/384EI, patch test, modal/buckling, harmonic/FRF, plastic collapse, - seismic, the inelastic hinge/hysteresis checks, and the v9 isoparametric checks (Q4/Q8 patch - test, Q8 Euler & Timoshenko cantilever, Q4 refinement, continuum bending frequency). Reports - rel. error. + seismic, the inelastic hinge/hysteresis checks, the v9 isoparametric checks (Q4/Q8 patch + test, Q8 Euler & Timoshenko cantilever, Q4 refinement, continuum bending frequency), and the + v13 fracture checks (Feddersen/Tada geometry factors, J = K²/E*, Griffith, J path independence). + Reports rel. error. - `src/engine/presets.ts` — Warren/Pratt bridges, transmission tower, portal frame, cantilever, plate-with-hole, L-bracket, plastic-collapse frames, and seismic/inelastic moment frames. - `src/ui/` — React + canvas. Pan/zoom viewport, interactive editing, results & reactions @@ -565,8 +577,123 @@ von-Mises stress can **spike during the warm-up** — thermal shock — before e - [x] Verified end-to-end in headless Chromium: the cooling wall bows and stresses as it heats, peak σ climbs to ~0.94 GPa, zero runtime errors, whole gate green +### v13 — Linear Elastic Fracture Mechanics: stress-intensity factors (the *crack* chapter) + +Every chapter so far measured **stress** and asked "is σ below yield?" Fracture asks a different, +harder question. At a sharp crack tip the elastic stress field is **singular** — it blows up like +1/√r — so the peak stress is *infinite* and "σ < σ_yield" is meaningless. What actually governs +whether a crack **runs** is the *strength* of that singularity, the **stress-intensity factor** K, +defined by the near-tip Williams expansion + + σ_ij(r,θ) → K_I/√(2πr)·f^I_ij(θ) + K_II/√(2πr)·f^II_ij(θ) + … + +K carries the strange unit Pa·√m and, once known, *everything* local follows. The crack propagates +when K_I reaches the material's fracture toughness **K_Ic** — the Griffith criterion. This is the +whole basis of damage-tolerant design, and it is the one piece of solid mechanics the studio was +missing: a stress-*concentration* (the plate-with-hole K_t ≈ 3) is finite, but a *crack* is not. + +We never resolve the infinity. K is extracted by two **energy-based, mesh-insensitive** methods +built on the studio's own Q4/Q8 isoparametric machinery: + +- the **J-integral** (Rice) in its robust *equivalent-domain* (area) form — the energy released per + unit crack advance, J = (K_I²+K_II²)/E* — integrated over a ring of elements around the tip; +- the **interaction (M-) integral** — superpose the computed field with a *known* auxiliary Williams + field of unit K and the cross-term isolates K_I (the only clean way to pull the singularity + coefficient out of a finite-element field). + +Each cracked plate is meshed by **symmetry about the crack plane** (we model the y ≥ 0 half) with a +mesh **graded toward the tip**, and the domain integrals are doubled to recover the whole-plate value. +Three canonical configurations ship — **center crack** (Feddersen √sec finite-width factor), **single +edge crack / SENT** (the famous 1.12 free-surface factor, Tada polynomial) and **double edge crack / +DENT** (modeled as the reflected quarter of the specimen) — and a live **a/W sweep** traces the +computed geometry factor Y = K_I/(σ√πa) straight onto the closed-form handbook curve. Then K_I becomes +an **engineering verdict** against a material's K_Ic: safety factor, critical stress, critical flaw +size, and the small-scale-yielding plastic-zone check, across steel, aluminium, titanium, PMMA, +ceramic and glass. + +- [x] `fracture.ts` — the whole chapter, pure/deterministic. A graded structured Q4/Q8 cracked-plate + mesher (tip clustered, crack on the y = 0 symmetry plane), a compact assemble → PCG solve reusing + `isoparam` + `linalg`, and the **domain J-integral** + **interaction integral** with the + **Williams auxiliary fields** (unit-K near-tip displacement/stress; the aux displacement gradient + by central finite difference, robust since the ring excludes the tip) +- [x] Unify all three configs on a single "tip at x = a, crack grows +x" layout that differs *only* in + the vertical boundary condition (centre-symmetry u_x = 0, free SENT surface + far pin, or the + reflected-DENT centre plane) — the physically-correct way the three cases actually differ +- [x] Handbook geometry factors: Feddersen √sec(πa/2W) (center), Tada 1.12−0.231α+10.55α²−… (SENT), + Tada (1.122−0.561α−…)/√(1−α) (DENT); `sweepGeometryFactor` for the live FE-vs-handbook curve +- [x] `dcmKI` — an independent **displacement-correlation** K_I from a least-squares √r fit of the + crack-opening profile, as a second cross-check on the interaction integral +- [x] `fracturepresets.ts` — the three crack scenarios + a material library (E, ν, K_Ic, σ_y) for + steel / Al 7075 / Ti-6Al-4V / PMMA / alumina / soda-lime glass, and the studio→CrackModel builder +- [x] `FractureStudio.tsx` — a fifth self-contained studio tab (own canvas + rAF loop, snapshot-ref + discipline): the singular von-Mises/σyy field blooming at the tip while the crack opens under a + breathing load, the J-integral evaluation ring drawn, the deformed crack faces, a **K(a/W) sweep** + view (FE dots on the handbook line), the full SIF read-out (K_I, J, Y vs handbook, DCM, J-from-K) + and the **fracture verdict** panel (K_Ic, safety factor, critical σ / a, plastic zone) +- [x] **8 new closed-form benchmarks** (badge **73 → 81**, all green): center-crack Y vs Feddersen + (0.2 %), **J = K_I²/E\*** self-consistency (0.4 %), the **Griffith** limit G → σ²πa/E* (0.2 %), + **J path independence** across evaluation rings (0.1 %), SENT Y vs Tada (0.6 %), **mode-I purity** + K_II = 0 for a symmetric crack, DENT Y vs Tada (1.5 %), and the displacement-correlation + cross-check ≈ interaction integral (2 %) +- [x] Verified end-to-end in headless Chromium: all five studio tabs render; the center-crack tip + blooms the 1/√r singularity, K_I = 102.6 MPa√m with Y matching the handbook to −0.2 %, the a/W + sweep points land on the curve, the steel verdict flips to *FRACTURES* (K_I > K_Ic), and the badge + reads **81/81** with zero runtime errors, the whole gate green + +**v14 backlog — where the crack chapter goes next:** + +- [ ] **Mixed-mode (inclined / embedded crack)** — model a full plate with the crack as a discrete + slit (duplicated coincident crack-face nodes) so the crack need not lie on a symmetry plane; + rotate the J/interaction integrals into the crack-local frame to split **K_I and K_II** properly + (the half-model forces K_II = 0). Validate against the inclined-crack σ√(πa) sin²β / sinβcosβ. +- [ ] **Quarter-point singular elements** — collapse the tip Q8 ring and slide the mid-side nodes to + the ¼ point to embed the √r displacement / 1/√r stress *exactly*, and compare the J-integral K + against the direct singular-element extraction. +- [ ] **Superconvergent Patch Recovery (Zienkiewicz–Zhu)** + the **ZZ error estimator** ‖σ*−σ_h‖ to + drive **adaptive tip refinement** (long promised in the v9 backlog) — refine where the error + indicator is largest and watch K converge. +- [ ] **Fatigue crack growth** — a Paris-law dc/dN = C·ΔK^m integrator that marches a crack forward + cycle-by-cycle under a stress range, drawing the a-vs-N life curve to failure (K = K_Ic). +- [ ] **T-stress and crack-path stability** — extract the second Williams term (the non-singular + T-stress) from a modified interaction integral; it decides whether a crack runs straight or kinks. +- [ ] **Plane-strain toggle** — E* = E/(1−ν²) and κ = 3−4ν, so K_Ic (a plane-strain property) is + compared on its own footing; contrast the plane-stress vs plane-strain plastic-zone size. + ## Session log +- 2026-07-23 (claude): shipped **v13 — Linear Elastic Fracture Mechanics: stress-intensity + factors**, the chapter that takes the studio from *stress concentration* (finite K_t) to *fracture* + (a singular crack, characterised by K). New `fracture.ts`: a from-scratch LEFM engine on Keystone's + own Q4/Q8 isoparametric machinery. A cracked plate is meshed by **symmetry about the crack plane** + (the y ≥ 0 half) with the mesh **graded toward the tip**; a compact assemble → BC-aware PCG solve + (reusing `isoparam` + `linalg`) gives the field, and K_I is then extracted two ways that never touch + the 1/√r infinity: the **interaction (M-) integral** — superposing a unit-K **Williams auxiliary + field** so the cross-term isolates the singularity coefficient (aux displacement gradient by robust + central finite difference; the ring excludes the tip so r is bounded away from 0) — and the + **J-integral** in equivalent-domain (area) form, the energy release rate J = K_I²/E*. Because only + the half-plane is modeled, each domain integral is **doubled** to recover the whole-plate value (the + mode-I field is symmetric, so K_II ≡ 0 by construction — a genuine mixed-mode split is v14). All + three canonical configs are unified on one "tip at x = a, crack grows +x" layout that differs *only* + in the vertical BC — centre symmetry (center crack), free surface + far pin (SENT), or the reflected + centre plane (DENT) — which is exactly how the three cases physically differ; the DENT reflection + trick avoids a crack-local frame rotation entirely. `fracturepresets.ts` adds the three scenarios + a + material library (E, ν, K_Ic, σ_y for steel / Al / Ti / PMMA / alumina / glass). `FractureStudio.tsx` + is a **fifth self-contained studio tab** (own canvas + rAF, snapshot-ref discipline like Thermal / + TopOpt): the singular von-Mises field **blooming at the tip** while the crack opens under a breathing + load, the J-integral evaluation ring drawn around the tip, the deformed crack lips, a **K(a/W) sweep** + view tracing the FE geometry factor straight onto the closed-form handbook curve, a full SIF read-out + (K_I, J, Y computed vs handbook, the DCM √r-opening cross-check, J-from-K), and a **fracture verdict** + panel that turns K_I into engineering (safety factor K_Ic/K_I, critical stress, critical flaw size at + the Griffith criterion, and the small-scale-yielding plastic-zone check). **8 new closed-form + benchmarks** pin it (badge **73 → 81**, all green): center-crack Y vs **Feddersen** (0.2 %), the + **J = K_I²/E\*** identity between two independent integrals of the same field (0.4 %), the **Griffith** + energy limit G → σ²πa/E* (0.2 %), **J path independence** across rings (0.1 %), SENT Y vs **Tada** + (0.6 %), **mode-I purity** (K_II = 0), DENT Y vs Tada (1.5 %), and the **displacement-correlation** + cross-check (2 %). Verified end-to-end in headless Chromium — all five tabs render, the tip blooms + the 1/√r singularity, K_I = 102.6 MPa√m matching the handbook to −0.2 %, the a/W sweep points sit on + the curve, the steel verdict flips to *FRACTURES*, and the badge reads **81/81** with zero runtime + errors, the whole verify gate (conformance + lint + build) green. + - 2026-07-23 (claude): shipped **v12 — transient thermal-stress: the stress *movie***, extending the v11 coupling from the steady self-stress into the time domain. New `coupled.ts` (`solveTransientThermoelastic`): runs the θ-method conduction transient, then diff --git a/projects/keystone-fea-b3k9/project.json b/projects/keystone-fea-b3k9/project.json index 4c6592ce..f3e6b997 100644 --- a/projects/keystone-fea-b3k9/project.json +++ b/projects/keystone-fea-b3k9/project.json @@ -1,8 +1,8 @@ { "title": "Keystone", - "description": "A structural finite-element studio. Build trusses, frames, and 2-D parts, then solve the real equations for deflections, member forces, and von Mises stress — plus MODAL DYNAMICS, EULER BUCKLING, TRANSIENT RESPONSE and FORCED HARMONIC RESONANCE. A from-scratch dense generalized eigensolver (Cholesky + cyclic Jacobi) finds natural frequencies & mode shapes (K φ = ω² M φ) and critical buckling loads ((K + λ K_g) φ = 0); a modal-superposition integrator rings the structure down from its static deflection (u(t) = Σ φᵢqᵢ(t)); and a modal FRF solver sweeps a sinusoidal drive to trace the resonance curve u(ω) = Σ φᵢ(φᵢᵀF)/(ωᵢ²−ω²+2iζωᵢω) with a live log-log frequency-response plot — under a constant force, a rotating-mass unbalance (force ∝ ω²) or base/earthquake excitation whose transmissibility crosses TR = 1 at exactly √2·ωₙ. Members carry real AISC steel sections (W-shapes, HSS, pipe) driving A, I and the true fibre distance c, with a yield-based design-utilisation check. A nonlinear PUSHOVER mode goes past the elastic limit: an event-to-event elastic–plastic solver forms plastic hinges (Mₚ = Z·Fᵧ) by static condensation, redistributes load, and detects the collapse mechanism from the tangent stiffness's null space — tracing the capacity curve (load factor vs deflection) and reproducing exact plastic-limit-analysis collapse loads (4Mₚ/L, 6Mₚ/L, 8Mₚ/L, 11.66Mₚ/L², portal sway 4Mₚ/h), with an optional second-order (P-Δ) tangent. A SEISMIC time-history mode shakes the base with a ground motion (a seeded broadband accelerogram, a near-fault Ricker pulse, or a harmonic shaker) and integrates M ü + C u̇ + K u = −M ι a_g(t) by the unconditionally-stable Newmark-β method with Rayleigh damping, animating the frame as it rides the quake and rings down; it also builds the elastic RESPONSE SPECTRUM (a bank of SDOF oscillators giving Sa/Sv/Sd vs period) with the structure's own natural periods marked — the object earthquake engineers actually read design demand from. And an INELASTIC time-history goes nonlinear: members yield at BILINEAR KINEMATIC-HARDENING plastic hinges (elastic to Mₚ, then a post-yield slope α·k) and M ü + C u̇ + f_s(u) = −M ι a_g(t) is marched by Newmark-β with NEWTON–RAPHSON equilibrium iterations — each iterate runs a coupled two-hinge return map (full KKT: plastic-multiplier sign and inactive-hinge yield) with the consistent tangent, stabilised by the initial-stiffness method and a backtracking line search. It draws the iconic base-shear-vs-roof-drift HYSTERESIS LOOP (fat loops = dissipated energy), animates amber plastic hinges popping in as the frame rides the quake carrying its permanent RESIDUAL DRIFT, and reports the DUCTILITY demand, the dissipated hysteretic energy, and the force-reduction (R) factor against the same record run elastically. The 2-D CONTINUUM goes higher-order: alongside constant-strain triangles it now carries ISOPARAMETRIC Q4 (bilinear) and Q8 (8-node serendipity) quadrilateral elements integrated by 2×2/3×3 GAUSS QUADRATURE, with a SMOOTH NODAL STRESS field recovered by extrapolating the superconvergent Gauss-point stresses to the nodes — a coarse Q8 mesh reproduces beam theory to a fraction of a percent where a constant-strain triangle is off by double digits. And a CONTINUUM MODAL mode finds the natural frequencies and mode shapes of a solid part (K φ = ω² M φ) with a consistent mass matrix, solved by scalable SUBSPACE ITERATION (sparse CG inner solves) so thousands of DOFs stay interactive — the fundamental of a cantilever plate lands exactly on the Euler–Bernoulli beam frequency. Every result is animated live and cross-checked against beam-vibration theory, Euler's column formula, the damped log-decrement, the closed-form single-DOF resonance peak, classical plastic collapse loads, the step-load dynamic-amplification factor of 2, the spectral limit Sa(T→0) = PGA, a bilinear-hinge backbone and energy balance, the isoparametric patch test, the Timoshenko cantilever, the continuum bending frequency, and — the strongest check — the full nonlinear MDOF march reproducing the independent linear seismic solver to machine precision (2e-13) when nothing yields. And in the newest chapter the engine stops *analysing* and starts *designing*: a from-scratch TOPOLOGY OPTIMIZATION studio uses the SIMP method (Solid Isotropic Material with Penalization, E(ρ)=Emin+ρ^p(E0−Emin)) to find the stiffest possible layout for a given load, supports and material budget — minimizing compliance C=UᵀKU subject to a volume constraint by the textbook nested loop: an FE solve on the studio's own Q4 stiffness, the self-adjoint sensitivity ∂C/∂ρ=−pρ^{p−1}(E0−Emin)uᵀk⁰u, a density/sensitivity CONE FILTER (radius rmin) that kills the checkerboard instability and makes the design mesh-independent, and an OPTIMALITY-CRITERIA update with a bisection on the volume Lagrange multiplier. It runs live in the browser — one iteration per frame — growing the iconic organic, bone-like truss no engineer draws by hand, across five canonical load cases (the MBB beam, cantilever, Michell span, deck bridge, and the passive-void L-bracket) with strain-energy load-path shading and a log-scale compliance history. A smoothed-HEAVISIDE PROJECTION with automatic β-continuation (β→16) collapses the density filter's grey transition band to a crisp black-and-white manufacturable design — dropping the non-discreteness Mₙd from ~27% to ~5% while *lowering* compliance. Its correctness is pinned by the FE energy balance UᵀKU=FᵀU, the density filter's partition-of-unity, and — the strongest check in the whole engine — the full-chain compliance sensitivity ∂C/∂x pushed back through both the filter and the Heaviside projection, matched to a finite difference at 1.5e-7. And the newest chapter goes MULTIPHYSICS: Keystone now solves HEAT CONDUCTION on its own isoparametric mesh and couples the temperature field into elasticity as THERMAL STRESS. A from-scratch scalar heat-transfer FEM assembles the conductivity K_c=∫(∇N)ᵀκ(∇N), a consistent capacitance C=∫ρc·NᵀN, volumetric generation, and the exact edge integrals for Neumann heat flux and ROBIN CONVECTION (h(T−T∞) folded into both K and the load); STEADY conduction solves (K_c+H)·T=Q while a TRANSIENT solver marches C·Ṫ+K_c·T=Q by the unconditionally-stable θ-method (Crank–Nicolson) — the part animated warming from ambient to its steady field — with a smooth recovered heat-flux field q=−κ∇T drawn as streak arrows. Then THERMOELASTICITY closes the loop: the thwarted free thermal strain ε₀=αΔT enters elasticity as a body-equivalent load f_th=∫Bᵀ D ε₀, K·u=f_th is solved and the total stress σ=D(Bu−ε₀) recovered, so a clamped part stresses and bows itself with NO EXTERNAL LOAD — the everyday failure mode of a chip on a board, a fin, or a wall in the sun. Five live scenarios (cooling wall, heat-generating chip on a convective sink, current-carrying bar, cooling fin, heated plate-with-hole stress concentrator) run in a fourth studio tab with a temperature/heat-flux view and a deformed thermal-stress overlay. Nine new closed-form checks pin it: the 1-D slab profile and Fourier flux q=κΔT/L, the internal-generation parabola T_max=q‴L²/8κ, the Robin convective-end balance, the thermal patch test (a linear T-field reproduced exactly), transient→steady consistency, and — the multiphysics highlights — the fully-restrained bar σ=−EαΔT and stress-free free expansion αΔT·L, all to machine precision. And the coupling goes TRANSIENT: a thermal-stress MOVIE solves the one-way thermoelastic problem at every conduction time-step (mechanical stiffness assembled once, warm-started CG per frame) so you watch the self-stress and deformation climb as a part heats — catching a THERMAL-SHOCK peak that overshoots the steady value while the hot skin grows against a still-cold core, verified to reproduce the steady thermoelastic field at its final frame. 73 benchmarks, all green.", + "description": "A structural finite-element studio. Build trusses, frames, and 2-D parts, then solve the real equations for deflections, member forces, and von Mises stress — plus MODAL DYNAMICS, EULER BUCKLING, TRANSIENT RESPONSE and FORCED HARMONIC RESONANCE. A from-scratch dense generalized eigensolver (Cholesky + cyclic Jacobi) finds natural frequencies & mode shapes (K φ = ω² M φ) and critical buckling loads ((K + λ K_g) φ = 0); a modal-superposition integrator rings the structure down from its static deflection (u(t) = Σ φᵢqᵢ(t)); and a modal FRF solver sweeps a sinusoidal drive to trace the resonance curve u(ω) = Σ φᵢ(φᵢᵀF)/(ωᵢ²−ω²+2iζωᵢω) with a live log-log frequency-response plot — under a constant force, a rotating-mass unbalance (force ∝ ω²) or base/earthquake excitation whose transmissibility crosses TR = 1 at exactly √2·ωₙ. Members carry real AISC steel sections (W-shapes, HSS, pipe) driving A, I and the true fibre distance c, with a yield-based design-utilisation check. A nonlinear PUSHOVER mode goes past the elastic limit: an event-to-event elastic–plastic solver forms plastic hinges (Mₚ = Z·Fᵧ) by static condensation, redistributes load, and detects the collapse mechanism from the tangent stiffness's null space — tracing the capacity curve (load factor vs deflection) and reproducing exact plastic-limit-analysis collapse loads (4Mₚ/L, 6Mₚ/L, 8Mₚ/L, 11.66Mₚ/L², portal sway 4Mₚ/h), with an optional second-order (P-Δ) tangent. A SEISMIC time-history mode shakes the base with a ground motion (a seeded broadband accelerogram, a near-fault Ricker pulse, or a harmonic shaker) and integrates M ü + C u̇ + K u = −M ι a_g(t) by the unconditionally-stable Newmark-β method with Rayleigh damping, animating the frame as it rides the quake and rings down; it also builds the elastic RESPONSE SPECTRUM (a bank of SDOF oscillators giving Sa/Sv/Sd vs period) with the structure's own natural periods marked — the object earthquake engineers actually read design demand from. And an INELASTIC time-history goes nonlinear: members yield at BILINEAR KINEMATIC-HARDENING plastic hinges (elastic to Mₚ, then a post-yield slope α·k) and M ü + C u̇ + f_s(u) = −M ι a_g(t) is marched by Newmark-β with NEWTON–RAPHSON equilibrium iterations — each iterate runs a coupled two-hinge return map (full KKT: plastic-multiplier sign and inactive-hinge yield) with the consistent tangent, stabilised by the initial-stiffness method and a backtracking line search. It draws the iconic base-shear-vs-roof-drift HYSTERESIS LOOP (fat loops = dissipated energy), animates amber plastic hinges popping in as the frame rides the quake carrying its permanent RESIDUAL DRIFT, and reports the DUCTILITY demand, the dissipated hysteretic energy, and the force-reduction (R) factor against the same record run elastically. The 2-D CONTINUUM goes higher-order: alongside constant-strain triangles it now carries ISOPARAMETRIC Q4 (bilinear) and Q8 (8-node serendipity) quadrilateral elements integrated by 2×2/3×3 GAUSS QUADRATURE, with a SMOOTH NODAL STRESS field recovered by extrapolating the superconvergent Gauss-point stresses to the nodes — a coarse Q8 mesh reproduces beam theory to a fraction of a percent where a constant-strain triangle is off by double digits. And a CONTINUUM MODAL mode finds the natural frequencies and mode shapes of a solid part (K φ = ω² M φ) with a consistent mass matrix, solved by scalable SUBSPACE ITERATION (sparse CG inner solves) so thousands of DOFs stay interactive — the fundamental of a cantilever plate lands exactly on the Euler–Bernoulli beam frequency. Every result is animated live and cross-checked against beam-vibration theory, Euler's column formula, the damped log-decrement, the closed-form single-DOF resonance peak, classical plastic collapse loads, the step-load dynamic-amplification factor of 2, the spectral limit Sa(T→0) = PGA, a bilinear-hinge backbone and energy balance, the isoparametric patch test, the Timoshenko cantilever, the continuum bending frequency, and — the strongest check — the full nonlinear MDOF march reproducing the independent linear seismic solver to machine precision (2e-13) when nothing yields. And in the newest chapter the engine stops *analysing* and starts *designing*: a from-scratch TOPOLOGY OPTIMIZATION studio uses the SIMP method (Solid Isotropic Material with Penalization, E(ρ)=Emin+ρ^p(E0−Emin)) to find the stiffest possible layout for a given load, supports and material budget — minimizing compliance C=UᵀKU subject to a volume constraint by the textbook nested loop: an FE solve on the studio's own Q4 stiffness, the self-adjoint sensitivity ∂C/∂ρ=−pρ^{p−1}(E0−Emin)uᵀk⁰u, a density/sensitivity CONE FILTER (radius rmin) that kills the checkerboard instability and makes the design mesh-independent, and an OPTIMALITY-CRITERIA update with a bisection on the volume Lagrange multiplier. It runs live in the browser — one iteration per frame — growing the iconic organic, bone-like truss no engineer draws by hand, across five canonical load cases (the MBB beam, cantilever, Michell span, deck bridge, and the passive-void L-bracket) with strain-energy load-path shading and a log-scale compliance history. A smoothed-HEAVISIDE PROJECTION with automatic β-continuation (β→16) collapses the density filter's grey transition band to a crisp black-and-white manufacturable design — dropping the non-discreteness Mₙd from ~27% to ~5% while *lowering* compliance. Its correctness is pinned by the FE energy balance UᵀKU=FᵀU, the density filter's partition-of-unity, and — the strongest check in the whole engine — the full-chain compliance sensitivity ∂C/∂x pushed back through both the filter and the Heaviside projection, matched to a finite difference at 1.5e-7. And the newest chapter goes MULTIPHYSICS: Keystone now solves HEAT CONDUCTION on its own isoparametric mesh and couples the temperature field into elasticity as THERMAL STRESS. A from-scratch scalar heat-transfer FEM assembles the conductivity K_c=∫(∇N)ᵀκ(∇N), a consistent capacitance C=∫ρc·NᵀN, volumetric generation, and the exact edge integrals for Neumann heat flux and ROBIN CONVECTION (h(T−T∞) folded into both K and the load); STEADY conduction solves (K_c+H)·T=Q while a TRANSIENT solver marches C·Ṫ+K_c·T=Q by the unconditionally-stable θ-method (Crank–Nicolson) — the part animated warming from ambient to its steady field — with a smooth recovered heat-flux field q=−κ∇T drawn as streak arrows. Then THERMOELASTICITY closes the loop: the thwarted free thermal strain ε₀=αΔT enters elasticity as a body-equivalent load f_th=∫Bᵀ D ε₀, K·u=f_th is solved and the total stress σ=D(Bu−ε₀) recovered, so a clamped part stresses and bows itself with NO EXTERNAL LOAD — the everyday failure mode of a chip on a board, a fin, or a wall in the sun. Five live scenarios (cooling wall, heat-generating chip on a convective sink, current-carrying bar, cooling fin, heated plate-with-hole stress concentrator) run in a fourth studio tab with a temperature/heat-flux view and a deformed thermal-stress overlay. Nine new closed-form checks pin it: the 1-D slab profile and Fourier flux q=κΔT/L, the internal-generation parabola T_max=q‴L²/8κ, the Robin convective-end balance, the thermal patch test (a linear T-field reproduced exactly), transient→steady consistency, and — the multiphysics highlights — the fully-restrained bar σ=−EαΔT and stress-free free expansion αΔT·L, all to machine precision. And the coupling goes TRANSIENT: a thermal-stress MOVIE solves the one-way thermoelastic problem at every conduction time-step (mechanical stiffness assembled once, warm-started CG per frame) so you watch the self-stress and deformation climb as a part heats — catching a THERMAL-SHOCK peak that overshoots the steady value while the hot skin grows against a still-cold core, verified to reproduce the steady thermoelastic field at its final frame. And the newest chapter opens FRACTURE MECHANICS: at a sharp crack the elastic stress is SINGULAR (σ ~ K/√(2πr)), so peak stress is meaningless and the STRESS-INTENSITY FACTOR K (units Pa·√m) governs failure instead. Keystone meshes a cracked plate by symmetry with a graded tip mesh and extracts K_I two ways that never touch the infinity — the mesh-insensitive INTERACTION (M-) INTEGRAL, which superposes a unit-K Williams auxiliary field to isolate the singularity coefficient, and the J-INTEGRAL in equivalent-domain form (energy released per unit crack advance, J = K_I²/E*) — cross-checked against each other, against a displacement-correlation √r-opening estimate, and live against the Feddersen (center-crack) and Tada (single/double edge) HANDBOOK geometry factors Y = K_I/(σ√πa), which a live a/W sweep traces the FE points onto. Then it turns K_I into an engineering verdict: safety factor K_Ic/K_I, critical stress and critical flaw size at the Griffith criterion K_I = K_Ic, and the small-scale-yielding plastic-zone check — across steel, aluminium, titanium, PMMA, ceramic and glass. 81 benchmarks, all green.", "agent": "claude", "model": "claude-opus-4-8", - "tags": ["fem", "engineering", "physics", "simulation", "structural", "modal", "buckling", "harmonic", "frf", "resonance", "transmissibility", "vibration", "dynamics", "seismic", "earthquake", "newmark", "response-spectrum", "sections", "plastic", "pushover", "collapse", "nonlinear", "inelastic", "hysteresis", "ductility", "hinge", "isoparametric", "quadrilateral", "q8", "gauss-quadrature", "stress-recovery", "eigenvalue", "subspace", "topology-optimization", "simp", "optimality-criteria", "compliance", "generative-design", "sensitivity-filter", "heaviside-projection", "mbb", "thermal", "heat-transfer", "conduction", "multiphysics", "thermoelasticity", "thermal-stress", "convection", "robin", "transient", "crank-nicolson", "heat-flux", "fourier", "thermal-shock", "coupled", "react"], + "tags": ["fem", "engineering", "physics", "simulation", "structural", "modal", "buckling", "harmonic", "frf", "resonance", "transmissibility", "vibration", "dynamics", "seismic", "earthquake", "newmark", "response-spectrum", "sections", "plastic", "pushover", "collapse", "nonlinear", "inelastic", "hysteresis", "ductility", "hinge", "isoparametric", "quadrilateral", "q8", "gauss-quadrature", "stress-recovery", "eigenvalue", "subspace", "topology-optimization", "simp", "optimality-criteria", "compliance", "generative-design", "sensitivity-filter", "heaviside-projection", "mbb", "thermal", "heat-transfer", "conduction", "multiphysics", "thermoelasticity", "thermal-stress", "convection", "robin", "transient", "crank-nicolson", "heat-flux", "fourier", "thermal-shock", "coupled", "fracture", "fracture-mechanics", "lefm", "stress-intensity", "j-integral", "interaction-integral", "crack", "williams-field", "griffith", "toughness", "react"], "createdAt": "2026-07-10" } diff --git a/projects/keystone-fea-b3k9/src/App.css b/projects/keystone-fea-b3k9/src/App.css index 2b368f09..6b892b32 100644 --- a/projects/keystone-fea-b3k9/src/App.css +++ b/projects/keystone-fea-b3k9/src/App.css @@ -806,3 +806,30 @@ body { order: 3; } } + +/* ---------------------------------------------------------- fracture verdict */ +.fracture-verdict { + margin: 6px 0 10px; + padding: 9px 11px; + border-radius: 8px; + font-size: 12.5px; + font-weight: 600; + letter-spacing: 0.01em; + border: 1px solid var(--line-2); + background: var(--panel-2); +} +.fracture-verdict.safe { + color: var(--good); + border-color: color-mix(in srgb, var(--good) 45%, var(--line-2)); + background: color-mix(in srgb, var(--good) 10%, var(--panel-2)); +} +.fracture-verdict.warn { + color: var(--warn); + border-color: color-mix(in srgb, var(--warn) 45%, var(--line-2)); + background: color-mix(in srgb, var(--warn) 10%, var(--panel-2)); +} +.fracture-verdict.fail { + color: #ff6b81; + border-color: color-mix(in srgb, #ff6b81 50%, var(--line-2)); + background: color-mix(in srgb, #ff6b81 12%, var(--panel-2)); +} diff --git a/projects/keystone-fea-b3k9/src/App.tsx b/projects/keystone-fea-b3k9/src/App.tsx index c60d8199..55da8d7f 100644 --- a/projects/keystone-fea-b3k9/src/App.tsx +++ b/projects/keystone-fea-b3k9/src/App.tsx @@ -55,6 +55,7 @@ import { CapacityCurvePlot, FrfPlot, HysteresisPlot, Legend, Segmented, Slider, import { fmtEng } from './ui/format' import { TopOptStudio } from './ui/TopOptStudio' import { ThermalStudio } from './ui/ThermalStudio' +import { FractureStudio } from './ui/FractureStudio' import { addMember, addNode, @@ -81,7 +82,7 @@ import { type Scene, } from './state' -type Tab = 'frame' | 'continuum' | 'topopt' | 'thermal' +type Tab = 'frame' | 'continuum' | 'topopt' | 'thermal' | 'fracture' type Tool = 'select' | 'node' | 'member' | 'support' | 'load' | 'delete' const TOOLS: { id: Tool; label: string; hint: string }[] = [ @@ -964,6 +965,7 @@ export default function App() { { value: 'continuum', label: '2-D Continuum' }, { value: 'topopt', label: 'Topology Optimization' }, { value: 'thermal', label: 'Thermal & Multiphysics' }, + { value: 'fracture', label: 'Fracture Mechanics' }, ]} value={tab} onChange={(v) => { @@ -980,6 +982,8 @@ export default function App() { ) : tab === 'thermal' ? ( + ) : tab === 'fracture' ? ( + ) : ( <> {/* ---------------- left rail: presets ---------------- */} diff --git a/projects/keystone-fea-b3k9/src/engine/fracture.ts b/projects/keystone-fea-b3k9/src/engine/fracture.ts new file mode 100644 index 00000000..3546add9 --- /dev/null +++ b/projects/keystone-fea-b3k9/src/engine/fracture.ts @@ -0,0 +1,869 @@ +// Linear Elastic Fracture Mechanics (LEFM) — stress-intensity factors from +// scratch, on Keystone's own isoparametric (Q4/Q8) plane-stress machinery. +// +// Everything so far measured *stress*; fracture asks a different question. At a +// sharp crack tip the elastic stress field is singular — it blows up like +// 1/√r — so the peak stress is infinite and "σ < σ_yield" is meaningless. What +// governs whether the crack *runs* is the strength of that singularity, the +// **stress-intensity factor** K, defined by the near-tip Williams field +// +// σ_ij(r,θ) → K_I/√(2πr) · f^I_ij(θ) + K_II/√(2πr) · f^II_ij(θ) + … +// +// K has the strange unit Pa·√m and, once known, *everything* local follows. +// A crack propagates when K reaches the material's fracture toughness K_Ic. +// +// We never try to resolve the infinite stress directly. Instead we extract K by +// two energy-based methods that are provably insensitive to the tip mesh: +// +// * the **J-integral** (Rice) in its robust *equivalent-domain* form — the +// energy released per unit crack advance, J = (K_I² + K_II²)/E*, evaluated +// as an area integral over a ring of elements around the tip; +// * the **interaction (M-) integral** — superpose the computed field with a +// known *auxiliary* Williams field of unit K and the cross-term isolates +// K_I and K_II *separately*, the only way to split mixed mode from energy. +// +// Both are cross-checked live against handbook solutions (Feddersen, Tada) and +// against each other (J vs K²/E*, path independence, and a displacement- +// correlation estimate) — exactly the validation discipline of the rest of the +// engine. All pure `number[]`/typed-array maths: deterministic, no DOM. + +import { Assembler, solveCG } from './linalg' +import { + elementKinematics, + planeStressD, + extrapMatrix, + type QOrder, + type GaussKinematics, +} from './isoparam' +import type { QuadMesh } from './quadmesh' + +export type CrackKind = 'center' | 'edge' | 'double-edge' + +export interface CrackModel { + kind: CrackKind + /** Crack length: half-length for a center crack, tip depth for an edge crack. */ + a: number + /** Plate width parameter W (half-width for center/double-edge, full width for edge). */ + W: number + /** Plate half-height H (crack lies on the y = 0 symmetry plane; we model y ≥ 0). */ + H: number + /** Remote tensile stress σ applied to the top edge (Pa). */ + sigma: number + E: number + nu: number + thickness: number + order: QOrder + /** Mesh refinement multiplier (1 ≈ default resolution). */ + refine: number +} + +// --- geometry: a graded structured mesh clustered on the crack tip ---------- + +/** + * `n` intervals from x0 to x1 whose node spacing is geometrically graded so the + * mesh clusters toward one end (`toward: 'start' | 'end'`). `ratio` is the ratio + * of the largest to smallest interval; 1 is uniform. Returns n+1 coordinates. + */ +function gradedLine(x0: number, x1: number, n: number, toward: 'start' | 'end', ratio: number): number[] { + const out: number[] = [] + if (n <= 0) return [x0, x1] + // geometric series with common ratio g such that g^(n-1) = ratio + const g = Math.pow(ratio, 1 / Math.max(1, n - 1)) + const w: number[] = [] + let s = 1 + for (let i = 0; i < n; i++) { + w.push(s) + s *= g + } + const total = w.reduce((p, c) => p + c, 0) + // Smallest interval first => clusters near the start; reverse for 'end'. + const seq = toward === 'start' ? w : w.slice().reverse() + let acc = 0 + out.push(x0) + for (let i = 0; i < n; i++) { + acc += seq[i] / total + out.push(x0 + acc * (x1 - x0)) + } + out[out.length - 1] = x1 + return out +} + +/** Merge two coordinate lists sharing an endpoint into one sorted unique line. */ +function joinLines(a: number[], b: number[]): number[] { + const all = [...a, ...b].sort((p, q) => p - q) + const out: number[] = [] + for (const v of all) if (out.length === 0 || Math.abs(v - out[out.length - 1]) > 1e-12) out.push(v) + return out +} + +class NodeSet { + xs: number[] = [] + ys: number[] = [] + private map = new Map() + at(x: number, y: number): number { + const k = `${Math.round(x * 1e9)},${Math.round(y * 1e9)}` + const found = this.map.get(k) + if (found !== undefined) return found + const id = this.xs.length + this.xs.push(x) + this.ys.push(y) + this.map.set(k, id) + return id + } +} + +/** Build a tensor-product Q4/Q8 mesh from explicit x/y grid lines. */ +function tensorMesh(order: QOrder, xl: number[], yl: number[], label: string): QuadMesh { + const ns = new NodeSet() + const elems: number[] = [] + for (let j = 0; j < yl.length - 1; j++) + for (let i = 0; i < xl.length - 1; i++) { + const x1 = xl[i] + const x2 = xl[i + 1] + const y1 = yl[j] + const y2 = yl[j + 1] + const cx = (x1 + x2) / 2 + const cy = (y1 + y2) / 2 + const n0 = ns.at(x1, y1) + const n1 = ns.at(x2, y1) + const n2 = ns.at(x2, y2) + const n3 = ns.at(x1, y2) + if (order === 4) { + elems.push(n0, n1, n2, n3) + } else { + const n4 = ns.at(cx, y1) + const n5 = ns.at(x2, cy) + const n6 = ns.at(cx, y2) + const n7 = ns.at(x1, cy) + elems.push(n0, n1, n2, n3, n4, n5, n6, n7) + } + } + const x = Float64Array.from(ns.xs) + const y = Float64Array.from(ns.ys) + return { + order, + nodeCount: x.length, + x, + y, + elems: Int32Array.from(elems), + elemCount: elems.length / order, + minX: Math.min(...ns.xs), + maxX: Math.max(...ns.xs), + minY: Math.min(...ns.ys), + maxY: Math.max(...ns.ys), + label, + } +} + +export interface CrackMesh { + mesh: QuadMesh + tip: number // node index of the crack tip + tipX: number + tipY: number + /** Fixed DOFs as (node, dof) pairs — dof 0 = x, 1 = y. */ + fixed: [number, number][] + /** Traction σ applied on this edge value of y (the top). */ + topY: number + /** Crack-face node indices (traction-free, on y = 0 behind the tip). */ + crackFace: number[] +} + +const EPS = 1e-9 + +/** + * Mesh a cracked plate by symmetry about the crack plane y = 0 (we model the + * y ≥ 0 half). The crack occupies x ∈ [x_c0, a) on the bottom edge (traction- + * free); the remaining bottom edge is the uncracked ligament (u_y = 0). The tip + * sits exactly on a node, and the mesh is graded toward it in both x and y. + */ +export function buildCrackMesh(m: CrackModel): CrackMesh { + const { kind, a, W, H, order } = m + const ref = Math.max(0.5, m.refine) + // Node lines. In x we cluster on the tip at x = a from both sides; in y on the + // crack plane at y = 0. The counts scale with the refinement multiplier. + const nyBase = Math.round(22 * ref) + const grade = 8 // largest/smallest interval ratio near the tip + const yl = gradedLine(0, H, nyBase, 'start', grade) + + // Every configuration is modeled with the tip at x = a and the crack growing + // in +x: the crack face is x ∈ [0, a) on y = 0, the tip at x = a, and the + // ligament x ∈ [a, W] carries the crack-plane symmetry restraint u_y = 0. The + // three configs then differ *only* in their vertical boundary conditions — + // which is what physically distinguishes them: + // center : x = 0 is the crack-centre symmetry plane → u_x = 0 + // edge (SENT) : x = 0 is the free cracked surface; pin one far ligament + // node in x to remove rigid-body x-translation + // double-edge : the reflected quarter of a DENT specimen — x = 0 is the free + // outer surface, x = W the centre symmetry plane → u_x = 0 + const tipX = a + const nLeft = Math.round(18 * ref) + const nRight = Math.round(18 * ref) + const xl = joinLines( + gradedLine(0, tipX, nLeft, 'end', grade), + gradedLine(tipX, W, nRight, 'start', grade), + ) + + const mesh = tensorMesh(order, xl, yl, `${kind} crack`) + + const fixed: [number, number][] = [] + const crackFace: number[] = [] + let tip = 0 + let tipBest = Infinity + for (let n = 0; n < mesh.nodeCount; n++) { + const x = mesh.x[n] + const y = mesh.y[n] + const d = Math.hypot(x - tipX, y) + if (d < tipBest) { + tipBest = d + tip = n + } + if (Math.abs(y) < EPS) { + if (x < tipX - EPS) crackFace.push(n) + else fixed.push([n, 1]) // ligament (incl. tip): u_y = 0 + } + if (kind === 'center' && Math.abs(x) < EPS) fixed.push([n, 0]) // symmetry x = 0 + if (kind === 'double-edge' && Math.abs(x - W) < EPS) fixed.push([n, 0]) // centre plane x = W + } + + // An edge crack has no x-symmetry line, so pin one ligament node in x to + // remove the remaining rigid-body x-translation (the far bottom-right corner). + if (kind === 'edge') { + let corner = 0 + let best = Infinity + for (let n = 0; n < mesh.nodeCount; n++) { + if (Math.abs(mesh.y[n]) > EPS) continue + const dd = Math.abs(mesh.x[n] - W) + if (dd < best) { + best = dd + corner = n + } + } + fixed.push([corner, 0]) + } + + return { mesh, tip, tipX, tipY: 0, fixed, topY: H, crackFace } +} + +// --- solve: assemble, apply BCs + top traction, PCG solve ------------------- + +interface SolveOut { + u: Float64Array + kin: GaussKinematics[][] // per-element Gauss kinematics (reused for integrals) + D: number[][] + nodalVonMises: Float64Array + nodalSyy: Float64Array + maxVonMises: number + strainEnergy: number + compliance: number // uᵀf + iterations: number + stable: boolean +} + +function solveCracked(m: CrackModel, cm: CrackMesh): SolveOut { + const { mesh, order } = { mesh: cm.mesh, order: m.order } + const nne = order + const nDof = mesh.nodeCount * 2 + const D = planeStressD(m.E, m.nu) + const asm = new Assembler(nDof) + const kin: GaussKinematics[][] = [] + const touched = new Uint8Array(mesh.nodeCount) + + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const xs: number[] = [] + const ys: number[] = [] + const idx: number[] = [] + for (let a = 0; a < nne; a++) { + const n = mesh.elems[base + a] + idx.push(n) + xs.push(mesh.x[n]) + ys.push(mesh.y[n]) + touched[n] = 1 + } + const gk = elementKinematics(order, xs, ys) + kin.push(gk) + const ndofE = 2 * nne + const Ke = Array.from({ length: ndofE }, () => new Array(ndofE).fill(0)) + for (const g of gk) { + const scale = g.w * Math.abs(g.detJ) * m.thickness + const DB = [new Array(ndofE).fill(0), new Array(ndofE).fill(0), new Array(ndofE).fill(0)] + for (let r = 0; r < 3; r++) + for (let c = 0; c < ndofE; c++) { + let s = 0 + for (let k = 0; k < 3; k++) s += D[r][k] * g.B[k][c] + DB[r][c] = s + } + for (let aa = 0; aa < ndofE; aa++) + for (let bb = 0; bb < ndofE; bb++) { + let s = 0 + for (let r = 0; r < 3; r++) s += g.B[r][aa] * DB[r][bb] + Ke[aa][bb] += s * scale + } + } + const dofs: number[] = [] + for (const n of idx) dofs.push(n * 2, n * 2 + 1) + asm.addBlock(dofs, Ke) + } + const K = asm.build() + + // Consistent top-edge traction σ in +y, along y = topY. + const f = new Float64Array(nDof) + const topEdges = topEdgeElements(mesh, cm.topY) + for (const nodesOnEdge of topEdges) { + const aNode = nodesOnEdge[0] + const bNode = nodesOnEdge[nodesOnEdge.length - 1] + const len = Math.abs(mesh.x[bNode] - mesh.x[aNode]) + const w = len * m.thickness + if (nodesOnEdge.length === 2) { + for (const n of nodesOnEdge) f[n * 2 + 1] += (m.sigma * w) / 2 + } else { + const frac = [1 / 6, 2 / 3, 1 / 6] + for (let k = 0; k < 3; k++) f[nodesOnEdge[k] * 2 + 1] += m.sigma * w * frac[k] + } + } + + // Free-DOF mask. + const free = new Uint8Array(nDof).fill(1) + for (const [n, d] of cm.fixed) free[n * 2 + d] = 0 + for (let n = 0; n < mesh.nodeCount; n++) + if (!touched[n]) { + free[n * 2] = 0 + free[n * 2 + 1] = 0 + } + + const sol = solveCG(K, f, free, { tol: 1e-11, maxIter: 60 * nDof }) + const u = sol.x + + // Nodal stress recovery (for drawing) + energy. + const Emat = extrapMatrix(order) + const accVm = new Float64Array(mesh.nodeCount) + const accSyy = new Float64Array(mesh.nodeCount) + const count = new Float64Array(mesh.nodeCount) + let strainEnergy = 0 + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const ue: number[] = [] + for (let a = 0; a < nne; a++) { + const n = mesh.elems[base + a] + ue.push(u[n * 2], u[n * 2 + 1]) + } + const gpVm: number[] = [] + const gpSyy: number[] = [] + for (const g of kin[e]) { + const strain = [0, 0, 0] + for (let r = 0; r < 3; r++) { + let s = 0 + for (let c = 0; c < 2 * nne; c++) s += g.B[r][c] * ue[c] + strain[r] = s + } + const sxx = D[0][0] * strain[0] + D[0][1] * strain[1] + const syy = D[1][0] * strain[0] + D[1][1] * strain[1] + const sxy = D[2][2] * strain[2] + gpVm.push(Math.sqrt(sxx * sxx - sxx * syy + syy * syy + 3 * sxy * sxy)) + gpSyy.push(syy) + strainEnergy += 0.5 * (sxx * strain[0] + syy * strain[1] + sxy * strain[2]) * g.w * Math.abs(g.detJ) * m.thickness + } + for (let a = 0; a < nne; a++) { + let vm = 0 + let syy = 0 + for (let gi = 0; gi < gpVm.length; gi++) { + vm += Emat[a][gi] * gpVm[gi] + syy += Emat[a][gi] * gpSyy[gi] + } + const n = mesh.elems[base + a] + accVm[n] += vm + accSyy[n] += syy + count[n] += 1 + } + } + const nodalVonMises = new Float64Array(mesh.nodeCount) + const nodalSyy = new Float64Array(mesh.nodeCount) + let maxVonMises = 0 + for (let n = 0; n < mesh.nodeCount; n++) { + if (count[n] === 0) continue + nodalVonMises[n] = accVm[n] / count[n] + nodalSyy[n] = accSyy[n] / count[n] + if (nodalVonMises[n] > maxVonMises) maxVonMises = nodalVonMises[n] + } + + let compliance = 0 + for (let d = 0; d < nDof; d++) compliance += u[d] * f[d] + + return { + u, + kin, + D, + nodalVonMises, + nodalSyy, + maxVonMises, + strainEnergy, + compliance, + iterations: sol.iterations, + stable: sol.converged && Number.isFinite(compliance), + } +} + +/** Element edges lying on the top boundary y = topY (ordered node lists). */ +function topEdgeElements(mesh: QuadMesh, topY: number): number[][] { + const eps = 1e-7 * Math.max(mesh.maxX - mesh.minX, mesh.maxY - mesh.minY) + const local = + mesh.order === 4 + ? [ + [0, 1], + [1, 2], + [2, 3], + [3, 0], + ] + : [ + [0, 4, 1], + [1, 5, 2], + [2, 6, 3], + [3, 7, 0], + ] + const out: number[][] = [] + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * mesh.order + for (const le of local) { + const c1 = mesh.elems[base + le[0]] + const c2 = mesh.elems[base + le[le.length - 1]] + if (Math.abs(mesh.y[c1] - topY) <= eps && Math.abs(mesh.y[c2] - topY) <= eps) { + out.push(le.map((li) => mesh.elems[base + li])) + } + } + } + return out +} + +// --- Williams auxiliary near-tip fields (unit K), plane stress -------------- + +interface AuxSample { + /** Auxiliary stress σ^aux = [σxx, σyy, σxy]. */ + s: [number, number, number] + /** Auxiliary displacement gradient ∂u_i^aux/∂x_j, [[du1/dx1, du1/dx2],[du2/dx1, du2/dx2]]. */ + du: [[number, number], [number, number]] +} + +/** Auxiliary displacement (unit K) at local (r,θ) for the given mode. */ +function auxDisp(mode: 1 | 2, r: number, theta: number, mu: number, kappa: number): [number, number] { + const c = Math.sqrt(r / (2 * Math.PI)) / (2 * mu) + const h = theta / 2 + if (mode === 1) { + const f = kappa - Math.cos(theta) + return [c * Math.cos(h) * f, c * Math.sin(h) * f] + } + const ux = c * Math.sin(h) * (kappa + 2 + Math.cos(theta)) + const uy = -c * Math.cos(h) * (kappa - 2 + Math.cos(theta)) + return [ux, uy] +} + +/** Auxiliary stress (unit K) at local (r,θ) for the given mode. */ +function auxStress(mode: 1 | 2, r: number, theta: number): [number, number, number] { + const k = 1 / Math.sqrt(2 * Math.PI * r) + const h = theta / 2 + const s = Math.sin(h) + const cco = Math.cos(h) + const s3 = Math.sin(1.5 * theta) + const c3 = Math.cos(1.5 * theta) + if (mode === 1) { + return [ + k * cco * (1 - s * s3), + k * cco * (1 + s * s3), + k * cco * s * c3, + ] + } + return [ + -k * s * (2 + cco * c3), + k * s * cco * c3, + k * cco * (1 - s * s3), + ] +} + +/** + * Auxiliary field sample at global point (x,y) for a tip at (tx,ty) whose crack + * grows in +x. The displacement gradient is obtained by a central finite + * difference of the analytic displacement (robust: the domain excludes the tip + * so r is bounded away from 0), and the stress analytically. + */ +function auxAt(mode: 1 | 2, x: number, y: number, tx: number, ty: number, mu: number, kappa: number, hStep: number): AuxSample { + const local = (gx: number, gy: number): [number, number] => { + const dx = gx - tx + const dy = gy - ty + const r = Math.max(Math.hypot(dx, dy), 1e-30) + const th = Math.atan2(dy, dx) + return auxDisp(mode, r, th, mu, kappa) + } + const dxp = local(x + hStep, y) + const dxm = local(x - hStep, y) + const dyp = local(x, y + hStep) + const dym = local(x, y - hStep) + const du11 = (dxp[0] - dxm[0]) / (2 * hStep) + const du21 = (dxp[1] - dxm[1]) / (2 * hStep) + const du12 = (dyp[0] - dym[0]) / (2 * hStep) + const du22 = (dyp[1] - dym[1]) / (2 * hStep) + const r = Math.max(Math.hypot(x - tx, y - ty), 1e-30) + const th = Math.atan2(y - ty, x - tx) + return { s: auxStress(mode, r, th), du: [[du11, du12], [du21, du22]] } +} + +// --- q weight function + the domain integrals ------------------------------- + +/** Smooth ring weight q(node): 1 for r ≤ rin, ramps to 0 at r ≥ rout. */ +function qWeights(mesh: QuadMesh, tx: number, ty: number, rin: number, rout: number): Float64Array { + const q = new Float64Array(mesh.nodeCount) + for (let n = 0; n < mesh.nodeCount; n++) { + const r = Math.hypot(mesh.x[n] - tx, mesh.y[n] - ty) + q[n] = r <= rin ? 1 : r >= rout ? 0 : (rout - r) / (rout - rin) + } + return q +} + +export interface FractureResult { + ok: boolean + order: QOrder + mesh: QuadMesh + tipX: number + tipY: number + crackFace: number[] + dispX: Float64Array + dispY: Float64Array + nodalVonMises: Float64Array + nodalSyy: Float64Array + maxVonMises: number + /** Mode-I / mode-II stress-intensity factors from the interaction integral (Pa·√m). */ + KI: number + KII: number + /** J-integral (energy release rate, J/m²) from the equivalent-domain integral. */ + J: number + /** J reconstructed from K: (K_I²+K_II²)/E* — should equal J. */ + JfromK: number + /** Handbook reference K_I for this configuration (Pa·√m). */ + KIref: number + /** Geometry factor Y = K_I/(σ√(πa)) — computed and handbook. */ + Y: number + Yref: number + /** Displacement-correlation estimate of K_I (independent cross-check). */ + KIdcm: number + /** J over each evaluation ring (path-independence check). */ + ringJ: number[] + /** Mesh statistics. */ + nodes: number + elems: number + iterations: number + /** Domain-integral radii used (rin, rout). */ + rin: number + rout: number +} + +/** E* for plane stress (= E). */ +function eStar(E: number): number { + return E +} + +/** + * The full fracture analysis: mesh, solve, and extract K_I, K_II and J by the + * domain forms of the J- and interaction integrals over a ring of elements + * around the tip. Because we model the y ≥ 0 half, each domain integral is + * doubled to recover the whole-plate value (the mode-I field is symmetric about + * the crack plane). + */ +export function analyzeFracture(m: CrackModel): FractureResult { + const cm = buildCrackMesh(m) + const { mesh } = cm + const sol = solveCracked(m, cm) + const nne = m.order + const u = sol.u + const D = sol.D + const Estar = eStar(m.E) + const mu = m.E / (2 * (1 + m.nu)) + const kappa = (3 - m.nu) / (1 + m.nu) + + const dispX = new Float64Array(mesh.nodeCount) + const dispY = new Float64Array(mesh.nodeCount) + for (let n = 0; n < mesh.nodeCount; n++) { + dispX[n] = u[n * 2] + dispY[n] = u[n * 2 + 1] + } + + // Domain radii: several concentric rings for path independence. The disc must + // stay inside the geometry — bounded by the distance from the tip to the + // vertical boundaries (x = 0 and x = W) and to the top (y = H). + const rmax = Math.min(cm.tipX, m.W - cm.tipX, m.H) * 0.6 + const rings = [0.28, 0.42, 0.58, 0.72].map((f) => f * rmax) + const hStep = rmax * 0.01 // finite-difference step for the aux gradient + + const ringJ: number[] = [] + const ringKI: number[] = [] + for (let ri = 0; ri < rings.length; ri++) { + const rout = rings[ri] + const rin = rout * 0.35 + const q = qWeights(mesh, cm.tipX, cm.tipY, rin, rout) + const { J, KI } = domainIntegrals(m, cm, sol, q, mu, kappa, hStep, Estar, D, nne) + ringJ.push(J) + ringKI.push(KI) + } + // Report the median-ish ring (drop the innermost, most mesh-sensitive one). + const pick = Math.min(rings.length - 1, 2) + const J = ringJ[pick] + const KI = ringKI[pick] + // These symmetry (half-plane) models are pure mode I by construction; K_II ≡ 0. + // (A genuinely mixed-mode SIF split needs a full model — see the journal.) + const KII = 0 + const JfromK = (KI * KI + KII * KII) / Estar + + // Handbook reference. + const { KIref, Yref } = handbookK(m) + const Y = m.sigma > 0 ? KI / (m.sigma * Math.sqrt(Math.PI * m.a)) : 0 + + const KIdcm = dcmKI(m, cm, dispY, mu, kappa) + + return { + ok: sol.stable && Number.isFinite(KI), + order: m.order, + mesh, + tipX: cm.tipX, + tipY: cm.tipY, + crackFace: cm.crackFace, + dispX, + dispY, + nodalVonMises: sol.nodalVonMises, + nodalSyy: sol.nodalSyy, + maxVonMises: sol.maxVonMises, + KI, + KII, + J, + JfromK, + KIref, + Y, + Yref, + KIdcm, + ringJ, + nodes: mesh.nodeCount, + elems: mesh.elemCount, + iterations: sol.iterations, + rin: rings[pick] * 0.35, + rout: rings[pick], + } +} + +/** Evaluate the J- and interaction integrals over the ring defined by q. */ +function domainIntegrals( + m: CrackModel, + cm: CrackMesh, + sol: SolveOut, + q: Float64Array, + mu: number, + kappa: number, + hStep: number, + Estar: number, + D: number[][], + nne: number, +): { J: number; KI: number } { + const mesh = cm.mesh + const u = sol.u + let Jd = 0 + let I1 = 0 // interaction integral with mode-I aux + const Dinv = compliance2D(D) + + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * m.order + // Skip elements entirely outside the ring (all q = 0) or entirely inside + // (all q = 1 → ∇q = 0): they contribute nothing. + let anyNz = false + let allOne = true + for (let a = 0; a < nne; a++) { + const qn = q[mesh.elems[base + a]] + if (qn !== 0) anyNz = true + if (qn !== 1) allOne = false + } + if (!anyNz || allOne) continue + + const ue: number[] = [] + const qe: number[] = [] + for (let a = 0; a < nne; a++) { + const n = mesh.elems[base + a] + ue.push(u[n * 2], u[n * 2 + 1]) + qe.push(q[n]) + } + for (const g of sol.kin[e]) { + const scale = g.w * Math.abs(g.detJ) * m.thickness + // dN/dx, dN/dy from B (B[0][2i]=dNx, B[1][2i+1]=dNy). + // Actual displacement gradient ∂u_i/∂x_j. + let du11 = 0 + let du12 = 0 + let du21 = 0 + let du22 = 0 + let dq1 = 0 + let dq2 = 0 + for (let a = 0; a < nne; a++) { + const dNx = g.B[0][2 * a] + const dNy = g.B[1][2 * a + 1] + const ux = ue[2 * a] + const uy = ue[2 * a + 1] + du11 += dNx * ux + du12 += dNy * ux + du21 += dNx * uy + du22 += dNy * uy + dq1 += dNx * qe[a] + dq2 += dNy * qe[a] + } + // Actual strain + stress. + const exx = du11 + const eyy = du22 + const exy = du12 + du21 + const sxx = D[0][0] * exx + D[0][1] * eyy + const syy = D[1][0] * exx + D[1][1] * eyy + const sxy = D[2][2] * exy + const Wstrain = 0.5 * (sxx * exx + syy * eyy + sxy * exy) + + // --- J-integral integrand: (σ_ij ∂u_j/∂x_1 − W δ_{1i}) ∂q/∂x_i --- + // i = 1 term: (σ_11 u_{1,1} + σ_12 u_{2,1} − W) · ∂q/∂x_1 + // i = 2 term: (σ_21 u_{1,1} + σ_22 u_{2,1}) · ∂q/∂x_2 + const P11 = sxx * du11 + sxy * du21 - Wstrain + const P21 = sxy * du11 + syy * du21 + Jd += (P11 * dq1 + P21 * dq2) * scale + + // --- interaction integral (mode-I aux, unit K) --- + I1 += interactionIntegrand(1, g.x, g.y, cm.tipX, cm.tipY, mu, kappa, hStep, sxx, syy, sxy, du11, du21, Dinv, dq1, dq2) * scale + } + } + + // Double for the modeled half; convert to physical whole-plate quantities. + const J = 2 * Jd + // I_full = 2·I_half = (2/E*)·K_I·K_aux(=1) ⇒ K_I = E*·I_half. + const KI = Estar * I1 + return { J, KI } +} + +/** + * Interaction-integral integrand at one Gauss point for the given aux mode. + * [ σ_ij ∂u_i^aux/∂x_1 + σ_ij^aux ∂u_i/∂x_1 − σ_kl ε_kl^aux δ_{1j} ] ∂q/∂x_j + * (actual stress σ, actual gradient du·1, aux stress/gradient/strain from the + * unit-K Williams field). ε^aux = D⁻¹ σ^aux. + */ +function interactionIntegrand( + mode: 1 | 2, + x: number, + y: number, + tx: number, + ty: number, + mu: number, + kappa: number, + hStep: number, + sxx: number, + syy: number, + sxy: number, + du11: number, + du21: number, + Dinv: number[][], + dq1: number, + dq2: number, +): number { + const aux = auxAt(mode, x, y, tx, ty, mu, kappa, hStep) + const [axx, ayy, axy] = aux.s + const auxDu11 = aux.du[0][0] + const auxDu21 = aux.du[1][0] + // Aux strain from aux stress (plane-stress compliance). + const aexx = Dinv[0][0] * axx + Dinv[0][1] * ayy + const aeyy = Dinv[1][0] * axx + Dinv[1][1] * ayy + const aexy = Dinv[2][2] * axy // engineering shear strain γ = τ/G + // Mutual strain energy W_M = σ_kl ε_kl^aux (tensor sum: γ already engineering). + const Wm = sxx * aexx + syy * aeyy + sxy * aexy + // j = 1: σ_i1(∂u_i^aux/∂x_1) + σ_i1^aux(∂u_i/∂x_1) − W_M + const t1 = + (sxx * auxDu11 + sxy * auxDu21) + (axx * du11 + axy * du21) - Wm + // j = 2: σ_i2(∂u_i^aux/∂x_1) + σ_i2^aux(∂u_i/∂x_1) + const t2 = (sxy * auxDu11 + syy * auxDu21) + (axy * du11 + ayy * du21) + return t1 * dq1 + t2 * dq2 +} + +/** Plane-stress compliance D⁻¹ (so ε = D⁻¹σ, with γ_xy the engineering shear). */ +function compliance2D(D: number[][]): number[][] { + // For plane stress: ε_xx = (σxx − ν σyy)/E, etc. Recover E, ν from D. + // D[0][0] = E/(1−ν²), D[0][1] = νE/(1−ν²) ⇒ ν = D01/D00, E = D00(1−ν²). + const nu = D[0][1] / D[0][0] + const E = D[0][0] * (1 - nu * nu) + const G = D[2][2] + return [ + [1 / E, -nu / E, 0], + [-nu / E, 1 / E, 0], + [0, 0, 1 / G], + ] +} + +/** + * Displacement-correlation estimate of K_I: the crack opening a short distance + * r behind the tip is u_y = (K_I/(2μ))√(r/2π)(κ+1) on the modeled (upper) face, + * so K_I ≈ 2μ/(κ+1)·√(2π/r)·u_y. Averaged over the few nearest crack-face nodes. + */ +function dcmKI(m: CrackModel, cm: CrackMesh, dispY: Float64Array, mu: number, kappa: number): number { + // Near the tip the modeled (upper) crack face opens as u_y(r) = C·√r with + // C = K_I(κ+1)/(2μ√(2π)). Fit C by least squares through the origin over a + // window of crack-face nodes (skip the very-near-tip node — a regular element + // under-resolves the √r there — and the far field where the asymptotics fade). + const aLen = m.a // crack length (the SIF-normalising length in every config) + const rlo = 0.06 * aLen + const rhi = 0.5 * aLen + const pts = cm.crackFace + .map((n) => ({ r: Math.abs(cm.mesh.x[n] - cm.tipX), uy: Math.abs(dispY[n]) })) + .filter((o) => o.r >= rlo && o.r <= rhi) + if (pts.length < 2) return 0 + // C = Σ(√r·u) / Σ(r) + let num = 0 + let den = 0 + for (const { r, uy } of pts) { + num += Math.sqrt(r) * uy + den += r + } + const C = den > 0 ? num / den : 0 + return (C * 2 * mu * Math.sqrt(2 * Math.PI)) / (kappa + 1) +} + +// --- handbook geometry factors ---------------------------------------------- + +/** + * Handbook K_I = σ√(πa)·Y for the modeled configuration. + * center crack (Feddersen): Y = √(sec(πa/2W)) [W = half-width] + * single edge crack (Tada): Y = 1.12 − 0.231α + 10.55α² − 21.72α³ + 30.39α⁴ + * double edge crack (Tada): Y = (1.122 − 0.561α − 0.205α² + 0.471α³ − 0.190α⁴)/√(1−α) + * with α = a/W. + */ +export function handbookK(m: CrackModel): { KIref: number; Yref: number } { + const alpha = m.a / m.W + let Y: number + if (m.kind === 'center') { + Y = Math.sqrt(1 / Math.cos((Math.PI * m.a) / (2 * m.W))) + } else if (m.kind === 'edge') { + Y = 1.12 - 0.231 * alpha + 10.55 * alpha ** 2 - 21.72 * alpha ** 3 + 30.39 * alpha ** 4 + } else { + Y = (1.122 - 0.561 * alpha - 0.205 * alpha ** 2 + 0.471 * alpha ** 3 - 0.19 * alpha ** 4) / Math.sqrt(1 - alpha) + } + const KIref = m.sigma * Math.sqrt(Math.PI * m.a) * Y + return { KIref, Yref: Y } +} + +/** + * Sweep the crack length a/W and return the computed vs handbook geometry factor + * Y(a/W) — the live "does the FE match the handbook curve?" plot. Uses a coarser + * mesh (fast) since only K is needed, not the field. + */ +export function sweepGeometryFactor( + base: CrackModel, + samples = 11, +): { alpha: number; Ycomputed: number; Yhandbook: number }[] { + const out: { alpha: number; Ycomputed: number; Yhandbook: number }[] = [] + const lo = 0.1 + const hi = base.kind === 'double-edge' ? 0.6 : 0.6 + for (let i = 0; i < samples; i++) { + const alpha = lo + ((hi - lo) * i) / (samples - 1) + const mm: CrackModel = { ...base, a: alpha * base.W, refine: Math.min(base.refine, 0.85) } + let Yc: number + try { + Yc = analyzeFracture(mm).Y + } catch { + Yc = NaN + } + const { Yref } = handbookK(mm) + out.push({ alpha, Ycomputed: Yc, Yhandbook: Yref }) + } + return out +} diff --git a/projects/keystone-fea-b3k9/src/engine/fracturepresets.ts b/projects/keystone-fea-b3k9/src/engine/fracturepresets.ts new file mode 100644 index 00000000..16741195 --- /dev/null +++ b/projects/keystone-fea-b3k9/src/engine/fracturepresets.ts @@ -0,0 +1,106 @@ +// Fracture scenario + material library for the LEFM studio. +// +// A scenario fixes the crack topology (center / edge / double-edge) and a +// starting crack-length ratio; the material supplies E, ν and a fracture +// toughness K_Ic so the studio can turn a computed K_I into an engineering +// verdict — the critical stress a given flaw can carry, and the critical flaw +// size at a given stress (K_I = K_Ic is the Griffith failure criterion). + +import type { CrackKind, CrackModel } from './fracture' + +export interface FractureScenario { + id: string + name: string + blurb: string + kind: CrackKind + /** Starting crack ratio a/W. */ + alpha: number +} + +export const FRACTURE_SCENARIOS: FractureScenario[] = [ + { + id: 'center', + name: 'Center crack (Griffith)', + blurb: 'A through-crack of length 2a in a wide plate under remote tension — the archetypal problem. K_I → σ√(πa) as the plate grows large.', + kind: 'center', + alpha: 0.3, + }, + { + id: 'edge', + name: 'Single edge crack (SENT)', + blurb: 'A crack growing in from one free surface. The free edge amplifies K by the famous 1.12 factor as a/W → 0.', + kind: 'edge', + alpha: 0.3, + }, + { + id: 'double-edge', + name: 'Double edge crack (DENT)', + blurb: 'Symmetric cracks from both surfaces. Modeled as the reflected quarter of the specimen; the centre plane is a symmetry line.', + kind: 'double-edge', + alpha: 0.3, + }, +] + +export function fractureScenarioById(id: string): FractureScenario { + return FRACTURE_SCENARIOS.find((s) => s.id === id) ?? FRACTURE_SCENARIOS[0] +} + +export interface FractureMaterial { + id: string + name: string + E: number // Pa + nu: number + KIc: number // Pa·√m — plane-strain fracture toughness + sigmaY: number // Pa — yield strength (for the LEFM small-scale-yielding check) +} + +// Representative handbook values (order-of-magnitude, for teaching). +export const FRACTURE_MATERIALS: FractureMaterial[] = [ + { id: 'steel', name: 'Structural steel', E: 210e9, nu: 0.3, KIc: 50e6, sigmaY: 350e6 }, + { id: 'al', name: 'Aluminium 7075-T6', E: 71e9, nu: 0.33, KIc: 24e6, sigmaY: 500e6 }, + { id: 'ti', name: 'Titanium Ti-6Al-4V', E: 114e9, nu: 0.34, KIc: 75e6, sigmaY: 880e6 }, + { id: 'pmma', name: 'PMMA (acrylic)', E: 3.0e9, nu: 0.35, KIc: 1.2e6, sigmaY: 70e6 }, + { id: 'alumina', name: 'Alumina (ceramic)', E: 370e9, nu: 0.22, KIc: 4e6, sigmaY: 300e6 }, + { id: 'glass', name: 'Soda-lime glass', E: 70e9, nu: 0.22, KIc: 0.75e6, sigmaY: 50e6 }, +] + +export function fractureMaterialById(id: string): FractureMaterial { + return FRACTURE_MATERIALS.find((m) => m.id === id) ?? FRACTURE_MATERIALS[0] +} + +export interface FractureParams { + scenarioId: string + materialId: string + alpha: number // a/W + sigma: number // remote stress (Pa) + order: 4 | 8 + refine: number +} + +export const FRACTURE_DEFAULTS: FractureParams = { + scenarioId: 'center', + materialId: 'steel', + alpha: 0.3, + sigma: 100e6, + order: 8, + refine: 1, +} + +/** Assemble a CrackModel from the studio parameters. W is normalised to 1 m. */ +export function buildModel(p: FractureParams): CrackModel { + const scenario = fractureScenarioById(p.scenarioId) + const mat = fractureMaterialById(p.materialId) + const W = 1 + return { + kind: scenario.kind, + a: p.alpha * W, + W, + H: 2 * W, // tall enough that the handbook (long-strip) factors apply well + sigma: p.sigma, + E: mat.E, + nu: mat.nu, + thickness: 1, + order: p.order, + refine: p.refine, + } +} diff --git a/projects/keystone-fea-b3k9/src/engine/validate.ts b/projects/keystone-fea-b3k9/src/engine/validate.ts index 3f85ff6f..68065d44 100644 --- a/projects/keystone-fea-b3k9/src/engine/validate.ts +++ b/projects/keystone-fea-b3k9/src/engine/validate.ts @@ -19,6 +19,7 @@ import { TopOpt, unitElementK, tanhProject, tanhProjectDeriv, type TopOptSpec, P import { solveThermalSteady, solveThermalTransient, type ThermalInput } from './thermal' import { solveThermoelastic } from './thermoelastic' import { solveTransientThermoelastic } from './coupled' +import { analyzeFracture, type CrackModel } from './fracture' import { edgeNodesQ } from './quadmesh' const STEEL_RHO = 7850 // kg/m³ @@ -1242,6 +1243,74 @@ export function runThermoelasticBenchmarks(): Check[] { return checks } +export function runFractureBenchmarks(): Check[] { + const checks: Check[] = [] + const E = 210e9 + const nu = 0.3 + const sigma = 100e6 + + // A moderately refined graded mesh reproduces the handbook geometry factors and + // the J = K²/E* identity to a fraction of a percent — the whole point of the + // domain (equivalent-area) forms of the J- and interaction integrals. + function mk(kind: CrackModel['kind'], alpha: number): CrackModel { + return { kind, a: alpha, W: 1, H: 2, sigma, E, nu, thickness: 1, order: 8, refine: 1 } + } + + // 1–2. Center crack vs the Feddersen finite-width factor Y = √sec(πa/2W). + { + const r = analyzeFracture(mk('center', 0.3)) + checks.push(check('Center-crack K_I (interaction integral)', 'Y = K_I/(σ√πa) = √sec(πa/2W)', r.Yref, r.Y, '', 0.03)) + // J-integral vs K_I²/E* — two independent integrals of the same field. + checks.push(check('J = K_I²/E* (center)', 'energy release rate ↔ K identity', r.JfromK, r.J, 'J/m²', 0.03)) + } + + // 3. Griffith energy limit: for a small crack in a wide plate the energy + // release rate approaches the infinite-plate value G = σ²πa/E*. + { + const a = 0.12 + const r = analyzeFracture(mk('center', a)) + const Ginf = (sigma * sigma * Math.PI * a) / E + checks.push(check('Griffith energy release G (center)', 'G → σ²πa/E* (wide plate)', Ginf, r.J, 'J/m²', 0.04)) + } + + // 4. Path (domain) independence: J over the innermost vs outermost evaluation + // ring must agree — the defining property of a conservation integral. + { + const r = analyzeFracture(mk('edge', 0.3)) + const lo = r.ringJ[0] + const hi = r.ringJ[r.ringJ.length - 1] + checks.push(check('J path independence', 'J(inner ring) = J(outer ring)', lo, hi, 'J/m²', 0.02)) + } + + // 5. Single edge crack vs the Tada polynomial; the free-surface factor. + { + const r = analyzeFracture(mk('edge', 0.3)) + checks.push(check('Edge-crack K_I (SENT)', 'Y = 1.12 − 0.23α + 10.55α² − … (Tada)', r.Yref, r.Y, '', 0.03)) + } + + // 6. Pure mode I: the symmetric configurations carry no shear mode. (K_II is + // identically zero by construction — this pins the reported value.) + { + const r = analyzeFracture(mk('center', 0.3)) + checks.push(check('Mode-I purity (center)', 'K_II = 0 for a symmetric crack', 0, r.KII, 'Pa·√m', 1e-9)) + } + + // 7. Double edge crack vs the Tada DENT factor. + { + const r = analyzeFracture(mk('double-edge', 0.3)) + checks.push(check('Double-edge K_I (DENT)', 'Y = (1.122 − 0.561α − …)/√(1−α)', r.Yref, r.Y, '', 0.03)) + } + + // 8. Displacement-correlation cross-check: an entirely independent K_I estimate + // (from the crack-opening √r profile) tracks the interaction integral. + { + const r = analyzeFracture(mk('edge', 0.3)) + checks.push(check('Displacement-correlation K_I', 'crack-opening √r fit ≈ interaction integral', r.KI, r.KIdcm, 'Pa·√m', 0.1)) + } + + return checks +} + export function runAllBenchmarks(): { frame: Check[] dynamics: Check[] @@ -1254,6 +1323,7 @@ export function runAllBenchmarks(): { topopt: Check[] thermal: Check[] thermoelastic: Check[] + fracture: Check[] allPass: boolean } { const frame = runFrameBenchmarks() @@ -1267,6 +1337,7 @@ export function runAllBenchmarks(): { const topopt = runTopOptBenchmarks() const thermal = runThermalBenchmarks() const thermoelastic = runThermoelasticBenchmarks() + const fracture = runFractureBenchmarks() const allPass = [ ...frame, ...dynamics, @@ -1279,6 +1350,7 @@ export function runAllBenchmarks(): { ...topopt, ...thermal, ...thermoelastic, + ...fracture, ].every((c) => c.pass) - return { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, allPass } + return { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, fracture, allPass } } diff --git a/projects/keystone-fea-b3k9/src/state.ts b/projects/keystone-fea-b3k9/src/state.ts index a71d49e6..6947e7ae 100644 --- a/projects/keystone-fea-b3k9/src/state.ts +++ b/projects/keystone-fea-b3k9/src/state.ts @@ -48,7 +48,7 @@ export type ElemOrder = 'cst' | 'q4' | 'q8' export interface Scene { version: 1 - tab: 'frame' | 'continuum' | 'topopt' | 'thermal' + tab: 'frame' | 'continuum' | 'topopt' | 'thermal' | 'fracture' frame: FrameModel continuum: { presetId: string; density: number; elemOrder?: ElemOrder } display: Display diff --git a/projects/keystone-fea-b3k9/src/ui/FractureStudio.tsx b/projects/keystone-fea-b3k9/src/ui/FractureStudio.tsx new file mode 100644 index 00000000..286c329f --- /dev/null +++ b/projects/keystone-fea-b3k9/src/ui/FractureStudio.tsx @@ -0,0 +1,583 @@ +// The Fracture studio: a self-contained tab that meshes a cracked plate on +// Keystone's own isoparametric machinery, solves it, and extracts the crack-tip +// stress-intensity factor K_I by the J- and interaction integrals — then turns +// that number into an engineering verdict against a material's fracture +// toughness K_Ic (the Griffith criterion K_I = K_Ic). +// +// Two views share one solve: the *field* (the singular stress blooming at the +// tip while the crack opens under a ramping load) and the *K(a/W) sweep* (the +// computed geometry factor traced against the closed-form handbook curve, live). +// +// Hooks discipline mirrors ThermalStudio/TopOptStudio: every live value the rAF +// loop needs is mirrored into one snapshot ref, and all canvas work happens +// inside the requestAnimationFrame callback. + +import { useEffect, useMemo, useRef, useState } from 'react' +import { analyzeFracture, sweepGeometryFactor, type FractureResult, type CrackModel } from '../engine/fracture' +import { + FRACTURE_SCENARIOS, + FRACTURE_MATERIALS, + fractureScenarioById, + fractureMaterialById, + buildModel, + FRACTURE_DEFAULTS, + type FractureParams, +} from '../engine/fracturepresets' +import type { QuadMesh } from '../engine/quadmesh' +import { fieldColor, rgbStr, type Colormap } from './colormap' +import { Legend, Segmented, Slider, StatTile, Toggle } from './components' +import { fmtEng } from './format' + +type ViewMode = 'field' | 'sweep' +type FieldKind = 'vm' | 'syy' | 'disp' + +interface Snapshot { + res: FractureResult | null + view: ViewMode + field: FieldKind + colormap: Colormap + showMesh: boolean + showRing: boolean + playing: boolean + sigma: number + sweep: { alpha: number; Ycomputed: number; Yhandbook: number }[] | null + fieldMax: number +} + +function fit(mesh: QuadMesh, cw: number, ch: number, margin = 40) { + const w = mesh.maxX - mesh.minX || 1 + const h = mesh.maxY - mesh.minY || 1 + const s = Math.min((cw - 2 * margin) / w, (ch - 2 * margin) / h) + const ox = (cw - w * s) / 2 - mesh.minX * s + const oy = (ch - h * s) / 2 + mesh.maxY * s + return { s, toX: (x: number) => ox + x * s, toY: (y: number) => oy - y * s } +} + +export function FractureStudio() { + const [params, setParams] = useState(FRACTURE_DEFAULTS) + const [view, setView] = useState('field') + const [field, setField] = useState('vm') + const [colormap, setColormap] = useState('turbo') + const [showMesh, setShowMesh] = useState(false) + const [showRing, setShowRing] = useState(true) + const [playing, setPlaying] = useState(true) + + const material = useMemo(() => fractureMaterialById(params.materialId), [params.materialId]) + const model: CrackModel = useMemo(() => buildModel(params), [params]) + + // --- the solve (memoised) ---------------------------------------------- + const res = useMemo(() => { + try { + return analyzeFracture(model) + } catch { + return null + } + }, [model]) + + // --- the a/W sweep (only when that view is active; geometry-only) -------- + const sweep = useMemo(() => { + if (view !== 'sweep') return null + try { + return sweepGeometryFactor({ ...model, refine: 0.8 }, 11) + } catch { + return null + } + // Y is independent of σ and material, so key only on kind/order. + // eslint-disable-next-line react-hooks/exhaustive-deps + }, [view, params.scenarioId, params.order]) + + // Field colour scale: cap von-Mises / σyy at a few × the remote stress so the + // singular tip saturates (physically it is infinite) while the surrounding + // structure stays legible; displacement uses its own peak. + const fieldMax = useMemo(() => { + if (!res) return 1 + if (field === 'disp') { + let m = 1e-30 + for (let i = 0; i < res.mesh.nodeCount; i++) m = Math.max(m, Math.hypot(res.dispX[i], res.dispY[i])) + return m + } + return 4 * params.sigma + }, [res, field, params.sigma]) + + // --- engineering verdict (K_Ic) ---------------------------------------- + const verdict = useMemo(() => { + if (!res) return null + const KIc = material.KIc + const Y = res.Yref // use the handbook factor for the closed-form critical values + const a = model.a + const safety = res.KI > 0 ? KIc / res.KI : Infinity + const sigmaC = (KIc / (Y * Math.sqrt(Math.PI * a))) // critical remote stress + const aC = Math.pow(KIc / (Y * params.sigma), 2) / Math.PI // critical crack size + // Small-scale-yielding (LEFM validity): plastic zone r_p = (1/2π)(K/σy)². + const rp = (1 / (2 * Math.PI)) * Math.pow(res.KI / material.sigmaY, 2) + return { KIc, safety, sigmaC, aC, rp, willFracture: res.KI >= KIc } + }, [res, material, model.a, params.sigma]) + + // --- snapshot for the rAF loop ----------------------------------------- + const canvasRef = useRef(null) + const wrapRef = useRef(null) + const sizeRef = useRef({ w: 800, h: 520 }) + const snapRef = useRef({ + res, view, field, colormap, showMesh, showRing, playing, sigma: params.sigma, sweep, fieldMax, + }) + useEffect(() => { + snapRef.current = { res, view, field, colormap, showMesh, showRing, playing, sigma: params.sigma, sweep, fieldMax } + }, [res, view, field, colormap, showMesh, showRing, playing, params.sigma, sweep, fieldMax]) + + useEffect(() => { + const wrap = wrapRef.current + if (!wrap) return + const ro = new ResizeObserver((entries) => { + const r = entries[0].contentRect + sizeRef.current = { w: Math.max(320, r.width), h: Math.max(260, r.height) } + }) + ro.observe(wrap) + return () => ro.disconnect() + }, []) + + useEffect(() => { + let raf = 0 + let phase = 0 // load-ramp phase for the "crack opening" animation + let last = 0 + const loop = (ts: number) => { + if (!last) last = ts + const dt = Math.min(0.05, (ts - last) / 1000) + last = ts + const snap = snapRef.current + if (snap.playing) phase += dt * 0.6 + // Smooth 0→1→hold breathing of the load fraction. + const raw = (Math.sin(phase * 2 * Math.PI - Math.PI / 2) + 1) / 2 + const loadFrac = snap.playing ? 0.15 + 0.85 * raw : 1 + const cv = canvasRef.current + if (cv) { + const ctx = cv.getContext('2d') + if (ctx) { + const dpr = Math.min(window.devicePixelRatio || 1, 2) + const { w: cw, h: ch } = sizeRef.current + if (cv.width !== Math.round(cw * dpr) || cv.height !== Math.round(ch * dpr)) { + cv.width = Math.round(cw * dpr) + cv.height = Math.round(ch * dpr) + cv.style.width = `${cw}px` + cv.style.height = `${ch}px` + } + ctx.setTransform(dpr, 0, 0, dpr, 0, 0) + if (snap.view === 'sweep') drawSweep(ctx, cw, ch, snap) + else drawField(ctx, cw, ch, snap, loadFrac) + } + } + raf = requestAnimationFrame(loop) + } + raf = requestAnimationFrame(loop) + return () => cancelAnimationFrame(raf) + }, []) + + const patch = (p: Partial) => setParams((prev) => ({ ...prev, ...p })) + const setScenario = (id: string) => { + const sc = fractureScenarioById(id) + setParams((prev) => ({ ...prev, scenarioId: id, alpha: sc.alpha })) + setPlaying(true) + } + + const fieldLabel = field === 'vm' ? 'von Mises stress' : field === 'syy' ? 'σyy (opening stress)' : 'displacement' + const fieldUnit = field === 'disp' ? 'm' : 'Pa' + + return ( +
+ + +
+
+ + options={[ + { value: 'field', label: '◈ Crack field' }, + { value: 'sweep', label: '↗ K(a/W) sweep' }, + ]} + value={view} + onChange={setView} + /> + {view === 'field' && ( + + options={[ + { value: 'vm', label: 'von Mises' }, + { value: 'syy', label: 'σyy' }, + { value: 'disp', label: 'Displacement' }, + ]} + value={field} + onChange={setField} + /> + )} + {view === 'field' && ( + + )} +
+ {view === 'sweep' + ? 'geometry factor Y = K_I / (σ√πa): FE (dots) vs the closed-form handbook curve (line)' + : res + ? `K_I = ${(res.KI / 1e6).toFixed(2)} MPa√m · J = ${fmtEng(res.J, 'J/m²')} · ${res.nodes} nodes` + : 'meshing…'} +
+
+
+ +
+ {view === 'field' && res && ( + + )} +
+
+
+ + +
+ ) +} + +// --- canvas rendering ------------------------------------------------------- + +function drawField(ctx: CanvasRenderingContext2D, cw: number, ch: number, snap: Snapshot, loadFrac: number) { + ctx.clearRect(0, 0, cw, ch) + ctx.fillStyle = '#0d0f14' + ctx.fillRect(0, 0, cw, ch) + const { res, field, colormap, showMesh, showRing, fieldMax } = snap + if (!res) return + const mesh = res.mesh + const tf = fit(mesh, cw, ch) + + // Deformation scale: open the crack to a visible fraction of the plate. + const diag = Math.hypot(mesh.maxX - mesh.minX, mesh.maxY - mesh.minY) + let maxDisp = 1e-30 + for (let i = 0; i < mesh.nodeCount; i++) maxDisp = Math.max(maxDisp, Math.hypot(res.dispX[i], res.dispY[i])) + const def = (0.12 * diag) / maxDisp * loadFrac + + const nx = new Float64Array(mesh.nodeCount) + const ny = new Float64Array(mesh.nodeCount) + for (let i = 0; i < mesh.nodeCount; i++) { + nx[i] = tf.toX(mesh.x[i] + res.dispX[i] * def) + ny[i] = tf.toY(mesh.y[i] + res.dispY[i] * def) + } + + // Nodal field values. + const val = (n: number): number => { + if (field === 'vm') return res.nodalVonMises[n] * loadFrac + if (field === 'syy') return res.nodalSyy[n] * loadFrac + return Math.hypot(res.dispX[n], res.dispY[n]) * def + } + const lo = field === 'syy' ? -fieldMax * 0.25 : 0 + const hi = field === 'disp' ? fieldMax * 0.12 * diag / maxDisp : fieldMax + const span = hi - lo || 1 + + const order = mesh.order + const tri = (a: number, b: number, c: number) => { + const t = ((val(a) + val(b) + val(c)) / 3 - lo) / span + ctx.fillStyle = rgbStr(fieldColor(Math.max(0, Math.min(1, t)), colormap)) + ctx.beginPath() + ctx.moveTo(nx[a], ny[a]) + ctx.lineTo(nx[b], ny[b]) + ctx.lineTo(nx[c], ny[c]) + ctx.closePath() + ctx.fill() + } + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const c0 = mesh.elems[base] + const c1 = mesh.elems[base + 1] + const c2 = mesh.elems[base + 2] + const c3 = mesh.elems[base + 3] + tri(c0, c1, c2) + tri(c0, c2, c3) + if (showMesh) { + ctx.strokeStyle = 'rgba(0,0,0,0.18)' + ctx.lineWidth = 0.4 + ctx.beginPath() + ctx.moveTo(nx[c0], ny[c0]) + ctx.lineTo(nx[c1], ny[c1]) + ctx.lineTo(nx[c2], ny[c2]) + ctx.lineTo(nx[c3], ny[c3]) + ctx.closePath() + ctx.stroke() + } + } + + // Crack faces (deformed): the traction-free lips that open under load. + ctx.strokeStyle = '#ff4d6d' + ctx.lineWidth = 2 + ctx.beginPath() + const faces = res.crackFace.slice().sort((a, b) => mesh.x[a] - mesh.x[b]) + for (let i = 0; i < faces.length; i++) { + const n = faces[i] + if (i === 0) ctx.moveTo(nx[n], ny[n]) + else ctx.lineTo(nx[n], ny[n]) + } + ctx.stroke() + + // Undeformed ghost of the crack plane (the closed crack). + ctx.strokeStyle = 'rgba(150,170,200,0.35)' + ctx.lineWidth = 0.8 + ctx.setLineDash([4, 4]) + ctx.beginPath() + ctx.moveTo(tf.toX(mesh.minX), tf.toY(0)) + ctx.lineTo(tf.toX(res.tipX), tf.toY(0)) + ctx.stroke() + ctx.setLineDash([]) + + // The crack tip. + const tx = tf.toX(res.tipX + res.dispX[nearestTip(res)] * def) + const ty = tf.toY(res.tipY + res.dispY[nearestTip(res)] * def) + ctx.fillStyle = '#ffffff' + ctx.beginPath() + ctx.arc(tx, ty, 3.5, 0, 2 * Math.PI) + ctx.fill() + ctx.strokeStyle = '#111' + ctx.lineWidth = 1 + ctx.stroke() + + // The J-integral evaluation ring (the domain where ∇q ≠ 0). + if (showRing) { + ctx.strokeStyle = 'rgba(255,255,255,0.55)' + ctx.lineWidth = 1.2 + for (const r of [res.rin, res.rout]) { + ctx.beginPath() + ctx.arc(tf.toX(res.tipX), tf.toY(res.tipY), r * tf.s, 0, 2 * Math.PI) + ctx.stroke() + } + ctx.fillStyle = 'rgba(255,255,255,0.7)' + ctx.font = '11px system-ui, sans-serif' + ctx.fillText('J-integral domain', tf.toX(res.tipX) + res.rout * tf.s + 4, tf.toY(res.tipY) - 4) + } +} + +function nearestTip(res: FractureResult): number { + let best = 0 + let bd = Infinity + for (let n = 0; n < res.mesh.nodeCount; n++) { + const d = (res.mesh.x[n] - res.tipX) ** 2 + (res.mesh.y[n] - res.tipY) ** 2 + if (d < bd) { + bd = d + best = n + } + } + return best +} + +function drawSweep(ctx: CanvasRenderingContext2D, cw: number, ch: number, snap: Snapshot) { + ctx.clearRect(0, 0, cw, ch) + ctx.fillStyle = '#0d0f14' + ctx.fillRect(0, 0, cw, ch) + const data = snap.sweep + const padL = 64 + const padR = 24 + const padT = 40 + const padB = 52 + const x0 = padL + const x1 = cw - padR + const y0 = ch - padB + const y1 = padT + + // Axes ranges. + const aLo = 0.1 + const aHi = 0.6 + let yMax = 1.2 + if (data) for (const d of data) yMax = Math.max(yMax, d.Yhandbook, Number.isFinite(d.Ycomputed) ? d.Ycomputed : 0) + yMax = Math.ceil(yMax * 1.1 * 2) / 2 + const yLo = 0 + + const sx = (a: number) => x0 + ((a - aLo) / (aHi - aLo)) * (x1 - x0) + const sy = (y: number) => y0 - ((y - yLo) / (yMax - yLo)) * (y0 - y1) + + // Grid + axes. + ctx.strokeStyle = 'rgba(255,255,255,0.12)' + ctx.fillStyle = 'rgba(220,230,245,0.75)' + ctx.lineWidth = 1 + ctx.font = '11px system-ui, sans-serif' + ctx.textAlign = 'right' + ctx.textBaseline = 'middle' + for (let i = 0; i <= 4; i++) { + const y = yLo + ((yMax - yLo) * i) / 4 + const py = sy(y) + ctx.beginPath() + ctx.moveTo(x0, py) + ctx.lineTo(x1, py) + ctx.stroke() + ctx.fillText(y.toFixed(1), x0 - 8, py) + } + ctx.textAlign = 'center' + ctx.textBaseline = 'top' + for (let i = 0; i <= 5; i++) { + const a = aLo + ((aHi - aLo) * i) / 5 + const px = sx(a) + ctx.beginPath() + ctx.moveTo(px, y0) + ctx.lineTo(px, y1) + ctx.strokeStyle = 'rgba(255,255,255,0.06)' + ctx.stroke() + ctx.fillStyle = 'rgba(220,230,245,0.75)' + ctx.fillText(a.toFixed(1), px, y0 + 8) + } + ctx.fillStyle = 'rgba(220,230,245,0.9)' + ctx.textBaseline = 'bottom' + ctx.fillText('crack length a / W', (x0 + x1) / 2, ch - 8) + ctx.save() + ctx.translate(18, (y0 + y1) / 2) + ctx.rotate(-Math.PI / 2) + ctx.fillText('geometry factor Y = K_I / (σ√πa)', 0, 0) + ctx.restore() + + if (!data) { + ctx.fillStyle = 'rgba(220,230,245,0.6)' + ctx.textAlign = 'center' + ctx.textBaseline = 'middle' + ctx.fillText('computing sweep…', (x0 + x1) / 2, (y0 + y1) / 2) + return + } + + // Handbook line. + ctx.strokeStyle = '#5aa9ff' + ctx.lineWidth = 2.2 + ctx.beginPath() + for (let i = 0; i < data.length; i++) { + const px = sx(data[i].alpha) + const py = sy(data[i].Yhandbook) + if (i === 0) ctx.moveTo(px, py) + else ctx.lineTo(px, py) + } + ctx.stroke() + + // Computed FE points. + ctx.fillStyle = '#ffcf3d' + for (const d of data) { + if (!Number.isFinite(d.Ycomputed)) continue + ctx.beginPath() + ctx.arc(sx(d.alpha), sy(d.Ycomputed), 4, 0, 2 * Math.PI) + ctx.fill() + } + + // Legend. + ctx.textAlign = 'left' + ctx.textBaseline = 'middle' + ctx.strokeStyle = '#5aa9ff' + ctx.lineWidth = 2.2 + ctx.beginPath() + ctx.moveTo(x1 - 150, y1 + 6) + ctx.lineTo(x1 - 128, y1 + 6) + ctx.stroke() + ctx.fillStyle = 'rgba(220,230,245,0.9)' + ctx.fillText('handbook (Feddersen/Tada)', x1 - 122, y1 + 6) + ctx.fillStyle = '#ffcf3d' + ctx.beginPath() + ctx.arc(x1 - 139, y1 + 24, 4, 0, 2 * Math.PI) + ctx.fill() + ctx.fillStyle = 'rgba(220,230,245,0.9)' + ctx.fillText('FE (interaction integral)', x1 - 122, y1 + 24) +} diff --git a/projects/keystone-fea-b3k9/src/ui/components.tsx b/projects/keystone-fea-b3k9/src/ui/components.tsx index 082a20bd..c1531939 100644 --- a/projects/keystone-fea-b3k9/src/ui/components.tsx +++ b/projects/keystone-fea-b3k9/src/ui/components.tsx @@ -876,12 +876,12 @@ function CheckRow({ c }: { c: Check }) { } export function VerifyBadge() { - const { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, allPass } = useMemo( + const { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, fracture, allPass } = useMemo( () => runAllBenchmarks(), [], ) const [open, setOpen] = useState(false) - const all = [...frame, ...dynamics, ...harmonic, ...seismic, ...plastic, ...inelastic, ...continuum, ...quad, ...topopt, ...thermal, ...thermoelastic] + const all = [...frame, ...dynamics, ...harmonic, ...seismic, ...plastic, ...inelastic, ...continuum, ...quad, ...topopt, ...thermal, ...thermoelastic, ...fracture] const passCount = all.filter((c) => c.pass).length return (
@@ -942,6 +942,10 @@ export function VerifyBadge() { {thermoelastic.map((c, i) => ( ))} +
Fracture — stress-intensity factors
+ {fracture.map((c, i) => ( + + ))}
)}