From 79944663d96fbbadf24237c1dc9707491eb8bbae Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 24 Jul 2026 03:39:21 +0000 Subject: [PATCH] Keystone v14: continuum linear buckling & stability studio MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the elastic-instability chapter to the Q4/Q8 continuum. A new from-scratch engine assembles the geometric (initial-stress) stiffness K_g from a reference plane-stress field and solves the generalized eigenproblem (K + λK_g)φ = 0 for the lowest critical load factors and buckled mode shapes. Because K_g is indefinite, the v9 SPD-mass subspace solver can't be reused; the eigensolver runs a K-orthonormal subspace iteration on A = K⁻¹(−K_g) (self-adjoint in the K-inner-product) with a standard symmetric Rayleigh–Ritz reduction and relative rank-revealing deflation that never trips on the indefinite/rank-deficient stress stiffness. - engine/buckling.ts: K, reference solve, stress stiffness, eigensolver, Euler reference load - engine/bucklingpresets.ts: materials + fixed-free / fixed-pinned / panel scenarios with section properties for the column curve - ui/BucklingStudio.tsx: sixth studio tab — animated buckled shape over axial load-path shading, and the Euler column curve (σ_cr=π²E/λ² with yield cut-off + FE design point); Q8/Q4 shear-locking toggle - 3 new live benchmarks (81 → 84): fixed-free Euler column to 0.05%, fixed-pinned to 0.5%, cantilever mode ladder λ₂/λ₁ = 9 to 0.25% Scope + conformance + lint + build all green; verified end-to-end in headless Chromium with 84/84 benchmarks and zero runtime errors. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01RGnTxnxKH7Fmj125T54jqg --- projects/keystone-fea-b3k9/JOURNAL.md | 137 +++- projects/keystone-fea-b3k9/project.json | 4 +- projects/keystone-fea-b3k9/src/App.css | 18 + projects/keystone-fea-b3k9/src/App.tsx | 6 +- .../keystone-fea-b3k9/src/engine/buckling.ts | 647 +++++++++++++++++ .../src/engine/bucklingpresets.ts | 198 ++++++ .../keystone-fea-b3k9/src/engine/validate.ts | 69 +- projects/keystone-fea-b3k9/src/state.ts | 2 +- .../src/ui/BucklingStudio.tsx | 670 ++++++++++++++++++ .../keystone-fea-b3k9/src/ui/components.tsx | 8 +- 10 files changed, 1749 insertions(+), 10 deletions(-) create mode 100644 projects/keystone-fea-b3k9/src/engine/buckling.ts create mode 100644 projects/keystone-fea-b3k9/src/engine/bucklingpresets.ts create mode 100644 projects/keystone-fea-b3k9/src/ui/BucklingStudio.tsx diff --git a/projects/keystone-fea-b3k9/JOURNAL.md b/projects/keystone-fea-b3k9/JOURNAL.md index 1eecb84d..9655906f 100644 --- a/projects/keystone-fea-b3k9/JOURNAL.md +++ b/projects/keystone-fea-b3k9/JOURNAL.md @@ -77,11 +77,25 @@ implemented from scratch, numerically validated in-app. - `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, +- `src/engine/buckling.ts` — **continuum linear (Euler) buckling** (v14): assembles the geometric + (initial-stress) stiffness K_g = ∫GᵀτG from a reference plane-stress field on the Q4/Q8 + isoparametric machinery, and solves the generalized eigenproblem (K + λK_g)φ = 0 for the lowest + critical load factors + buckled mode shapes by a **K-orthonormal subspace iteration** on + A = K⁻¹(−K_g) — a standard-symmetric Rayleigh–Ritz reduction with relative rank-revealing deflation + that never trips on the indefinite stress stiffness. Plus `eulerLoad` = π²EI/(KL)². Pure, testable. +- `src/engine/bucklingpresets.ts` — the buckling scenario library (fixed–free flagpole K=2, + fixed–pinned braced strut K≈0.699, compression panel) + a material table with yield stress, and the + builder that returns section properties (I, r, Euler load, slenderness KL/r) for the analytical curve. +- `src/ui/BucklingStudio.tsx` — the Buckling & Stability tab: own canvas + rAF loop, the buckled shape + animated as a breathing lateral bow over the axial load-path shading, and the Euler **column curve** + (σ_cr = π²E/λ² with the yield cut-off + the FE design point), a mode-pill selector and a Q8/Q4 + shear-locking toggle. +- `src/engine/validate.ts` — analytical benchmarks that run live (84 of them): truss statics, cantilever PL³/3EI, 5wL⁴/384EI, patch test, modal/buckling, harmonic/FRF, plastic collapse, 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). + test, Q8 Euler & Timoshenko cantilever, Q4 refinement, continuum bending frequency), the + v13 fracture checks (Feddersen/Tada geometry factors, J = K²/E*, Griffith, J path independence), + and the v14 continuum-buckling checks (fixed–free & fixed–pinned Euler columns, the 1:9 mode ladder). 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. @@ -659,8 +673,125 @@ ceramic and glass. - [ ] **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. +### v14 — Continuum linear buckling & stability (the *column* chapter) + +The frame solver has answered "when does it go unstable?" since v2 — Euler buckling of beam-columns +via a geometric stiffness built from the axial force in each member. But the **continuum** — the Q4/Q8 +plane-stress world where v9–v13 live — could only ever answer "how much does it sag?" A slender 2-D +part under compression has a second, non-negotiable failure mode the static solver is blind to: it +**buckles**, snapping sideways at a load *far* below the stress that would yield it. That is a +stiffness collapse, not a strength one, and it is exactly the question a plate, a web, a wall or a +strut modeled as a continuum has to answer. + +The physics is a generalized eigenproblem. Under a reference in-plane stress field σ₀ (from a unit +reference load) each element gains a **geometric (initial-stress) stiffness** + + K_g = ∫ Gᵀ τ G dΩ , τ = [[σxx, σxy], [σxy, σyy]] (block-diagonal in the u/v gradients) + +where G gathers the displacement gradients ∇u, ∇v. The structure loses stability when the tangent +stiffness first goes singular: + + (K + λ K_g) φ = 0 ⇔ K φ = λ (−K_g) φ , + +i.e. the smallest positive load multiplier λ_cr that admits a non-trivial mode φ. λ_cr × (reference +load) is the buckling load; φ is the buckled shape. Everything reuses the studio's own isoparametric +kinematics (`isoparam`), sparse assembler + BC-aware PCG (`linalg`) and the dense reduced eigensolver +(`eigen`) — no new linear-algebra primitives. + +The eigenproblem is indefinite (K_g is not sign-definite) so the v9 modal subspace machinery, which +assumes an SPD mass, cannot be reused verbatim. Instead the solver runs a **K-orthonormal subspace +iteration** on the operator A = K⁻¹(−K_g): it is self-adjoint in the K-inner-product ⟨a,b⟩_K = aᵀK b, +so each pass solves a handful of sparse systems, K-orthonormalises the iterate block with a +**rank-revealing (relative) deflation** that drops the rigid / zero-stress directions, and reads the +dominant Ritz pairs from a *standard* symmetric reduced problem Hᵣ = Qᵀ(−K_g)Q — which, unlike a +Cholesky-based generalized reduction, never trips on the indefinite stress stiffness. The dominant +eigenvalues of A are the smallest |λ|: exactly the modes that buckle first. + +- [x] `buckling.ts` — the whole chapter, pure/deterministic. Assemble K, run the reference static + solve for the pre-buckling stress field, recover Gauss-point stresses, assemble the **stress + stiffness** S = −K_g (block-diagonal identical u/v blocks, from the B-matrix gradients), and solve + (K + λK_g)φ = 0 by the K-orthonormal subspace iteration above +- [x] Robustness: the reduced problem is a **standard** symmetric eigenproblem (`jacobiEig`), not a + Cholesky generalized one, so the indefinite/rank-deficient S never breaks it; deflation is + **relative** (residual K-norm ÷ original), because in SI units the working K-norms are ~1e-11 and + an absolute floor wrongly deflates everything (the bug that made the first cut return no modes) +- [x] `eulerLoad(E,I,L,K)` = π²EI/(KL)² — the closed-form reference the studio and benchmarks compare to +- [x] `bucklingpresets.ts` — a material table (steel / Al 6061 / Ti-6Al-4V, each with σ_y), three + scenarios (fixed–free **flagpole** K=2, fixed–pinned **braced strut** K≈0.699, a **compression + panel** with no 1-D analogue) and a builder that also returns the section properties (I, r, Euler + load, slenderness KL/r) so the studio can draw the analytical column +- [x] `BucklingStudio.tsx` — a sixth self-contained studio tab (own canvas + rAF loop, snapshot-ref + discipline). Two views share one solve: the **buckled shape** (the pre-buckling straight part + shaded by its axial load path, the eigenmode animated as a breathing lateral bow, base hatching + + load arrows) and the **column curve** (the FE critical stress dropped onto the Euler hyperbola + σ_cr=π²E/λ² with its yield cut-off and transition slenderness — read buckling-vs-squash at a + glance). Mode-pill selector, a Q8/Q4 toggle that *shows shear locking*, a debounced slenderness + slider, and a full FE-vs-Euler read-out +- [x] **3 new closed-form benchmarks** (badge **81 → 84**, all green): the **fixed–free Euler column** + P_cr = π²EI/(2L)² to **0.05 %**, the **fixed–pinned** column P_cr = π²EI/(0.699L)² to **0.5 %**, + and the **mode ladder** λ₂/λ₁ = 3²/1² = 9 to 0.25 % — the odd-quarter-wave signature of a + cantilever column +- [x] Q8 vs Q4: the benchmark uses Q8 deliberately — fully-integrated Q4 **shear-locks** in bending and + overstiffens the buckling load by ~16 % at the same mesh, so the studio defaults to Q8 and offers + Q4 as a live "watch it lock" toggle +- [x] Performance: the panel's clustered eigenvalues made the first cut take 56 s; capped subspace + iterations, loosened the inner-CG tolerance (the Ritz step recovers the eigenvalue accuracy) and + coarsened the meshes to keep every solve interactive, with the slenderness slider debounced so a + drag re-solves once, on release +- [x] Verified end-to-end in headless Chromium: all six studio tabs render; the flagpole bows into its + first mode with FE-vs-Euler **0.08 %** and the 1 : 8.97 : 24.79 ladder, the braced strut shows the + fixed–pinned half-wave at 0.63 %, the panel buckles into a 2-D bulge flagged as *squash*, the + column curve seats the design point on the Euler hyperbola, and the badge reads **84/84** with the + whole gate (scope + conformance + lint + build) green and zero runtime errors + +**v15 backlog — where the stability chapter goes next:** + +- [ ] **True pinned–pinned continuum column** — a clean Euler K=1 case needs the base free to *rotate* + (only the centre point pinned in y), which makes the reference stress field non-uniform; resolve it + with a St-Venant load-introduction block or a prescribed uniform prestress so the classic 1:4:9 + ladder validates too. +- [ ] **Nonlinear (geometric) buckling** — replace the eigenvalue estimate with an arc-length / + Newton path that follows the load–deflection curve *through* the bifurcation, capturing + post-buckling stiffening and imperfection sensitivity (the knock-down real shells suffer). +- [ ] **Classic plate-buckling coefficients** — a properly simply-supported rectangular panel under + uniaxial / biaxial / shear load, validated against the handbook k in σ_cr = k·π²E/(12(1−ν²))·(t/b)² + and its half-wave count vs aspect ratio. +- [ ] **Lateral–torsional & local buckling of the real AISC sections** (marry v4's section library to + the geometric stiffness) so a wide-flange beam's LTB capacity falls out. +- [ ] **Buckling-constrained topology optimization** — feed λ_cr back as a constraint in the v10 SIMP + loop so the generated layout is stiff *and* stable, not just compliant-optimal. + ## Session log +- 2026-07-24 (claude): shipped **v14 — continuum linear buckling & stability**, the chapter that gives + the Q4/Q8 continuum the second failure mode the static solver is blind to: elastic instability. New + `buckling.ts` assembles the **geometric (initial-stress) stiffness** K_g = ∫GᵀτG from the + pre-buckling stress field (a reference static solve on the studio's own isoparametric machinery), + recovers the Gauss-point stresses, and solves the generalized eigenproblem (K + λK_g)φ = 0 for the + lowest **critical load factors** and their buckled mode shapes. Because K_g is indefinite, the v9 + SPD-mass subspace solver can't be reused; instead a **K-orthonormal subspace iteration** on the + operator A = K⁻¹(−K_g) — self-adjoint in the K-inner-product — projects onto a *standard* symmetric + reduced problem Hᵣ = Qᵀ(−K_g)Q solved by cyclic Jacobi, with **relative** rank-revealing deflation of + the rigid directions (an absolute deflation floor was the first-cut bug: SI-unit K-norms sit at ~1e-11 + and everything got wrongly dropped, returning no modes). `bucklingpresets.ts` adds a material table + (steel / Al 6061 / Ti-6Al-4V, with σ_y), three scenarios — fixed–free **flagpole** (K=2), fixed–pinned + **braced strut** (K≈0.699) and a **compression panel** — and returns the section properties so the + studio can draw the analytical column. `BucklingStudio.tsx` is a sixth studio tab: a **buckled-shape** + view (the straight part shaded by its axial load path, the eigenmode animated as a breathing lateral + bow, base hatching + load arrows) and a **column-curve** view that drops the FE critical stress onto + the Euler hyperbola σ_cr = π²E/λ² with its yield cut-off, so buckling-vs-squash reads at a glance — + plus a mode-pill selector, a debounced slenderness slider, and a **Q8/Q4 toggle that visibly shows + shear locking** (fully-integrated Q4 overstiffens the buckling load ~16 %, so Q8 is the default). + Three new closed-form benchmarks pin it (badge **81 → 84**, all green): the fixed–free Euler column + P_cr = π²EI/(2L)² to **0.05 %**, the fixed–pinned column (K≈0.699) to 0.5 %, and the cantilever mode + ladder λ₂/λ₁ = 9 (3²/1²) to 0.25 %. Verified end-to-end in headless Chromium — all six tabs render, + the flagpole reproduces Euler to 0.08 % and the 1:8.97:24.79 ladder, the braced strut shows the + fixed–pinned half-wave, the panel buckles into a 2-D bulge, the column curve seats the design point on + the hyperbola — and the whole gate (scope + conformance + lint + build) is green with **84/84** + benchmarks and zero runtime errors. First cut had the panel taking 56 s (clustered eigenvalues) and + the eigensolver returning nothing (absolute-deflation bug); both fixed — every solve is now + interactive and the math is validated to a fraction of a percent against Euler. + - 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 diff --git a/projects/keystone-fea-b3k9/project.json b/projects/keystone-fea-b3k9/project.json index f3e6b997..44f6450b 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. 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.", + "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. And the newest chapter carries STABILITY into the continuum: a CONTINUUM LINEAR-BUCKLING studio assembles the GEOMETRIC (initial-stress) STIFFNESS K_g=∫Gᵀτ G from a reference plane-stress field and solves the generalized eigenproblem (K + λ K_g) φ = 0 by K-orthonormal SUBSPACE ITERATION (sparse CG inner solves with a Rayleigh–Ritz reduction that never trips on the indefinite stress stiffness) for the lowest critical load factors and their buckled mode shapes — a flagpole (fixed–free, K=2), a braced strut (fixed–pinned, K≈0.699) and a compression panel, each animated as a breathing lateral bow and dropped onto the live EULER COLUMN CURVE (σ_cr=π²E/λ² with its yield cut-off) so you read buckling-vs-squash at a glance; the FE critical load reproduces Euler's π²EI/(KL)² to 0.05% and the fixed–free mode ladder to its exact 1 : 9 : 25 (odd² quarter-waves). 84 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", "fracture", "fracture-mechanics", "lefm", "stress-intensity", "j-integral", "interaction-integral", "crack", "williams-field", "griffith", "toughness", "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", "continuum-buckling", "stability", "geometric-stiffness", "euler", "column-curve", "eigenvalue-buckling", "react"], "createdAt": "2026-07-10" } diff --git a/projects/keystone-fea-b3k9/src/App.css b/projects/keystone-fea-b3k9/src/App.css index 6b892b32..74f27784 100644 --- a/projects/keystone-fea-b3k9/src/App.css +++ b/projects/keystone-fea-b3k9/src/App.css @@ -833,3 +833,21 @@ body { border-color: color-mix(in srgb, #ff6b81 50%, var(--line-2)); background: color-mix(in srgb, #ff6b81 12%, var(--panel-2)); } + +/* ---------------------------------------------------- buckling studio extras */ +.mode-row .mode-val { + text-align: right; + font-weight: 600; +} +.mode-row .mode-ratio { + text-align: right; + color: var(--faint); +} +.mode-pills { + display: flex; + gap: 4px; + align-items: center; +} +.seg-row { + margin: 4px 0 2px; +} diff --git a/projects/keystone-fea-b3k9/src/App.tsx b/projects/keystone-fea-b3k9/src/App.tsx index 55da8d7f..5b8e1e80 100644 --- a/projects/keystone-fea-b3k9/src/App.tsx +++ b/projects/keystone-fea-b3k9/src/App.tsx @@ -56,6 +56,7 @@ import { fmtEng } from './ui/format' import { TopOptStudio } from './ui/TopOptStudio' import { ThermalStudio } from './ui/ThermalStudio' import { FractureStudio } from './ui/FractureStudio' +import { BucklingStudio } from './ui/BucklingStudio' import { addMember, addNode, @@ -82,7 +83,7 @@ import { type Scene, } from './state' -type Tab = 'frame' | 'continuum' | 'topopt' | 'thermal' | 'fracture' +type Tab = 'frame' | 'continuum' | 'topopt' | 'thermal' | 'fracture' | 'buckling' type Tool = 'select' | 'node' | 'member' | 'support' | 'load' | 'delete' const TOOLS: { id: Tool; label: string; hint: string }[] = [ @@ -966,6 +967,7 @@ export default function App() { { value: 'topopt', label: 'Topology Optimization' }, { value: 'thermal', label: 'Thermal & Multiphysics' }, { value: 'fracture', label: 'Fracture Mechanics' }, + { value: 'buckling', label: 'Buckling & Stability' }, ]} value={tab} onChange={(v) => { @@ -984,6 +986,8 @@ export default function App() { ) : tab === 'fracture' ? ( + ) : tab === 'buckling' ? ( + ) : ( <> {/* ---------------- left rail: presets ---------------- */} diff --git a/projects/keystone-fea-b3k9/src/engine/buckling.ts b/projects/keystone-fea-b3k9/src/engine/buckling.ts new file mode 100644 index 00000000..69a2a52d --- /dev/null +++ b/projects/keystone-fea-b3k9/src/engine/buckling.ts @@ -0,0 +1,647 @@ +// Continuum linear (Euler) buckling — the stability chapter. +// +// The static solver answers "how much does it sag?" Buckling answers a +// different, non-negotiable question: "at what load does it stop being stable +// and snap sideways?" A slender column fails this way *far* below its yield +// stress — a stiffness collapse, not a strength one — and no amount of static +// stress checking sees it coming. +// +// The physics is a generalized eigenproblem. Under a reference in-plane stress +// field σ₀ (from a unit reference load) an element gains a *geometric* (a.k.a. +// initial-stress) stiffness K_g = ∫ Gᵀ τ G dΩ, where G gathers the displacement +// gradients and τ is the block-diagonal stress matrix. The structure loses +// stability when +// +// (K + λ K_g) φ = 0 ⇔ K φ = λ (−K_g) φ , +// +// i.e. the smallest positive load multiplier λ that makes the tangent stiffness +// singular. λ_cr · (reference load) is the buckling (critical) load and φ is the +// buckling mode shape. +// +// We reuse the isoparametric Q4/Q8 kinematics (`isoparam.ts`), the sparse CSR +// assembler + BC-aware PCG (`linalg.ts`) for K and the reference solve, and the +// dense reduced eigensolver (`eigen.ts`). The eigenproblem itself is solved by a +// Bathe-style **subspace iteration** driven by the operator A = K⁻¹(−K_g): each +// pass solves a handful of sparse systems, projects onto the iterate subspace, +// and reads the dominant Ritz pairs — the *lowest* buckling factors — from a +// tiny reduced generalized problem whose "mass" matrix Kᵣ = X̄ᵀK X̄ is guaranteed +// SPD, so the reduction never breaks on the indefinite stress stiffness. Pure, +// testable. + +import { Assembler, matVec, solveCG, dot, type CSR } from './linalg' +import { jacobiEig } from './eigen' +import { + elementKinematics, + extrapMatrix, + planeStressD, + type QOrder, +} from './isoparam' +import { + edgeNodesQ, + boundaryElementEdges, + type EdgeName, + type QuadMesh, +} from './quadmesh' + +export interface FixGroupB { + edge?: EdgeName + nodes?: number[] + dofs: ('x' | 'y')[] +} + +export interface BucklingInput { + mesh: QuadMesh + E: number + nu: number + thickness: number + fix: FixGroupB[] + /** Reference (unit) load whose λ multiplier the eigenproblem reports. */ + traction?: { edge: EdgeName; tx: number; ty: number } + pointLoads?: { node: number; fx: number; fy: number }[] + bodyForce?: { gx: number; gy: number } + /** How many buckling modes to extract (default 4). */ + nModes?: number +} + +export interface BucklingMode { + /** Critical load multiplier λ (× reference load = buckling load). */ + loadFactor: number + dispX: Float64Array + dispY: Float64Array + /** Peak nodal amplitude (for auto-scaling the drawn mode). */ + maxAmp: number +} + +export interface BucklingResult { + order: QOrder + modes: BucklingMode[] + freeDofCount: number + /** Pre-buckling (reference) displacement field. */ + refDispX: Float64Array + refDispY: Float64Array + refMaxDisp: number + /** Smooth nodal reference stress field (for the load-path shading). */ + refSxx: Float64Array + refSyy: Float64Array + refSxy: Float64Array + refSvm: Float64Array + /** Mean axial (σ_yy) reference stress over the mesh — a scalar handle. */ + refMeanSyy: number + stable: boolean + eigIterations: number + eigConverged: boolean +} + +const STEEL_RHO = 7850 + +/** Free-DOF mask: fixed groups + orphan (element-less) nodes clamped. */ +function freeMask(mesh: QuadMesh, fix: FixGroupB[], touched: Uint8Array): Uint8Array { + const nDof = mesh.nodeCount * 2 + const free = new Uint8Array(nDof).fill(1) + for (const g of fix) { + const nodes = g.edge ? edgeNodesQ(mesh, g.edge) : (g.nodes ?? []) + for (const n of nodes) for (const d of g.dofs) free[n * 2 + (d === 'x' ? 0 : 1)] = 0 + } + for (let n = 0; n < mesh.nodeCount; n++) + if (!touched[n]) { + free[n * 2] = 0 + free[n * 2 + 1] = 0 + } + return free +} + +/** Reference load vector: consistent edge traction + point/body loads. */ +function buildLoad(input: BucklingInput): Float64Array { + const { mesh, thickness } = input + const order = mesh.order + const nne = order + const f = new Float64Array(mesh.nodeCount * 2) + if (input.bodyForce) { + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const xs: number[] = [] + const ys: number[] = [] + for (let a = 0; a < nne; a++) { + const n = mesh.elems[base + a] + xs.push(mesh.x[n]) + ys.push(mesh.y[n]) + } + for (const g of elementKinematics(order, xs, ys)) { + const w = g.w * Math.abs(g.detJ) * thickness + for (let a = 0; a < nne; a++) { + const n = mesh.elems[base + a] + f[n * 2] += input.bodyForce.gx * g.N[a] * w + f[n * 2 + 1] += input.bodyForce.gy * g.N[a] * w + } + } + } + } + if (input.traction) { + const { edge, tx, ty } = input.traction + for (const nodesOnEdge of boundaryElementEdges(mesh, edge)) { + const a = nodesOnEdge[0] + const b = nodesOnEdge[nodesOnEdge.length - 1] + const len = Math.hypot(mesh.x[b] - mesh.x[a], mesh.y[b] - mesh.y[a]) + const w = len * thickness + if (nodesOnEdge.length === 2) { + for (const n of nodesOnEdge) { + f[n * 2] += (tx * w) / 2 + f[n * 2 + 1] += (ty * w) / 2 + } + } else { + const frac = [1 / 6, 2 / 3, 1 / 6] + for (let k = 0; k < 3; k++) { + const n = nodesOnEdge[k] + f[n * 2] += tx * w * frac[k] + f[n * 2 + 1] += ty * w * frac[k] + } + } + } + } + for (const pl of input.pointLoads ?? []) { + f[pl.node * 2] += pl.fx + f[pl.node * 2 + 1] += pl.fy + } + return f +} + +/** + * Assemble the elastic stiffness K and cache the per-element Gauss kinematics + * (so the reference-stress and geometric-stiffness passes never recompute the + * Jacobian / B matrices). + */ +function assembleElastic(input: BucklingInput) { + const { mesh, E, nu, thickness } = input + const order = mesh.order + const nne = order + const nDof = mesh.nodeCount * 2 + const D = planeStressD(E, nu) + const asmK = new Assembler(nDof) + const touched = new Uint8Array(mesh.nodeCount) + const elemKin: ReturnType[] = [] + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const idx: number[] = [] + const xs: number[] = [] + const ys: 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 kin = elementKinematics(order, xs, ys) + elemKin.push(kin) + const ndofE = 2 * nne + const Ke = Array.from({ length: ndofE }, () => new Array(ndofE).fill(0)) + for (const g of kin) { + const scale = g.w * Math.abs(g.detJ) * 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 m = 0; m < 3; m++) s += D[r][m] * g.B[m][c] + DB[r][c] = s + } + for (let a = 0; a < ndofE; a++) + for (let b = 0; b < ndofE; b++) { + let s = 0 + for (let r = 0; r < 3; r++) s += g.B[r][a] * DB[r][b] + Ke[a][b] += s * scale + } + } + const dofs: number[] = [] + for (const n of idx) dofs.push(n * 2, n * 2 + 1) + asmK.addBlock(dofs, Ke) + } + return { K: asmK.build(), D, touched, elemKin } +} + +/** Gauss-point reference stresses [σxx, σyy, σxy] per element, from u_ref. */ +function gaussStresses( + mesh: QuadMesh, + D: number[][], + elemKin: ReturnType[], + u: Float64Array, +): number[][][] { + const order = mesh.order + const nne = order + const out: number[][][] = [] + 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 perGP: number[][] = [] + for (const g of elemKin[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] + perGP.push([sxx, syy, sxy]) + } + out.push(perGP) + } + return out +} + +/** + * Assemble the **stress stiffness** S = −K_g from the reference Gauss stresses, + * so the eigenproblem reads K φ = λ S φ with the critical multiplier λ > 0 in + * compression. For plane elasticity S is block-diagonal in the x/y DOFs with an + * identical scalar block, because the stress matrix τ couples only ∇u with ∇u + * and ∇v with ∇v: + * + * S_e[2i][2j] = S_e[2i+1][2j+1] = + * −∫ (Nᵢ,ₓ(σxx Nⱼ,ₓ + σxy Nⱼ,ᵧ) + Nᵢ,ᵧ(σxy Nⱼ,ₓ + σyy Nⱼ,ᵧ)) t dΩ + */ +function assembleStressStiffness( + input: BucklingInput, + elemKin: ReturnType[], + gpS: number[][][], +): CSR { + const { mesh, thickness } = input + const order = mesh.order + const nne = order + const nDof = mesh.nodeCount * 2 + const asm = new Assembler(nDof) + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const kin = elemKin[e] + const ndofE = 2 * nne + const Se = Array.from({ length: ndofE }, () => new Array(ndofE).fill(0)) + for (let gi = 0; gi < kin.length; gi++) { + const g = kin[gi] + const [sxx, syy, sxy] = gpS[e][gi] + const scale = g.w * Math.abs(g.detJ) * thickness + // Recover the physical shape-function gradients from the B matrix: + // dNx_i = B[0][2i] (∂u/∂x row), dNy_i = B[1][2i+1] (∂v/∂y row). + const dNx: number[] = new Array(nne) + const dNy: number[] = new Array(nne) + for (let i = 0; i < nne; i++) { + dNx[i] = g.B[0][2 * i] + dNy[i] = g.B[1][2 * i + 1] + } + for (let i = 0; i < nne; i++) + for (let j = 0; j < nne; j++) { + const t = + dNx[i] * (sxx * dNx[j] + sxy * dNy[j]) + + dNy[i] * (sxy * dNx[j] + syy * dNy[j]) + const v = -t * scale // S = −K_g + Se[2 * i][2 * j] += v + Se[2 * i + 1][2 * j + 1] += v + } + } + const dofs: number[] = [] + for (let a = 0; a < nne; a++) { + const n = mesh.elems[base + a] + dofs.push(n * 2, n * 2 + 1) + } + asm.addBlock(dofs, Se) + } + return asm.build() +} + +/** + * Lowest `p` buckling factors of K φ = λ S φ (S = −K_g, symmetric indefinite), + * by subspace iteration on the operator A = K⁻¹ S. The operator A is + * self-adjoint in the K-inner-product ⟨a,b⟩_K = aᵀK b, so we work entirely in + * that metric. Each pass: + * + * 1. Z = A X = K⁻¹ (S X) — q sparse PCG solves, + * 2. K-orthonormalise Z → Q — modified Gram–Schmidt in ⟨·,·⟩_K, with + * rank-revealing deflation of null-stress + * (rigid) directions, + * 3. reduced Rayleigh–Ritz operator Hᵣ = Qᵀ S Q (r×r symmetric; equals + * Qᵀ K A Q because Q is K-orthonormal), solved by dense Jacobi, + * 4. rotate X ← Q V (Ritz vectors), refill to q with fresh seeds. + * + * The subspace converges to the dominant eigenvectors of A (largest |μ|), whose + * μ = 1/λ are the *smallest* |λ| — exactly the modes that buckle first. Because + * the reduced problem is a *standard* symmetric eigenproblem, it never trips on + * the indefinite / rank-deficient stress stiffness. We keep the positive λ + * (compression instability) and sort ascending. + */ +function subspaceBuckling( + K: CSR, + S: CSR, + free: Uint8Array, + p: number, +): { factors: number[]; vectors: Float64Array[]; iterations: number; converged: boolean } { + const n = K.n + let nf = 0 + for (let i = 0; i < n; i++) if (free[i]) nf++ + const q = Math.min(nf, p + 4) + if (q <= 0) return { factors: [], vectors: [], iterations: 0, converged: false } + + const freeIdx: number[] = [] + for (let i = 0; i < n; i++) if (free[i]) freeIdx.push(i) + + const maskFree = (y: Float64Array) => { + for (let i = 0; i < n; i++) if (!free[i]) y[i] = 0 + } + const mul = (A: CSR, x: Float64Array) => { + const y = new Float64Array(n) + matVec(A, x, y) + maskFree(y) + return y + } + const seed = (c: number) => { + const v = new Float64Array(n) + for (let k = 0; k < freeIdx.length; k++) + v[freeIdx[k]] = Math.sin((k + 1) * 0.6180339 + (c + 1) * 1.31459) + return v + } + + /** + * K-orthonormalise a block of vectors (modified Gram–Schmidt in the + * K-inner-product). Returns the accepted basis Q and the cached K·qⱼ. Vectors + * whose residual K-norm falls below a relative floor are deflated (dropped) — + * these are the rigid / zero-stress directions the stress stiffness ignores. + */ + const kOrthonormalise = (block: Float64Array[]) => { + const Q: Float64Array[] = [] + const KQ: Float64Array[] = [] + for (const raw of block) { + const v = raw.slice() + maskFree(v) + const Kv = mul(K, v) + const nrm0 = Math.sqrt(Math.max(dot(v, Kv), 0)) // K-norm before projecting out + for (let j = 0; j < Q.length; j++) { + const coef = dot(v, KQ[j]) + if (coef === 0) continue + for (let i = 0; i < n; i++) { + v[i] -= coef * Q[j][i] + Kv[i] -= coef * KQ[j][i] + } + } + const nrm = Math.sqrt(Math.max(dot(v, Kv), 0)) + // Rank-revealing deflation is *relative*: a direction that all but + // vanishes against the already-accepted basis (or was numerically zero to + // begin with) is dropped. Absolute floors fail because the working K-norms + // here are ~1e-11 in SI units. + if (nrm0 <= 0 || nrm / nrm0 < 1e-8) continue // deflate + const inv = 1 / nrm + for (let i = 0; i < n; i++) { + v[i] *= inv + Kv[i] *= inv + } + Q.push(v) + KQ.push(Kv) + } + return Q + } + + const dominant = (mu: number[]) => + mu + .slice() + .sort((a, b) => Math.abs(b) - Math.abs(a)) + .slice(0, p) + + // Initial K-orthonormal subspace from deterministic seeds. + let X = kOrthonormalise(Array.from({ length: q }, (_, c) => seed(c))) + + let prev: number[] = [] + let iterations = 0 + let converged = false + const maxIter = 40 + let lastMu: number[] = [] + let lastQ: Float64Array[] = X + let lastV: number[][] = [] + + for (let iter = 0; iter < maxIter; iter++) { + iterations = iter + 1 + // 1. Z = A X = K⁻¹ (S X). CG tolerance is loose here — the subspace only + // needs eigenvector *directions*, and the reduced Rayleigh–Ritz step + // recovers the eigenvalues to full accuracy from those. + const Z: Float64Array[] = X.map((x) => { + const rhs = mul(S, x) + const sol = solveCG(K, rhs, free, { tol: 1e-8, maxIter: Math.max(1500, 20 * n) }) + return sol.x + }) + // 2. K-orthonormalise → Q + const Q = kOrthonormalise(Z) + if (Q.length === 0) break + const r = Q.length + // 3. reduced Rayleigh–Ritz: Hᵣ = Qᵀ S Q + const SQ = Q.map((qv) => mul(S, qv)) + const H = Array.from({ length: r }, () => new Array(r).fill(0)) + for (let a = 0; a < r; a++) + for (let b = a; b < r; b++) { + const v = dot(Q[a], SQ[b]) + H[a][b] = H[b][a] = v + } + const eig = jacobiEig(H) + lastMu = eig.values + lastQ = Q + lastV = eig.vectors + // 4. new block = Ritz vectors (Q·V), refilled to q with fresh seeds. + const ritz: Float64Array[] = [] + for (let k = 0; k < r; k++) { + const v = new Float64Array(n) + for (let l = 0; l < r; l++) { + const c = eig.vectors[l][k] + const qv = Q[l] + for (let i = 0; i < n; i++) v[i] += c * qv[i] + } + ritz.push(v) + } + const refill = [...ritz] + for (let c = 0; refill.length < q; c++) refill.push(seed(q + iter * 7 + c)) + X = kOrthonormalise(refill) + + const cur = dominant(lastMu) + if (prev.length === cur.length && cur.length > 0) { + let maxRel = 0 + for (let i = 0; i < cur.length; i++) + maxRel = Math.max(maxRel, Math.abs(cur[i] - prev[i]) / (Math.abs(cur[i]) + 1e-300)) + if (maxRel < 1e-7) { + converged = true + break + } + } + prev = cur + } + + // Assemble final Ritz pairs from the last reduced solve. + const pairs: { lambda: number; vec: Float64Array }[] = [] + const r = lastQ.length + for (let k = 0; k < lastMu.length; k++) { + const mu = lastMu[k] + if (!Number.isFinite(mu) || Math.abs(mu) < 1e-300) continue + const lambda = 1 / mu + if (!Number.isFinite(lambda) || lambda <= 0) continue + const vec = new Float64Array(n) + for (let l = 0; l < r; l++) { + const c = lastV[l]?.[k] ?? 0 + const qv = lastQ[l] + for (let i = 0; i < n; i++) vec[i] += c * qv[i] + } + maskFree(vec) + pairs.push({ lambda, vec }) + } + pairs.sort((a, b) => a.lambda - b.lambda) + const top = pairs.slice(0, p) + return { + factors: top.map((t) => t.lambda), + vectors: top.map((t) => t.vec), + iterations, + converged, + } +} + +/** Smooth nodal stress recovery (Gauss→node extrapolation + averaging). */ +function recoverNodalStress( + mesh: QuadMesh, + gpS: number[][][], +): { sxx: Float64Array; syy: Float64Array; sxy: Float64Array; svm: Float64Array; meanSyy: number } { + const order = mesh.order + const nne = order + const Emat = extrapMatrix(order) + const accSxx = new Float64Array(mesh.nodeCount) + const accSyy = new Float64Array(mesh.nodeCount) + const accSxy = new Float64Array(mesh.nodeCount) + const count = new Float64Array(mesh.nodeCount) + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + const perGP = gpS[e] + for (let a = 0; a < nne; a++) { + let sxx = 0 + let syy = 0 + let sxy = 0 + for (let gi = 0; gi < perGP.length; gi++) { + sxx += Emat[a][gi] * perGP[gi][0] + syy += Emat[a][gi] * perGP[gi][1] + sxy += Emat[a][gi] * perGP[gi][2] + } + const n = mesh.elems[base + a] + accSxx[n] += sxx + accSyy[n] += syy + accSxy[n] += sxy + count[n] += 1 + } + } + const sxx = new Float64Array(mesh.nodeCount) + const syy = new Float64Array(mesh.nodeCount) + const sxy = new Float64Array(mesh.nodeCount) + const svm = new Float64Array(mesh.nodeCount) + let sumSyy = 0 + let nn = 0 + for (let n = 0; n < mesh.nodeCount; n++) { + if (count[n] === 0) continue + const a = accSxx[n] / count[n] + const b = accSyy[n] / count[n] + const c = accSxy[n] / count[n] + sxx[n] = a + syy[n] = b + sxy[n] = c + svm[n] = Math.sqrt(a * a - a * b + b * b + 3 * c * c) + sumSyy += b + nn++ + } + return { sxx, syy, sxy, svm, meanSyy: nn ? sumSyy / nn : 0 } +} + +/** + * Full linear-buckling pipeline: assemble K, solve the reference static problem + * for the pre-buckling stress field, build the stress stiffness S = −K_g, and + * solve the generalized eigenproblem for the lowest critical load factors and + * their mode shapes. + */ +export function solveBuckling(input: BucklingInput): BucklingResult { + const order = input.mesh.order + const nModes = input.nModes ?? 4 + const { mesh } = input + const nDof = mesh.nodeCount * 2 + + const { K, D, touched, elemKin } = assembleElastic(input) + const free = freeMask(mesh, input.fix, touched) + let freeDofCount = 0 + for (let i = 0; i < nDof; i++) if (free[i]) freeDofCount++ + + // Reference static solve → pre-buckling stress field. + const f = buildLoad(input) + const sol = solveCG(K, f, free, { tol: 1e-11, maxIter: Math.max(2000, 60 * nDof) }) + const u = sol.x + + const refDispX = new Float64Array(mesh.nodeCount) + const refDispY = new Float64Array(mesh.nodeCount) + let refMaxDisp = 0 + for (let n = 0; n < mesh.nodeCount; n++) { + refDispX[n] = u[n * 2] + refDispY[n] = u[n * 2 + 1] + refMaxDisp = Math.max(refMaxDisp, Math.hypot(refDispX[n], refDispY[n])) + } + + const gpS = gaussStresses(mesh, D, elemKin, u) + const nodal = recoverNodalStress(mesh, gpS) + const S = assembleStressStiffness(input, elemKin, gpS) + + const eig = subspaceBuckling(K, S, free, nModes) + + const modes: BucklingMode[] = eig.factors.map((loadFactor, k) => { + const vec = eig.vectors[k] + const dispX = new Float64Array(mesh.nodeCount) + const dispY = new Float64Array(mesh.nodeCount) + let maxAmp = 0 + for (let n = 0; n < mesh.nodeCount; n++) { + dispX[n] = vec[n * 2] + dispY[n] = vec[n * 2 + 1] + maxAmp = Math.max(maxAmp, Math.hypot(dispX[n], dispY[n])) + } + // Normalise the mode shape to unit peak amplitude for stable drawing. + if (maxAmp > 0) + for (let n = 0; n < mesh.nodeCount; n++) { + dispX[n] /= maxAmp + dispY[n] /= maxAmp + } + return { loadFactor, dispX, dispY, maxAmp: 1 } + }) + + const stable = + sol.converged && + Number.isFinite(refMaxDisp) && + modes.length > 0 && + Number.isFinite(modes[0].loadFactor) + + return { + order, + modes, + freeDofCount, + refDispX, + refDispY, + refMaxDisp, + refSxx: nodal.sxx, + refSyy: nodal.syy, + refSxy: nodal.sxy, + refSvm: nodal.svm, + refMeanSyy: nodal.meanSyy, + stable, + eigIterations: eig.iterations, + eigConverged: eig.converged, + } +} + +/** + * Euler critical load for an ideal prismatic column: + * P_cr = π² E I / (K_eff L)² + * K_eff is the effective-length factor for the end conditions (2.0 cantilever + * fixed-free, 1.0 pinned-pinned, 0.5 fixed-fixed, 0.699… fixed-pinned). + */ +export function eulerLoad(E: number, I: number, L: number, kEff: number): number { + return (Math.PI * Math.PI * E * I) / (kEff * L * (kEff * L)) +} + +export { STEEL_RHO } diff --git a/projects/keystone-fea-b3k9/src/engine/bucklingpresets.ts b/projects/keystone-fea-b3k9/src/engine/bucklingpresets.ts new file mode 100644 index 00000000..4d5cf978 --- /dev/null +++ b/projects/keystone-fea-b3k9/src/engine/bucklingpresets.ts @@ -0,0 +1,198 @@ +// Buckling scenario library: prismatic columns with a range of textbook end +// conditions plus a compression panel, each expressed as a mesh + boundary / +// load recipe for the continuum buckling solver, together with the closed-form +// Euler reference (effective-length factor, second moment, radius of gyration) +// so the studio can plot the FE critical load against the analytical column. + +import { rectPlateQ, type QuadMesh } from './quadmesh' +import { eulerLoad, type BucklingInput } from './buckling' +import type { QOrder } from './isoparam' + +export interface BMaterial { + id: string + name: string + E: number // Pa + nu: number + rho: number // kg/m³ + sigmaY: number // Pa (yield / squash stress) +} + +export const BUCKLING_MATERIALS: BMaterial[] = [ + { id: 'steel', name: 'Structural steel', E: 210e9, nu: 0.3, rho: 7850, sigmaY: 250e6 }, + { id: 'aluminium', name: 'Aluminium 6061', E: 69e9, nu: 0.33, rho: 2700, sigmaY: 240e6 }, + { id: 'titanium', name: 'Titanium Ti-6Al-4V', E: 116e9, nu: 0.32, rho: 4500, sigmaY: 830e6 }, +] + +export function bucklingMaterialById(id: string): BMaterial { + return BUCKLING_MATERIALS.find((m) => m.id === id) ?? BUCKLING_MATERIALS[0] +} + +export interface BucklingScenario { + id: string + name: string + blurb: string + /** Textbook end-condition label. */ + ends: string + /** + * Effective-length factor K for the Euler reference P_cr = π²EI/(KL)². + * `null` for cases (panels) without a 1-D Euler analogue. + */ + kEff: number | null + /** Aspect the studio uses to pick a sensible default slenderness. */ + kind: 'column' | 'panel' +} + +export const BUCKLING_SCENARIOS: BucklingScenario[] = [ + { + id: 'cantilever', + name: 'Flagpole column', + blurb: 'Fixed base, free tip under axial compression — the fixed–free cantilever (K = 2).', + ends: 'fixed–free', + kEff: 2.0, + kind: 'column', + }, + { + id: 'fixed-pinned', + name: 'Braced column', + blurb: 'Clamped base, laterally-braced top under compression — the fixed–pinned strut (K ≈ 0.699).', + ends: 'fixed–pinned', + kEff: 0.6992, + kind: 'column', + }, + { + id: 'panel', + name: 'Compression panel', + blurb: 'A wide plate clamped along its loaded edges — buckles into a lateral bulge, not a 1-D column.', + ends: 'plate', + kEff: null, + kind: 'panel', + }, +] + +export function bucklingScenarioById(id: string): BucklingScenario { + return BUCKLING_SCENARIOS.find((s) => s.id === id) ?? BUCKLING_SCENARIOS[0] +} + +export interface BucklingParams { + scenarioId: string + materialId: string + /** Slenderness handle — column length as a multiple of its width b. */ + aspect: number + order: QOrder + nModes: number +} + +export const BUCKLING_DEFAULTS: BucklingParams = { + scenarioId: 'cantilever', + materialId: 'steel', + aspect: 45, + order: 8, + nModes: 3, +} + +/** Section width b (in-plane) and out-of-plane thickness t are fixed handles. */ +const SECTION_B = 0.02 // m +const SECTION_T = 0.01 // m +/** Reference compressive traction magnitude (unit stress), so λ reads directly. */ +const REF_TRACTION = 1 + +export interface BucklingGeometry { + L: number + b: number + t: number + area: number + I: number + radiusGyration: number + E: number + nu: number + kEff: number | null + /** Euler critical load for the ideal column (N), null for panels. */ + eulerLoad: number | null + /** Reference axial load actually applied (N) — λ × this = the buckling load. */ + refLoad: number + slenderness: number | null +} + +export interface BuiltBuckling { + input: BucklingInput + geom: BucklingGeometry + mesh: QuadMesh +} + +/** Column meshes are tall (axis = +y); panels are wide. Density tracks aspect. */ +function meshFor(scenario: BucklingScenario, order: QOrder, aspect: number) { + const b = SECTION_B + if (scenario.kind === 'panel') { + const W = b * 6 + const H = b * 4 + // Panels have clustered eigenvalues; keep the mesh coarse so the subspace + // solve stays interactive. + return { mesh: rectPlateQ(order, W, H, order === 8 ? 10 : 18, order === 8 ? 6 : 12), L: H, width: W } + } + const L = b * aspect + // Two elements across the width resolve the bending mode with Q8 (which does + // not shear-lock); a slender column needs several along its length. Q8 stays + // accurate to <0.3% with far fewer elements than Q4, so it can be much + // coarser — which keeps the eigenproblem interactive. + const ny = order === 8 ? Math.min(30, Math.max(16, Math.round(aspect * 0.6))) : Math.min(56, Math.max(28, aspect)) + const nx = order === 8 ? 2 : 6 + return { mesh: rectPlateQ(order, b, L, nx, ny), L, width: b } +} + +export function buildBuckling(params: BucklingParams): BuiltBuckling { + const scenario = bucklingScenarioById(params.scenarioId) + const mat = bucklingMaterialById(params.materialId) + const { mesh, L, width } = meshFor(scenario, params.order, params.aspect) + const t = SECTION_T + + const fix: BucklingInput['fix'] = [] + if (scenario.id === 'cantilever') { + fix.push({ edge: 'bottom', dofs: ['x', 'y'] }) + } else if (scenario.id === 'fixed-pinned') { + fix.push({ edge: 'bottom', dofs: ['x', 'y'] }) + fix.push({ edge: 'top', dofs: ['x'] }) // lateral brace → pin at the loaded end + } else { + // panel: clamp the bottom loaded edge, brace the top loaded edge laterally. + fix.push({ edge: 'bottom', dofs: ['x', 'y'] }) + fix.push({ edge: 'top', dofs: ['x'] }) + } + + const input: BucklingInput = { + mesh, + E: mat.E, + nu: mat.nu, + thickness: t, + fix, + traction: { edge: 'top', tx: 0, ty: -REF_TRACTION }, + nModes: params.nModes, + } + + // Section properties for in-plane (about-z) bending of a b×t rectangle. + const b = width + const area = b * t + const I = (t * b * b * b) / 12 + const radiusGyration = Math.sqrt(I / area) + const refLoad = REF_TRACTION * b * t // |traction| × loaded-edge area + const kEff = scenario.kEff + const euler = kEff != null ? eulerLoad(mat.E, I, L, kEff) : null + const slenderness = kEff != null ? (kEff * L) / radiusGyration : null + + return { + input, + mesh, + geom: { + L, + b, + t, + area, + I, + radiusGyration, + E: mat.E, + nu: mat.nu, + kEff, + eulerLoad: euler, + refLoad, + slenderness, + }, + } +} diff --git a/projects/keystone-fea-b3k9/src/engine/validate.ts b/projects/keystone-fea-b3k9/src/engine/validate.ts index 68065d44..1c9dfcd7 100644 --- a/projects/keystone-fea-b3k9/src/engine/validate.ts +++ b/projects/keystone-fea-b3k9/src/engine/validate.ts @@ -21,6 +21,7 @@ import { solveThermoelastic } from './thermoelastic' import { solveTransientThermoelastic } from './coupled' import { analyzeFracture, type CrackModel } from './fracture' import { edgeNodesQ } from './quadmesh' +import { solveBuckling as solveContinuumBuckling, eulerLoad, type BucklingInput } from './buckling' const STEEL_RHO = 7850 // kg/m³ @@ -1311,6 +1312,69 @@ export function runFractureBenchmarks(): Check[] { return checks } +/** + * Continuum linear-buckling benchmarks (v14). Each solves the plane-stress + * eigenproblem (K + λK_g)φ = 0 for a slender Q8 column and compares the FE + * critical load against the exact Euler value P_cr = π²EI/(KL)². The mode ladder + * (1 : 9 : 25 for a fixed–free column) checks the higher eigenpairs too. + */ +export function runBucklingBenchmarks(): Check[] { + const checks: Check[] = [] + const E = STEEL_E + const nu = 0.3 + const t = 0.01 + const b = 0.02 + const L = 1.0 + const I = (t * b * b * b) / 12 + const Pref = b * t // unit compressive traction × loaded-edge area + + // Fixed–free cantilever column (K = 2): the flagpole. Q8 kills shear locking. + { + const mesh = rectPlateQ(8, b, L, 2, 40) + const input: BucklingInput = { + mesh, + E, + nu, + thickness: t, + fix: [{ edge: 'bottom', dofs: ['x', 'y'] }], + traction: { edge: 'top', tx: 0, ty: -1 }, + nModes: 3, + } + const res = solveContinuumBuckling(input) + const Pcr = (res.modes[0]?.loadFactor ?? 0) * Pref + const Peuler = eulerLoad(E, I, L, 2.0) + checks.push(check('Euler column (fixed–free)', 'P_cr = π²EI/(2L)²', Peuler, Pcr, 'N', 1e-2)) + // Mode ladder: the second cantilever buckling load is 9× the first (3²/1²). + if (res.modes.length >= 2) { + const ratio = res.modes[1].loadFactor / res.modes[0].loadFactor + checks.push(check('Cantilever mode ratio', 'λ₂/λ₁ = 3²/1² = 9', 9, ratio, '', 3e-2)) + } + } + + // Fixed–pinned strut (K ≈ 0.6992): base clamped, top laterally braced. + { + const mesh = rectPlateQ(8, b, L, 2, 40) + const input: BucklingInput = { + mesh, + E, + nu, + thickness: t, + fix: [ + { edge: 'bottom', dofs: ['x', 'y'] }, + { edge: 'top', dofs: ['x'] }, + ], + traction: { edge: 'top', tx: 0, ty: -1 }, + nModes: 2, + } + const res = solveContinuumBuckling(input) + const Pcr = (res.modes[0]?.loadFactor ?? 0) * Pref + const Peuler = eulerLoad(E, I, L, 0.6992) + checks.push(check('Euler column (fixed–pinned)', 'P_cr = π²EI/(0.699L)²', Peuler, Pcr, 'N', 2e-2)) + } + + return checks +} + export function runAllBenchmarks(): { frame: Check[] dynamics: Check[] @@ -1324,6 +1388,7 @@ export function runAllBenchmarks(): { thermal: Check[] thermoelastic: Check[] fracture: Check[] + buckling: Check[] allPass: boolean } { const frame = runFrameBenchmarks() @@ -1338,6 +1403,7 @@ export function runAllBenchmarks(): { const thermal = runThermalBenchmarks() const thermoelastic = runThermoelasticBenchmarks() const fracture = runFractureBenchmarks() + const buckling = runBucklingBenchmarks() const allPass = [ ...frame, ...dynamics, @@ -1351,6 +1417,7 @@ export function runAllBenchmarks(): { ...thermal, ...thermoelastic, ...fracture, + ...buckling, ].every((c) => c.pass) - return { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, fracture, allPass } + return { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, fracture, buckling, allPass } } diff --git a/projects/keystone-fea-b3k9/src/state.ts b/projects/keystone-fea-b3k9/src/state.ts index 6947e7ae..d44618a7 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' | 'fracture' + tab: 'frame' | 'continuum' | 'topopt' | 'thermal' | 'fracture' | 'buckling' frame: FrameModel continuum: { presetId: string; density: number; elemOrder?: ElemOrder } display: Display diff --git a/projects/keystone-fea-b3k9/src/ui/BucklingStudio.tsx b/projects/keystone-fea-b3k9/src/ui/BucklingStudio.tsx new file mode 100644 index 00000000..fa9d5a0b --- /dev/null +++ b/projects/keystone-fea-b3k9/src/ui/BucklingStudio.tsx @@ -0,0 +1,670 @@ +// The Buckling studio: a self-contained tab that meshes a prismatic column (or +// panel) on Keystone's isoparametric machinery, solves the linear-buckling +// generalized eigenproblem (K + λK_g)φ = 0, and turns the lowest critical load +// factor into the two pictures engineers actually reason about: +// +// • the *buckled shape* — the pre-buckling straight column shaded by its axial +// load path, with the eigenmode animated as a breathing lateral bow, and +// • the *column curve* — the FE critical stress dropped onto the Euler +// hyperbola σ_cr = π²E/λ² with its yield cut-off, so you can read at a glance +// whether the member fails by buckling or by squashing. +// +// Hooks discipline mirrors the other studios: 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 { solveBuckling, type BucklingResult } from '../engine/buckling' +import { + BUCKLING_SCENARIOS, + BUCKLING_MATERIALS, + bucklingScenarioById, + buildBuckling, + BUCKLING_DEFAULTS, + type BucklingParams, + type BucklingGeometry, +} from '../engine/bucklingpresets' +import type { QuadMesh } from '../engine/quadmesh' +import { fieldColor, rgbStr, type Colormap } from './colormap' +import { Segmented, Slider, StatTile, Toggle } from './components' +import { fmtEng } from './format' + +type ViewMode = 'shape' | 'curve' + +interface Snapshot { + res: BucklingResult | null + mesh: QuadMesh + geom: BucklingGeometry + view: ViewMode + modeSel: number + colormap: Colormap + showMesh: boolean + playing: boolean +} + +function fit(minX: number, maxX: number, minY: number, maxY: number, cw: number, ch: number, margin = 46) { + const w = maxX - minX || 1 + const h = maxY - minY || 1 + const s = Math.min((cw - 2 * margin) / w, (ch - 2 * margin) / h) + const ox = (cw - w * s) / 2 - minX * s + const oy = (ch - h * s) / 2 + maxY * s + return { s, toX: (x: number) => ox + x * s, toY: (y: number) => oy - y * s } +} + +/** Corner-node local indices for a Q4/Q8 element (first four are the corners). */ +const CORNERS = [0, 1, 2, 3] + +export function BucklingStudio() { + const [params, setParams] = useState(BUCKLING_DEFAULTS) + const [view, setView] = useState('shape') + const [modeSel, setModeSel] = useState(0) + const [colormap] = useState('viridis') + const [showMesh, setShowMesh] = useState(true) + const [playing, setPlaying] = useState(true) + + // The slenderness slider drags continuously, but each solve is a full + // eigenproblem — debounce the commit so we solve once the drag settles, not + // 50 times mid-drag. + const [aspectDraft, setAspectDraft] = useState(params.aspect) + useEffect(() => { + if (aspectDraft === params.aspect) return + const id = setTimeout(() => setParams((prev) => ({ ...prev, aspect: aspectDraft })), 220) + return () => clearTimeout(id) + }, [aspectDraft, params.aspect]) + + const scenario = useMemo(() => bucklingScenarioById(params.scenarioId), [params.scenarioId]) + const built = useMemo(() => buildBuckling(params), [params]) + + const res = useMemo(() => { + try { + return solveBuckling(built.input) + } catch { + return null + } + }, [built]) + + // Clamp the selected mode to what actually came back. + const effMode = res && res.modes.length ? Math.min(modeSel, res.modes.length - 1) : 0 + const mode = res?.modes[effMode] ?? null + + // --- FE vs Euler numbers ------------------------------------------------ + const stats = useMemo(() => { + if (!res || !mode) return null + const g = built.geom + const Pcr = mode.loadFactor * g.refLoad + const sigmaCr = Pcr / g.area + const euler = g.eulerLoad + const sigmaEuler = euler != null ? euler / g.area : null + const relErr = euler != null && euler > 0 ? Math.abs(Pcr - euler) / euler : null + return { Pcr, sigmaCr, euler, sigmaEuler, relErr } + }, [res, mode, built.geom]) + + // --- snapshot for the rAF loop ----------------------------------------- + const canvasRef = useRef(null) + const wrapRef = useRef(null) + const sizeRef = useRef({ w: 820, h: 540 }) + const snapRef = useRef({ + res, + mesh: built.mesh, + geom: built.geom, + view, + modeSel: effMode, + colormap, + showMesh, + playing, + }) + useEffect(() => { + snapRef.current = { + res, + mesh: built.mesh, + geom: built.geom, + view, + modeSel: effMode, + colormap, + showMesh, + playing, + } + }, [res, built.mesh, built.geom, view, effMode, colormap, showMesh, playing]) + + 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 + 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.7 + 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 === 'curve') drawCurve(ctx, cw, ch, snap) + else drawShape(ctx, cw, ch, snap, phase) + } + } + raf = requestAnimationFrame(loop) + } + raf = requestAnimationFrame(loop) + return () => cancelAnimationFrame(raf) + }, []) + + const patch = (p: Partial) => setParams((prev) => ({ ...prev, ...p })) + const setScenario = (id: string) => { + const sc = bucklingScenarioById(id) + setModeSel(0) + setParams((prev) => { + const aspect = sc.kind === 'panel' ? prev.aspect : Math.max(prev.aspect, 30) + setAspectDraft(aspect) + return { ...prev, scenarioId: id, aspect } + }) + if (sc.kind === 'panel') setView('shape') + setPlaying(true) + } + + const nModes = res?.modes.length ?? 0 + + return ( +
+ + +
+
+ + options={[ + { value: 'shape', label: '⟂ Buckled shape' }, + { value: 'curve', label: '↘ Column curve' }, + ]} + value={view} + onChange={setView} + /> + {nModes > 0 && view === 'shape' && ( +
+ {res!.modes.map((m, i) => ( + + ))} +
+ )} + {view === 'shape' && ( + + )} +
+ {!res + ? 'solving eigenproblem…' + : view === 'curve' + ? 'critical stress vs slenderness — Euler hyperbola, yield cut-off, and this design point' + : mode + ? `mode ${effMode + 1} · λ_cr = ${mode.loadFactor.toExponential(3)} · ${res.freeDofCount} DOF · ${res.eigIterations} subspace iters` + : 'no buckling mode found'} +
+
+
+ +
+
+ + +
+ ) +} + +function materialYield(id: string): number { + return BUCKLING_MATERIALS.find((m) => m.id === id)?.sigmaY ?? 250e6 +} + +// --- drawing --------------------------------------------------------------- + +/** Element corner polygons for a deformed configuration. */ +function drawMeshPolys( + ctx: CanvasRenderingContext2D, + mesh: QuadMesh, + X: (n: number) => number, + Y: (n: number) => number, + fill: (e: number) => string | null, + stroke: string | null, + lw: number, +) { + const order = mesh.order + for (let e = 0; e < mesh.elemCount; e++) { + const base = e * order + ctx.beginPath() + for (let k = 0; k < CORNERS.length; k++) { + const n = mesh.elems[base + CORNERS[k]] + const px = X(n) + const py = Y(n) + if (k === 0) ctx.moveTo(px, py) + else ctx.lineTo(px, py) + } + ctx.closePath() + const f = fill(e) + if (f) { + ctx.fillStyle = f + ctx.fill() + } + if (stroke) { + ctx.strokeStyle = stroke + ctx.lineWidth = lw + ctx.stroke() + } + } +} + +function drawShape(ctx: CanvasRenderingContext2D, cw: number, ch: number, snap: Snapshot, phase: number) { + ctx.clearRect(0, 0, cw, ch) + ctx.fillStyle = '#0a0e17' + ctx.fillRect(0, 0, cw, ch) + const { mesh, res } = snap + const md = res?.modes[snap.modeSel] ?? null + + // Amplitude of the animated bow. Scale so the peak lateral deflection is a + // healthy fraction of the model size; breathe with the animation phase. + const span = Math.max(mesh.maxX - mesh.minX, mesh.maxY - mesh.minY) + const amp = md ? Math.sin(phase * 2 * Math.PI) * span * 0.16 : 0 + + // Fit the *deformed* extent so the bow never clips. + const pad = md ? span * 0.18 : 0 + const t = fit(mesh.minX - pad, mesh.maxX + pad, mesh.minY, mesh.maxY, cw, ch) + + const dX = (n: number) => mesh.x[n] + (md ? amp * md.dispX[n] : 0) + const dY = (n: number) => mesh.y[n] + (md ? amp * md.dispY[n] : 0) + + // Reference axial-stress load path (magnitude of compressive σyy). + let smax = 1e-30 + if (res) for (let n = 0; n < mesh.nodeCount; n++) smax = Math.max(smax, Math.abs(res.refSyy[n])) + + // Undeformed ghost. + drawMeshPolys( + ctx, + mesh, + (n) => t.toX(mesh.x[n]), + (n) => t.toY(mesh.y[n]), + () => 'rgba(120,140,175,0.05)', + 'rgba(120,140,175,0.14)', + 0.6, + ) + + // Deformed (buckled) shape, shaded by the reference axial load path. + const order = mesh.order + drawMeshPolys( + ctx, + mesh, + (n) => t.toX(dX(n)), + (n) => t.toY(dY(n)), + (e) => { + if (!res) return 'rgba(90,150,230,0.5)' + let s = 0 + for (let k = 0; k < CORNERS.length; k++) s += Math.abs(res.refSyy[mesh.elems[e * order + CORNERS[k]]]) + s /= CORNERS.length + const c = fieldColor(Math.min(1, s / smax), snap.colormap) + return rgbStr(c, 0.9) + }, + snap.showMesh ? 'rgba(10,14,23,0.55)' : null, + 0.5, + ) + + // Load arrows on the loaded (top) edge. + ctx.fillStyle = '#ff6b6b' + ctx.strokeStyle = '#ff6b6b' + ctx.lineWidth = 1.5 + const topY = t.toY(mesh.maxY + pad * 0.15) + const nArr = 7 + for (let i = 0; i <= nArr; i++) { + const fx = mesh.minX + ((mesh.maxX - mesh.minX) * i) / nArr + const sx = t.toX(fx) + ctx.beginPath() + ctx.moveTo(sx, topY - 22) + ctx.lineTo(sx, topY - 4) + ctx.stroke() + ctx.beginPath() + ctx.moveTo(sx - 3, topY - 9) + ctx.lineTo(sx, topY - 3) + ctx.lineTo(sx + 3, topY - 9) + ctx.fill() + } + + // Support hatching at the base. + ctx.strokeStyle = 'rgba(200,210,230,0.6)' + ctx.lineWidth = 1 + const baseY = t.toY(mesh.minY) + const x0 = t.toX(mesh.minX) + const x1 = t.toX(mesh.maxX) + ctx.beginPath() + ctx.moveTo(x0, baseY) + ctx.lineTo(x1, baseY) + ctx.stroke() + for (let x = x0; x < x1; x += 9) { + ctx.beginPath() + ctx.moveTo(x, baseY) + ctx.lineTo(x - 7, baseY + 7) + ctx.stroke() + } + + // Caption. + ctx.fillStyle = 'rgba(210,220,240,0.85)' + ctx.font = '12px ui-monospace, monospace' + ctx.textAlign = 'left' + if (md) { + ctx.fillText( + `buckling mode ${snap.modeSel + 1} λ_cr = ${md.loadFactor.toExponential(3)}`, + 14, + 22, + ) + } +} + +function drawCurve(ctx: CanvasRenderingContext2D, cw: number, ch: number, snap: Snapshot) { + ctx.clearRect(0, 0, cw, ch) + ctx.fillStyle = '#0a0e17' + ctx.fillRect(0, 0, cw, ch) + const { geom, res } = snap + const E = geom.E + + const pad = { l: 64, r: 22, t: 26, b: 46 } + const plotW = cw - pad.l - pad.r + const plotH = ch - pad.t - pad.b + + // Axes: slenderness 0..λmax, stress 0..σy·1.18. Stretch the slenderness axis + // so the current design point always lands inside the plot. + const sy = materialYieldFromE(geom) + const lambdaMax = Math.max(220, (geom.slenderness ?? 0) * 1.15) + const sMax = sy * 1.18 + + const xOf = (lam: number) => pad.l + (Math.min(lam, lambdaMax) / lambdaMax) * plotW + const yOf = (s: number) => pad.t + (1 - Math.min(s, sMax) / sMax) * plotH + + // Grid + axes. + ctx.strokeStyle = 'rgba(120,140,175,0.16)' + ctx.lineWidth = 1 + ctx.fillStyle = 'rgba(180,195,220,0.7)' + ctx.font = '11px ui-monospace, monospace' + ctx.textAlign = 'right' + for (let i = 0; i <= 4; i++) { + const s = (sMax * i) / 4 + const y = yOf(s) + ctx.beginPath() + ctx.moveTo(pad.l, y) + ctx.lineTo(cw - pad.r, y) + ctx.stroke() + ctx.fillText(`${(s / 1e6).toFixed(0)}`, pad.l - 8, y + 4) + } + ctx.textAlign = 'center' + for (let i = 0; i <= 4; i++) { + const lam = (lambdaMax * i) / 4 + const x = xOf(lam) + ctx.beginPath() + ctx.moveTo(x, pad.t) + ctx.lineTo(x, ch - pad.b) + ctx.stroke() + ctx.fillText(`${lam.toFixed(0)}`, x, ch - pad.b + 18) + } + ctx.fillStyle = 'rgba(210,220,240,0.85)' + ctx.fillText('slenderness λ = KL / r', pad.l + plotW / 2, ch - 10) + ctx.save() + ctx.translate(16, pad.t + plotH / 2) + ctx.rotate(-Math.PI / 2) + ctx.fillText('critical stress σ_cr (MPa)', 0, 0) + ctx.restore() + + // Yield plateau. + ctx.strokeStyle = 'rgba(255,180,90,0.85)' + ctx.setLineDash([6, 4]) + ctx.lineWidth = 1.5 + ctx.beginPath() + ctx.moveTo(xOf(0), yOf(sy)) + ctx.lineTo(xOf(lambdaMax), yOf(sy)) + ctx.stroke() + ctx.setLineDash([]) + ctx.fillStyle = 'rgba(255,180,90,0.9)' + ctx.textAlign = 'left' + ctx.fillText('σ_y (squash)', xOf(lambdaMax) - 96, yOf(sy) - 6) + + // Euler hyperbola σ_cr = π²E/λ². + ctx.strokeStyle = 'rgba(90,180,255,0.95)' + ctx.lineWidth = 2 + ctx.beginPath() + let started = false + for (let lam = 1; lam <= lambdaMax; lam += 1) { + const s = (Math.PI * Math.PI * E) / (lam * lam) + if (s > sMax) continue + const x = xOf(lam) + const y = yOf(s) + if (!started) { + ctx.moveTo(x, y) + started = true + } else ctx.lineTo(x, y) + } + ctx.stroke() + ctx.fillStyle = 'rgba(90,180,255,0.95)' + ctx.fillText('Euler π²E/λ²', xOf(lambdaMax * 0.5), yOf((Math.PI * Math.PI * E) / (lambdaMax * 0.5) ** 2) - 8) + + // Transition slenderness λ₁ = π√(E/σy). + const lam1 = Math.PI * Math.sqrt(E / sy) + ctx.strokeStyle = 'rgba(200,210,230,0.35)' + ctx.setLineDash([3, 4]) + ctx.beginPath() + ctx.moveTo(xOf(lam1), pad.t) + ctx.lineTo(xOf(lam1), ch - pad.b) + ctx.stroke() + ctx.setLineDash([]) + + // The design point(s). + if (res && geom.slenderness != null && geom.eulerLoad != null) { + const lam = geom.slenderness + const md = res.modes[Math.min(snap.modeSel, res.modes.length - 1)] + const sigmaFE = md ? (md.loadFactor * geom.refLoad) / geom.area : 0 + const sigmaEuler = geom.eulerLoad / geom.area + // Euler point (hollow). + ctx.strokeStyle = 'rgba(90,180,255,1)' + ctx.lineWidth = 2 + ctx.beginPath() + ctx.arc(xOf(lam), yOf(sigmaEuler), 6, 0, 2 * Math.PI) + ctx.stroke() + // FE point (filled). + ctx.fillStyle = '#7CFFB2' + ctx.beginPath() + ctx.arc(xOf(lam), yOf(sigmaFE), 5, 0, 2 * Math.PI) + ctx.fill() + ctx.strokeStyle = 'rgba(10,14,23,0.8)' + ctx.lineWidth = 1 + ctx.stroke() + // Callout. + ctx.fillStyle = 'rgba(220,255,235,0.95)' + ctx.textAlign = 'left' + ctx.font = '12px ui-monospace, monospace' + const lx = Math.min(xOf(lam) + 12, cw - pad.r - 150) + ctx.fillText(`FE σ_cr = ${(sigmaFE / 1e6).toFixed(1)} MPa`, lx, yOf(sigmaFE) - 10) + ctx.fillStyle = 'rgba(150,200,255,0.9)' + ctx.fillText(`Euler ${(sigmaEuler / 1e6).toFixed(1)} MPa`, lx, yOf(sigmaFE) + 6) + } else { + ctx.fillStyle = 'rgba(180,195,220,0.7)' + ctx.textAlign = 'center' + ctx.fillText('panel scenario — no single-slenderness column point', pad.l + plotW / 2, pad.t + plotH / 2) + } +} + +/** The yield stress belongs to the material; the curve module only has geometry, + * so we recover σ_y from the material via the studio-level lookup embedded in E. */ +function materialYieldFromE(geom: BucklingGeometry): number { + const m = BUCKLING_MATERIALS.find((mm) => Math.abs(mm.E - geom.E) < 1) + return m?.sigmaY ?? 250e6 +} diff --git a/projects/keystone-fea-b3k9/src/ui/components.tsx b/projects/keystone-fea-b3k9/src/ui/components.tsx index c1531939..75313543 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, fracture, allPass } = useMemo( + const { frame, dynamics, harmonic, seismic, plastic, inelastic, continuum, quad, topopt, thermal, thermoelastic, fracture, buckling, allPass } = useMemo( () => runAllBenchmarks(), [], ) const [open, setOpen] = useState(false) - const all = [...frame, ...dynamics, ...harmonic, ...seismic, ...plastic, ...inelastic, ...continuum, ...quad, ...topopt, ...thermal, ...thermoelastic, ...fracture] + const all = [...frame, ...dynamics, ...harmonic, ...seismic, ...plastic, ...inelastic, ...continuum, ...quad, ...topopt, ...thermal, ...thermoelastic, ...fracture, ...buckling] const passCount = all.filter((c) => c.pass).length return (
@@ -946,6 +946,10 @@ export function VerifyBadge() { {fracture.map((c, i) => ( ))} +
Buckling — Euler critical loads
+ {buckling.map((c, i) => ( + + ))}
)}