/ proof-of-concept / plate-mechanics-sim
Sign in
proof-of-concept / plate-mechanics-sim
main
Branches
main
plate-mechanics-sim / packages / core / src / methods / fd.ts
90 lines · 4.8 KBCodeBlameHistory
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 1import type { MethodDefinition } from "./types";
78d40e9Free edges for both finite-difference methodsJeremy Magland 4 FREE_EDGE_EQUATIONS,
5 GOVERNING_EQUATIONS,
78d40e9Free edges for both finite-difference methodsJeremy Magland 7 RADIATION_EQUATIONS,
8 STRESS_EQUATIONS,
9} from "./common";
11export const FD_METHOD: MethodDefinition = {
12 id: "fd",
13 name: "Explicit finite differences",
14 shortName: "Explicit FD",
15 summary: "Explicit second-order time stepping of the 13-point biharmonic stencil.",
16 solverFunction: "plate_fd",
17 planFunction: "plate_fd_plan",
18 supports: {
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 19 shape: ["rectangle", "ellipse"],
20 bc: ["simply-supported", "clamped", "free", "pinned-center"],
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 21 excitation: ["impulse", "harmonic"],
22 audio: "clip",
23 },
24 readingOrder: [
25 "plate_fd.m",
27 "plate_fd_loop.m",
28 "plate_fd_output.m",
30 "plate_setup.m",
33 "plate_laplacian.m",
34 "plate_curvature.m",
78d40e9Free edges for both finite-difference methodsJeremy Magland 35 "plate_matrices.m",
36 "plate_rigid_modes.m",
37 "plate_curvature_free.m",
40 "plate_force_footprint.m",
41 "plate_stress.m",
42 "plate_listener.m",
43 "plate_spectrum.m",
44 "plate_result_schema.m",
45 ],
46 equations: [
47 ...GOVERNING_EQUATIONS,
48 {
49 tex: String.raw`(1+\sigma_0\Delta t)\,w^{n+1} = 2w^n - (1-\sigma_0\Delta t)\,w^{n-1} - \kappa^2\Delta t^2\, \mathsf L\mathsf L\, w^n + 2\sigma_1\Delta t\,\mathsf L\!\left(w^n - w^{n-1}\right) + \frac{\Delta t^2}{\mu} F^n s`,
50 caption:
51 "Explicit scheme, second order in space and time; L is the five-point Laplacian and LL the resulting 13-point biharmonic stencil.",
52 },
53 {
54 tex: String.raw`\mathsf L w_{ij} = \frac{w_{i,j+1} - 2w_{ij} + w_{i,j-1}}{\Delta x^2} + \frac{w_{i+1,j} - 2w_{ij} + w_{i-1,j}}{\Delta y^2}`,
55 caption:
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 56 "On a rectangle, supported and clamped edges enter through one ghost layer: odd reflection for simply supported (w = w_nn = 0), even reflection for clamped (w = w_n = 0). Free edges, and every condition on an ellipse, take the sparse operators below.",
58 {
59 tex: String.raw`\Delta t \le \frac{2}{\lambda_{\max}\left(\sigma_1 + \sqrt{\kappa^2 + \sigma_1^2}\right)}, \qquad \lambda_{\max} = \frac{4}{\Delta x^2} + \frac{4}{\Delta y^2}`,
60 caption:
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 61 "Stability bound from von Neumann analysis; for a square grid without sigma1 it reads dt <= dx^2 / (4 kappa). Because dt scales with dx^2, halving the spacing costs eight times more work. The sparse operators of the free plate and of the ellipse stay just under the same bound.",
78d40e9Free edges for both finite-difference methodsJeremy Magland 63 ...FREE_EDGE_EQUATIONS,
65 ...CURVED_EDGE_EQUATIONS,
67 ...RADIATION_EQUATIONS,
68 ],
69 validity: [
70 "Second-order truncation error in dx, dy and dt; the top modes on the grid are dispersive and their frequencies are too low.",
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 71 "The probe is snapped to the nearest grid node with a full stencil inside the plate.",
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 72 "Stress fields are evaluated only at the saved frames, so the envelope samples the true peak.",
73 "The time-step limit makes long simulated times expensive: cost grows as (1/dx)^4 for a fixed duration.",
74 "Audio is a short clip covering only the simulated window, box-averaged onto audio samples.",
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 75 "With the GPU option the time loop runs in single precision (a WebGPU port of plate_fd_loop.m); the setup, the stresses, the audio summation and the spectra stay in double precision in numbl. Against the double-precision loop the displacement differs by a relative L2 error of about 1e-6 over 50 ms on a 61 x 41 grid, the stress and the sound by about 1e-4.",
78d40e9Free edges for both finite-difference methodsJeremy Magland 76 "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
81f3b30Elliptical plates, a WebGPU time loop for the explicit method, and WAV exportJeremy Magland 77 "Pinned centre: the support is one grid node, so the stress right at it is resolution-dependent, and the two rigid rotations about it are removed rather than restrained by a modelled washer.",
78 "Ellipse: the edge is a staircase of grid nodes, so the geometry, the edge conditions and the stress along the edge are first order in the grid spacing. The stress on the outermost ring of nodes is copied from the nearest node with a full stencil, which puts the peak stress of a clamped edge about one spacing inside. For simply supported edges the curvature term that the staircase loses is added back (the plate paradox), so the disk's frequencies converge to the exact ones.",
80 durationOf: (spec) => spec.fd.duration,
81 gridOf: (spec) => [spec.fd.nx, spec.fd.ny],
82 settingsKey: "fd",
83 usPerCellStep: 0.35,
85 loopFile: "plate_fd_loop.m",
86 setupFunction: "plate_fd_setup",
87 outputFunction: "plate_fd_output",
88 usPerCellStep: 0.02,
89 },
↑↓ move↵openescclose