1// Physical and numerical parameters of the simulation, and the derived
2// quantities (grid, material constants, geometry of the excitation).
3// All values are SI.
5export type PlateShape = "rectangle" | "ellipse";
6export type BowShape = "none" | "width" | "height" | "dish";
7export type Excitation = "push" | "strike";
9export interface PhysParams {
10 // Sheet
11 shape: PlateShape; // outline of the sheet within its bounding box
12 Lx: number; // width along x [m] (bounding box for an ellipse)
13 Ly: number; // height along y [m]
14 h: number; // thickness [m]
15 E: number; // Young's modulus [Pa]
16 nu: number; // Poisson ratio
17 rho: number; // density [kg/m^3]
18 holeInset: number; // distance of the thread hole from the edges near the top-left corner [m]
19 // Stress-free initial shape
20 bowShape: BowShape;
21 bowAmplitude: number; // sag of the initial shape at the center [m]
22 // Losses
23 sigma0: number; // frequency-independent loss of all motion [1/s]
24 sigma1: number; // loss proportional to the squared wavenumber, transverse motion [m^2/s]
25 sigmaIn: number; // additional frequency-independent loss of the in-plane motion [1/s]
26 // Excitation
27 excitation: Excitation;
28 pushAmount: number; // corner displacement [m]; negative pushes against the bow
29 fingerRadius: number; // radius of the pushed patch [m]
30 rampTime: number; // duration of the push ramp [s]
31 holdMax: number; // maximum settle time before the release [s]
32 releaseTime: number; // time over which the finger force falls to zero [s]
33 settleSigma0: number; // artificial damping during push and settle [1/s]
34 settleKE: number; // non-rigid kinetic energy below which the sheet counts as settled [J]
35 strikeFracX: number; // strike point as a fraction of the width from the left edge
36 strikeFracY: number; // strike point as a fraction of the height from the bottom edge
37 strikeVelocity: number; // stick tip speed at impact [m/s]
38 stickMass: number; // effective mass of the stick [kg]
39 contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
40 contactExponent: number; // Hertz exponent alpha (1.5 for a rounded tip)
41 contactRadius: number; // radius of the contact patch [m]
42 // Listener
43 listenerDist: number; // [m]
44 listenerAzimuth: number; // [deg], rotation about the sheet's vertical axis (0 = on the normal)
45 listenerElevation: number; // [deg]
46 // Numerics and output
47 spacing: number; // requested grid spacing [m]
48 fs: number; // audio sample rate [Hz]
49 fps: number; // movie frame rate, also the batch rate [Hz]
50 duration: number; // recorded duration after the release or the impact [s]
51}
53export const defaultParams: PhysParams = {
54 shape: "rectangle",
55 Lx: 12 * 0.0254,
56 Ly: 18 * 0.0254,
57 h: 0.79e-3,
58 E: 193e9,
59 nu: 0.29,
60 rho: 7900,
61 holeInset: 0.012,
62 bowShape: "width",
63 bowAmplitude: 3e-3,
64 sigma0: 1.15,
65 sigma1: 5e-4,
66 sigmaIn: 20,
67 excitation: "push",
68 pushAmount: 0.01,
69 fingerRadius: 0.01,
70 rampTime: 0.1,
71 holdMax: 0.5,
72 releaseTime: 0.02,
73 settleSigma0: 300,
74 settleKE: 1e-8,
75 strikeFracX: 0.7,
76 strikeFracY: 0.35,
77 strikeVelocity: 3,
78 stickMass: 0.03,
79 contactStiffness: 5e8,
80 contactExponent: 1.5,
81 contactRadius: 5e-3,
82 listenerDist: 1.0,
83 listenerAzimuth: 0,
84 listenerElevation: 0,
85 spacing: 0.008,
86 fs: 48000,
87 fps: 120,
88 duration: 4,
89};
91export const SPEED_OF_SOUND = 343; // [m/s]
92export const RHO_AIR = 1.2; // [kg/m^3]
93export const HIST_LEN = 256; // ring buffer length of the retarded-time listener sum [samples]
94export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
95export const RB_INTERVAL_HELD = 4; // while the finger constraint is active
97// ---------- Grid ----------
99export interface Grid {
100 nx: number;
101 ny: number;
102 hx: number;
103 hy: number;
104 n: number;
105 /** 1 for nodes that belong to the sheet, 0 otherwise; absent for a full rectangle. */
106 mask?: Uint8Array;
107}
109export function makeGrid(p: PhysParams): Grid {
110 const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
111 const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
112 const g: Grid = { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
113 return p.shape === "ellipse" ? withEllipseMask(g) : g;
114}
116/**
117 * Restricts a grid to the ellipse inscribed in its bounding box. The edge is a
118 * staircase of whole cells: a cell belongs to the sheet when its center lies
119 * inside the ellipse, and a node when it is a corner of such a cell. Deciding
120 * by cell centers rather than by nodes keeps the covered area unbiased; a
121 * node-based mask loses about half a cell all around the edge, which raises
122 * every frequency by a few percent. A node with no cell of its own would carry
123 * no mass, and the trailing loop removes any such node.
124 */
125export function withEllipseMask(g: Grid): Grid {
126 const mask = new Uint8Array(g.n);
127 const a = ((g.nx - 1) * g.hx) / 2, b = ((g.ny - 1) * g.hy) / 2;
128 for (let cj = 0; cj < g.ny - 1; cj++) {
129 for (let ci = 0; ci < g.nx - 1; ci++) {
130 const X = (nodeX(g, ci) + g.hx / 2) / a, Y = (nodeY(g, cj) + g.hy / 2) / b;
131 if (X * X + Y * Y <= 1) {
132 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
133 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
134 }
135 }
136 }
137 const gm: Grid = { ...g, mask };
138 for (let changed = true; changed; ) {
139 changed = false;
140 for (let j = 0; j < g.ny; j++) {
141 for (let i = 0; i < g.nx; i++) {
142 if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
143 mask[j * g.nx + i] = 0;
144 changed = true;
145 }
146 }
147 }
148 }
149 return gm;
150}
152/** Node coordinates relative to the center of the bounding box. */
153export function nodeX(g: Grid, i: number): number {
154 return (i - (g.nx - 1) / 2) * g.hx;
155}
156export function nodeY(g: Grid, j: number): number {
157 return (j - (g.ny - 1) / 2) * g.hy;
158}
160/** True when node (i, j) exists and belongs to the sheet. */
161export function insideNode(g: Grid, i: number, j: number): boolean {
162 if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
163 return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
164}
166/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the sheet. */
167export function insideCell(g: Grid, ci: number, cj: number): boolean {
168 return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
169}
171/**
172 * Quadrature and lumped-mass weight of a node: the fraction of its four
173 * surrounding cells that belong to the sheet. On a full rectangle this is the
174 * trapezoidal rule (1 inside, 1/2 on an edge, 1/4 at a corner). Zero outside.
175 */
176export function nodeWeight(g: Grid, i: number, j: number): number {
177 if (!insideNode(g, i, j)) return 0;
178 let c = 0;
179 if (insideCell(g, i, j)) c++;
180 if (insideCell(g, i - 1, j)) c++;
181 if (insideCell(g, i, j - 1)) c++;
182 if (insideCell(g, i - 1, j - 1)) c++;
183 return c / 4;
184}
186/** Per-node weights as a flat array. */
187export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
188 const w = new Float32Array(g.n);
189 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);
190 return w;
191}
193/** Number of nodes in the sheet. */
194export function countInside(g: Grid): number {
195 if (!g.mask) return g.n;
196 let c = 0;
197 for (let k = 0; k < g.n; k++) c += g.mask[k];
198 return c;
199}
201// ---------- Material ----------
203export interface Material {
204 D: number; // bending stiffness E h^3 / (12 (1 - nu^2)) [N m]
205 C: number; // membrane stiffness E h / (1 - nu^2) [N/m]
206 G: number; // membrane shear stiffness E h / (2 (1 + nu)) [N/m]
207 nu: number;
208 rhoh: number; // mass per unit area [kg/m^2]
209}
211export function makeMaterial(p: PhysParams): Material {
212 return {
213 D: (p.E * p.h ** 3) / (12 * (1 - p.nu * p.nu)),
214 C: (p.E * p.h) / (1 - p.nu * p.nu),
215 G: (p.E * p.h) / (2 * (1 + p.nu)),
216 nu: p.nu,
217 rhoh: p.rho * p.h,
218 };
219}
221/**
222 * 60 dB decay time [s] of a bending mode of angular frequency omega under the
223 * loss model sigma0 + sigma1 k^2, using the flat-plate dispersion k^2 = omega / kappa.
224 */
225export function t60(m: Material, p: PhysParams, omega: number): number {
226 const kappa = Math.sqrt(m.D / m.rhoh);
227 const sigma = p.sigma0 + p.sigma1 * (omega / kappa);
228 return (3 * Math.LN10) / sigma;
229}
231/**
232 * Frequency [Hz] at which the three-point second difference lowers a bending
233 * frequency by the relative amount `err` on a grid of spacing hmax. A plane
234 * bending wave along a grid axis with wavenumber k has the discrete frequency
235 * ratio (2 - 2 cos kh) / (kh)^2.
236 */
237export function dispersionFrequency(m: Material, g: Grid, err: number): number {
238 const h = Math.max(g.hx, g.hy);
239 let lo = 0, hi = Math.PI;
240 for (let it = 0; it < 60; it++) {
241 const x = (lo + hi) / 2;
242 if ((2 - 2 * Math.cos(x)) / (x * x) > 1 - err) lo = x; else hi = x;
243 }
244 const k = lo / h;
245 const kappa = Math.sqrt(m.D / m.rhoh);
246 return (kappa * k * k) / (2 * Math.PI);
247}
249// ---------- Geometry of the excitation ----------
251/** Strike point in sheet coordinates. */
252export function strikePoint(p: PhysParams): [number, number] {
253 return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
254}
256/**
257 * Center of the thread hole. Rectangle: inset from the two edges of the
258 * top-left corner. Ellipse: inset along the radius from the rim point in the
259 * top-left direction (parametric angle 135 degrees).
260 */
261export function holePoint(p: PhysParams): [number, number] {
262 if (p.shape === "ellipse") {
263 const c = Math.SQRT1_2;
264 return [-(p.Lx / 2 - p.holeInset) * c, (p.Ly / 2 - p.holeInset) * c];
265 }
266 return [-p.Lx / 2 + p.holeInset, p.Ly / 2 - p.holeInset];
267}
269/** Center of the pushed patch: the bottom-right corner, or the rim point opposite the hole. */
270export function pushPoint(p: PhysParams): [number, number] {
271 if (p.shape === "ellipse") {
272 const c = Math.SQRT1_2;
273 return [(p.Lx / 2) * c, -(p.Ly / 2) * c];
274 }
275 return [p.Lx / 2, -p.Ly / 2];
276}
278/**
279 * Stress-free initial shape w0 at every node. All three shapes are parabolic
280 * with the given sag at the center of the bounding box: "width" and "height"
281 * are cylindrical with zero sag along the two opposite edges, and "dish" has
282 * the same curvature in both directions, so its sag is measured relative to
283 * the corners (the edge midpoints sit at half the sag).
284 */
285export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
286 const w0 = new Float32Array(g.n);
287 if (p.bowShape === "none" || p.bowAmplitude === 0) return w0;
288 for (let j = 0; j < g.ny; j++) {
289 for (let i = 0; i < g.nx; i++) {
290 const X = (2 * nodeX(g, i)) / p.Lx; // -1..1
291 const Y = (2 * nodeY(g, j)) / p.Ly;
292 let f = 0;
293 if (p.bowShape === "width") f = 1 - X * X;
294 else if (p.bowShape === "height") f = 1 - Y * Y;
295 else f = 1 - 0.5 * (X * X + Y * Y);
296 w0[j * g.nx + i] = p.bowAmplitude * f;
297 }
298 }
299 return w0;
300}
302/**
303 * Nodes of the drumstick contact patch: the sheet node nearest the strike
304 * point, plus every sheet node within the contact radius. The contact force
305 * is spread over these nodes in proportion to their mass, and the contact
306 * coordinate is their mass-weighted mean displacement, so the two are
307 * conjugate and the contact has a potential energy.
308 */
309export function contactPatch(p: PhysParams, g: Grid): { nodes: number[]; nearest: number; mass: number; box: [number, number, number, number] } {
310 const m = makeMaterial(p);
311 const [sx, sy] = strikePoint(p);
312 let nearest = -1, best = Infinity;
313 for (let j = 0; j < g.ny; j++) {
314 for (let i = 0; i < g.nx; i++) {
315 if (nodeWeight(g, i, j) === 0) continue;
316 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
317 if (d2 < best) { best = d2; nearest = j * g.nx + i; }
318 }
319 }
320 const nodes: number[] = [];
321 let mass = 0;
322 const box: [number, number, number, number] = [g.nx, -1, g.ny, -1]; // i0, i1, j0, j1 (inclusive)
323 for (let j = 0; j < g.ny; j++) {
324 for (let i = 0; i < g.nx; i++) {
325 const k = j * g.nx + i;
326 const wgt = nodeWeight(g, i, j);
327 if (wgt === 0) continue;
328 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
329 if (d2 <= p.contactRadius ** 2 || k === nearest) {
330 nodes.push(k);
331 mass += m.rhoh * g.hx * g.hy * wgt;
332 box[0] = Math.min(box[0], i); box[1] = Math.max(box[1], i);
333 box[2] = Math.min(box[2], j); box[3] = Math.max(box[3], j);
334 }
335 }
336 }
337 return { nodes, nearest, mass, box };
338}