f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 1// Physical and numerical parameters of the simulation, plus derived quantities.
3export interface PhysParams {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 4 shape: PlateShape; // outline of the sheet within its bounding box
5 Lx: number; // plate width along x [m] (bounding box for an ellipse)
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 6 Ly: number; // plate height along y [m]
7 h: number; // thickness [m]
8 E: number; // Young's modulus [Pa]
9 nu: number; // Poisson ratio
10 rho: number; // density [kg/m^3]
11 sigma0: number; // frequency-independent loss [1/s]
12 sigma1: number; // frequency-dependent loss [m^2/s]
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 13 sigmaIn: number; // additional loss of the in-plane motion [1/s]
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 14 holeInset: number; // distance of the thread hole from the two edges of the top-left corner [m]
15 pushAmount: number; // corner displacement [m]
16 fingerRadius: number; // radius of the pushed patch [m]
17 rampTime: number; // duration of the push ramp [s]
18 holdMax: number; // maximum settle time before release [s]
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 19 releaseTime: number; // time over which the finger force falls to zero [s]
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 20 settleSigma0: number; // artificial damping during push and settle [1/s]
21 settleKE: number; // kinetic energy threshold that counts as settled [J]
22 spacing: number; // requested grid spacing [m]
23 fs: number; // audio sample rate [Hz]
24 fps: number; // movie frame rate (also the batch rate) [Hz]
25 duration: number; // recorded duration after release [s]
26 listenerDist: number; // [m]
27 listenerAzimuth: number; // [deg], rotation about the plate's vertical axis (0 = on the normal)
28 listenerElevation: number; // [deg]
29 radiation: "rayleigh" | "dipole";
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 30 bowShape: BowShape; // stress-free initial shape of the sheet
31 bowAmplitude: number; // sag of the initial shape, center relative to the edges [m]
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 32 excitation: Excitation;
33 strikeFracX: number; // strike point as a fraction of the width from the left edge
34 strikeFracY: number; // strike point as a fraction of the height from the bottom edge
35 strikeVelocity: number; // stick tip speed at impact [m/s]
36 stickMass: number; // effective mass of the stick [kg]
37 contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
38 contactExponent: number; // Hertz exponent alpha (1.5 for a sphere)
39 contactRadius: number; // radius of the contact patch [m]
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 42export type Excitation = "push" | "strike";
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 44/** A rectangle Lx by Ly, or the ellipse inscribed in that rectangle. */
45export type PlateShape = "rectangle" | "ellipse";
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 47export type BowShape = "none" | "width" | "height" | "dish";
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 49export const defaultParams: PhysParams = {
52 Ly: 18 * 0.0254,
54 E: 193e9,
55 nu: 0.29,
56 rho: 7900,
57 sigma0: 1.15,
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 58 sigma1: 5e-4,
59 sigmaIn: 20,
61 pushAmount: 0.01,
62 fingerRadius: 0.01,
63 rampTime: 0.1,
64 holdMax: 0.5,
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 65 releaseTime: 0.02,
67 settleKE: 1e-8,
70 fps: 120,
71 duration: 4,
72 listenerDist: 1.0,
73 listenerAzimuth: 0,
74 listenerElevation: 0,
75 radiation: "rayleigh",
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 76 bowShape: "width",
77 bowAmplitude: 3e-3,
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 78 excitation: "push",
79 strikeFracX: 0.7,
80 strikeFracY: 0.35,
81 strikeVelocity: 3,
82 stickMass: 0.03,
83 contactStiffness: 5e8,
84 contactExponent: 1.5,
85 contactRadius: 5e-3,
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 88/** Strike point in plate coordinates (relative to the center). */
89export function strikePoint(p: PhysParams): [number, number] {
90 return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
91}
94 * Center of the thread hole. Rectangle: inset from the two edges of the
95 * top-left corner. Ellipse: inset along the radius from the boundary point
96 * in the top-left direction (parametric angle 135 degrees).
97 */
98export function holePoint(p: PhysParams): [number, number] {
99 if (p.shape === "ellipse") {
100 const c = Math.SQRT1_2;
101 return [-(p.Lx / 2 - p.holeInset) * c, (p.Ly / 2 - p.holeInset) * c];
102 }
103 return [-p.Lx / 2 + p.holeInset, p.Ly / 2 - p.holeInset];
104}
106/** Center of the pushed patch: the bottom-right corner, or the boundary point opposite the hole. */
107export function pushPoint(p: PhysParams): [number, number] {
108 if (p.shape === "ellipse") {
109 const c = Math.SQRT1_2;
110 return [(p.Lx / 2) * c, -(p.Ly / 2) * c];
111 }
112 return [p.Lx / 2, -p.Ly / 2];
113}
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 115/** Stress-free initial shape w0 at every node. */
116export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
117 const w0 = new Float32Array(g.n);
118 if (p.bowShape === "none" || p.bowAmplitude === 0) return w0;
119 for (let j = 0; j < g.ny; j++) {
120 for (let i = 0; i < g.nx; i++) {
121 const X = (2 * nodeX(g, i)) / p.Lx; // -1..1
122 const Y = (2 * nodeY(g, j)) / p.Ly;
123 let f = 0;
124 if (p.bowShape === "width") f = 1 - X * X;
125 else if (p.bowShape === "height") f = 1 - Y * Y;
126 else f = 1 - 0.5 * (X * X + Y * Y);
127 w0[j * g.nx + i] = p.bowAmplitude * f;
128 }
129 }
130 return w0;
131}
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 133export const SPEED_OF_SOUND = 343; // [m/s]
134export const RHO_AIR = 1.2; // [kg/m^3]
135export const HIST_LEN = 256; // ring buffer length for the retarded-time audio sum [samples]
136export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
137export const RB_INTERVAL_HELD = 4; // while the finger constraint is active
139export interface Grid {
140 nx: number;
141 ny: number;
142 hx: number;
143 hy: number;
144 n: number;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 145 /** 1 for nodes that belong to the sheet, 0 otherwise; absent for a full rectangle. */
146 mask?: Uint8Array;
149export interface Material {
150 D: number; // bending stiffness [N m]
151 C: number; // membrane stiffness E h / (1 - nu^2) [N/m]
152 G: number; // membrane shear stiffness E h / (2 (1 + nu)) [N/m]
153 nu: number;
154 rhoh: number; // mass per unit area [kg/m^2]
155}
157export function makeGrid(p: PhysParams): Grid {
158 const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
159 const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 160 const g: Grid = { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
161 return p.shape === "ellipse" ? withEllipseMask(g) : g;
162}
164/**
165 * Restricts a grid to the ellipse inscribed in its bounding box. The edge is a
166 * staircase of whole cells: a cell belongs to the sheet when its center lies
167 * inside the ellipse, and a node when it is a corner of such a cell. Deciding
168 * by cell centers rather than by nodes keeps the covered area unbiased (a
169 * node-based mask loses about half a cell all around the edge, which raises
170 * every frequency by a few percent). Nodes without a cell of their own would
171 * carry no mass; the trailing loop removes any such node.
172 */
173export function withEllipseMask(g: Grid): Grid {
174 const mask = new Uint8Array(g.n);
175 const a = ((g.nx - 1) * g.hx) / 2, b = ((g.ny - 1) * g.hy) / 2;
176 for (let cj = 0; cj < g.ny - 1; cj++) {
177 for (let ci = 0; ci < g.nx - 1; ci++) {
178 const X = (nodeX(g, ci) + g.hx / 2) / a, Y = (nodeY(g, cj) + g.hy / 2) / b;
179 if (X * X + Y * Y <= 1) {
180 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
181 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
182 }
183 }
184 }
185 const gm: Grid = { ...g, mask };
186 for (let changed = true; changed; ) {
187 changed = false;
188 for (let j = 0; j < g.ny; j++) {
189 for (let i = 0; i < g.nx; i++) {
190 if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
191 mask[j * g.nx + i] = 0;
192 changed = true;
193 }
194 }
195 }
196 }
197 return gm;
198}
200/** True when node (i, j) exists and belongs to the sheet. */
201export function insideNode(g: Grid, i: number, j: number): boolean {
202 if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
203 return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
204}
206/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the sheet. */
207export function insideCell(g: Grid, ci: number, cj: number): boolean {
208 return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
209}
211/** Number of inside nodes (the degrees of freedom of the sheet). */
212export function countInside(g: Grid): number {
213 if (!g.mask) return g.n;
214 let c = 0;
215 for (let k = 0; k < g.n; k++) c += g.mask[k];
216 return c;
219export function makeMaterial(p: PhysParams): Material {
220 return {
221 D: (p.E * p.h ** 3) / (12 * (1 - p.nu * p.nu)),
222 C: (p.E * p.h) / (1 - p.nu * p.nu),
223 G: (p.E * p.h) / (2 * (1 + p.nu)),
224 nu: p.nu,
225 rhoh: p.rho * p.h,
226 };
227}
229/** Node coordinates relative to the plate center. */
230export function nodeX(g: Grid, i: number): number {
231 return (i - (g.nx - 1) / 2) * g.hx;
232}
233export function nodeY(g: Grid, j: number): number {
234 return (j - (g.ny - 1) / 2) * g.hy;
235}
238 * Quadrature (and lumped-mass) weight of a node: the fraction of its four
239 * surrounding cells that belong to the sheet. On a full rectangle this is the
240 * trapezoidal rule: 1 in the interior, 1/2 on an edge, 1/4 at a corner. Zero
241 * for nodes outside the sheet.
242 */
243export function nodeWeight(g: Grid, i: number, j: number): number {
244 if (!insideNode(g, i, j)) return 0;
245 let c = 0;
246 if (insideCell(g, i, j)) c++;
247 if (insideCell(g, i - 1, j)) c++;
248 if (insideCell(g, i, j - 1)) c++;
249 if (insideCell(g, i - 1, j - 1)) c++;
250 return c / 4;
251}
253/** Per-node weights as a flat array (0 outside the sheet). */
254export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
255 const w = new Float32Array(g.n);
256 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w[j * g.nx + i] = nodeWeight(g, i, j);
257 return w;
260/** Decay time T60 [s] of a mode with angular frequency omega under the two-parameter loss model. */
261export function t60(m: Material, p: PhysParams, omega: number): number {
262 const kappa = Math.sqrt(m.D / m.rhoh);
263 const k2 = omega / kappa;
264 const sigma = p.sigma0 + p.sigma1 * k2;
265 return (3 * Math.LN10) / sigma; // amplitude factor 1000 = 60 dB
266}