/ proof-of-concept / plate-mechanics-sim
Sign in
proof-of-concept / plate-mechanics-sim
plate-mechanics-sim / README.md
75 lines · 9.0 KBPreviewCodeBlameHistoryRaw
1# Plate Mechanics Lab
3An interactive lab for the mechanics of thin metal plates. You specify geometry (a rectangle, or an ellipse such as a disk), material, edge conditions and an excitation; the app solves the damped Kirchhoff-Love equation and shows displacement, surface stress (von Mises and principal) over time, a probe time series and spectrum, and the radiated sound at a listener point, which you can play or save as a WAV file.
5The numerics run in [numbl](https://github.com/flatironinstitute/numbl), a MATLAB-syntax environment executing in the browser. The MATLAB source shown on screen is the source that ran: edit any file and rerun. No physics is computed in JavaScript; the host only draws arrays and plays samples.
7## Layout
9```
10packages/core problem spec, .m generation, numbl runner, result typing, field and colormap primitives (no UI)
11packages/core/src/m the MATLAB library: plate_modal.m, plate_fd.m, plate_stress.m, ...
12packages/components React components: code view, parameter panel, field view, plots, equations, run status, audio
13apps/lab the app: one arrangement of the components
14apps/tutorials MDX documents with embedded interactive views (one tutorial: the CFL limit)
15```
17The app and the tutorials are peers over the same core and components.
19## Methods
21Three methods solve the same problem specification, and the app can run any two and compare them directly:
23- **Modal superposition** (`plate_modal.m`): closed-form modes of the simply supported rectangle, each integrated exactly in time for piecewise-linear forcing. No time-step limit; the approximation is the modal truncation.
24- **Explicit finite differences** (`plate_fd.m`): second-order scheme with the 13-point biharmonic stencil; simply supported, clamped or free edges, or free edges pinned at the centre like a cymbal on its stand (the last two use the energy-method operators of `plate_matrices.m`). Stable only below the CFL bound from `plate_fd_plan.m`; the slider goes past it on purpose. The work is split into `plate_fd_setup.m`, `plate_fd_loop.m` and `plate_fd_output.m`, and the loop has a WebGPU port (below).
25- **Implicit finite differences** (`plate_implicit.m`): the same stencils assembled as sparse matrices (`plate_matrices.m`), stiffness averaged over three time levels (theta scheme, theta = 1/4 is the trapezoidal rule). One sparse Cholesky factorization, two triangular solves per step, no stability limit; the step is chosen for accuracy. Needs numbl 0.4.21 or later for its sparse direct solvers.
27A method is a data record (`packages/core/src/methods/*.ts`) plus a `.m` file. Adding one means adding both; nothing else in core changes.
29## Shapes
31The plate is a rectangle or the ellipse inscribed in the same Lx by Ly box (a disk when Lx = Ly). The finite-difference methods embed the ellipse in the rectangular grid: `plate_domain.m` marks the nodes within half a spacing of the inside, so the edge is a staircase that straddles the true curve, and `plate_matrices.m` assembles every edge condition on that footprint as a sparse operator (for the rectangle it reproduces the old operators exactly). Supported and clamped edges hold the outermost ring of nodes at zero, the clamped ring with the mirrored Laplacian that gives the 7/dx^4 coefficient on a straight edge; free edges use the energy method on every node of the plate. The geometry, the edge conditions and the stress along the edge are then first order in the spacing, against second order inside, and the run notes say so.
33A simply supported curved edge needs one more term. With w = 0 on the edge the bending energy is D/2 [∫(∇²w)² dA − (1 − ν)∮ κ w_n² ds], with κ the curvature of the edge; a polygon or a staircase has κ = 0, so a scheme that only imposes w = 0 converges to a plate with zero Laplacian on the edge rather than zero moment (the Babuška, or Sapondzhyan, paradox). For a disk that puts the fundamental at λ² = 5.78 instead of 4.94 (ν = 0.3). `plate_matrices.m` adds the boundary term with the curvature of the true ellipse at each edge node, and the test suite checks the disk against the Bessel-function frequency equations: clamped 10.216, simply supported 4.935, free 5.358 and 9.003, all within a few percent on a 25-node grid and converging under refinement. The modal method has closed-form modes only for the simply supported rectangle and refuses the ellipse.
35## The GPU option
37The explicit scheme is many cheap steps, which is what a GPU does well, so its time loop has a WebGPU port: `packages/core/src/gpu/fdLoop.ts` is `plate_fd_loop.m` as compute shaders, three dispatches per step for a rectangle with supported or clamped edges (Laplacian, update, reductions for the probe, the audio bins and the blow-up test), three with the sparse operators of an ellipse, and five when the rigid-body projection of a free plate is added. A checkbox in the explicit method's settings turns it on when the browser exposes WebGPU. The setup (`plate_fd_setup.m`) and the output (`plate_fd_output.m`, stresses, sound, spectra) still run in numbl, unchanged and in double precision; the worker hands the setup struct to the GPU and the loop's record back. The generated script says which way it ran, the run details name the adapter, and editing `plate_fd_loop.m` sends the loop back to numbl, where the edit is live. This follows the pattern of the cymbal2 and turing-surface projects: the `.m` files remain the source of truth, and the port must compute the same record from the same struct.
39The GPU works in single precision. The state is carried as (w, w - w_prev) rather than (w, w_prev), the same arithmetic rearranged so that the rounding of w does not feed back into the velocity; measured against the double-precision loop over 50 ms on the 61 x 41 grid, the displacement differs by a relative L2 error of about 1e-6 and the stress and the sound by about 1e-4. The cost is about 20 µs per step on an integrated Intel GPU, nearly independent of the grid at these sizes: the 61 x 41 grid for 50 ms (27,000 steps) takes 0.8 s in total against 22 s in numbl, and a 121 x 81 grid (107,000 steps) takes 3 s. The implicit method is not ported: its step is two sparse triangular solves, which have no parallelism to offer at this size.
41## The live app
43The lab is deployed at <https://plate-mechanics-sim.vault1-sites.magland.org/>. `.mochi/workflows/deploy.yml`
44builds `apps/lab` on every push to `main` and publishes it; the tutorials are not part of that artifact. The
45site sits at the root of an origin of its own, so the build's base path comes from `configure-pages` rather
46than being written into the Vite config.
48## Running
50```bash
51npm install
52npm run dev # the lab, http://localhost:5173
53npm run dev:tutorials # the tutorials
54npm test # vitest: runs the .m files in numbl under Node and checks modal vs FD, CFL, frequencies, editing
55npm run typecheck
56npm run build
57```
59`packages/core/test/gpu.test.ts` runs the WebGPU loop from Node against desktop WebGPU through the optional `webgpu` package (prebuilt Google Dawn) and compares it with the numbl loop for every edge condition on both shapes; `npm i --no-save webgpu` installs it. Without the package or a GPU adapter the suite skips itself.
61The `.m` files also run unchanged from the numbl CLI, Octave or MATLAB:
63```bash
64npx numbl run --path packages/core/src/m my_script.m
65```
67## Honest limits
69- Kirchhoff-Love thin-plate theory: no shear deformation or rotary inertia, so thick plates and very high modes are wrong.
70- Edge conditions are idealizations. Free edges are handled by a discrete bending energy whose natural boundary conditions are first-order accurate at the edge, and the rigid-body motion of a free plate is removed rather than modelled, as the supports of a real free-free rig would do. The centre-pinned plate is held at one grid node, so the stress right at the support is resolution-dependent and the rocking a real stand's washers would restrain is removed instead.
71- An ellipse is a staircase on the grid. Its edge is first order in the grid spacing, the stress on the outermost ring of nodes is copied from the nearest node with a full stencil (so the peak stress of a clamped edge sits about one spacing inside), and degenerate mode pairs of a disk are split slightly by the grid's anisotropy. Refine the grid to check.
72- Damping is a two-parameter model, chosen, not measured.
73- Stress is a second derivative of the displacement and converges more slowly than displacement in both methods; near the strike it depends on the load footprint.
74- Radiated sound is Rayleigh's integral in an infinite baffle: an estimate of the waveform, not a calibrated level.
75- Explicit finite differences cost about 0.3 µs per cell-step in numbl, and the stable step scales with dx². The default 41 x 27 grid gives 20 ms in a few seconds; refining to 61 x 41 costs about five times more. The GPU option brings this to about 20 µs per step whatever the grid, in single precision. The implicit scheme removes the step limit (about 1 ms per step on 61 x 41) at the price of detuning the highest modes; long simulated times still belong to the modal method.
↑↓ move↵openescclose