1// Simulation parameters (strike, listener, numerics), the grid with its mask,
2// the plate as evaluated on the grid (material, thickness field, stress-free
3// shape, stand, loss constants), and the derived quantities both solvers need. All
4// values are SI. The design itself is produced by the script in design.ts.
6export interface SimParams {
7 // Strike
8 strikeFracX: number; // strike point as a fraction of the bounding-box width from the left edge
9 strikeFracY: number; // strike point as a fraction of the bounding-box height from the bottom edge
10 strikeVelocity: number; // stick tip speed at impact [m/s]
11 stickMass: number; // effective mass of the stick [kg]
12 contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
13 contactExponent: number; // Hertz exponent alpha (1.5 for a rounded tip)
14 contactRadius: number; // radius of the contact patch [m]
15 // Listener
16 listenerDist: number; // [m]
17 listenerAngle: number; // [deg] angle from the plate's normal (0 = straight above the plate)
18 listenerAzimuth: number; // [deg] rotation about the normal (0 = toward +x)
19 // Stand (a parameter of the setup, not of the design script)
20 standOn: boolean; // false for a free plate
21 standX: number; // position of the rod relative to the center of the bounding box [m]
22 standY: number;
23 padInnerRadius: number; // the felt annulus [m]
24 padOuterRadius: number;
25 padThickness: number; // [m]
26 padE: number; // Young's modulus of the felt in compression [Pa]
27 padLossFactor: number;
28 topPad: boolean; // a second pad above the plate under the wing nut
29 preload: number; // compression of the top pad at rest [m]
30 // Numerics and output
31 spacing: number; // requested grid spacing [m]
32 modeMaxFreq: number; // modes up to this frequency are kept in the linear model [Hz]
33 padLiftoff: boolean; // full simulation: the pads push only while compressed
34 fs: number; // audio sample rate [Hz]
35 fps: number; // movie frame rate, also the batch rate of the full simulation [Hz]
36 duration: number; // recorded duration after the impact [s]
37}
39export const defaultParams: SimParams = {
40 strikeFracX: 0.9,
41 strikeFracY: 0.5,
42 strikeVelocity: 3,
43 stickMass: 0.03,
44 contactStiffness: 1e7,
45 contactExponent: 1.5,
46 contactRadius: 5e-3,
47 listenerDist: 0.5,
48 listenerAngle: 30,
49 listenerAzimuth: 0,
50 standOn: true,
51 standX: 0,
52 standY: 0,
53 padInnerRadius: 7e-3,
54 padOuterRadius: 20e-3,
55 padThickness: 12e-3,
56 padE: 0.5e6,
57 padLossFactor: 0.3,
58 topPad: true,
59 preload: 1e-3,
60 spacing: 0.008,
61 modeMaxFreq: 16000,
62 padLiftoff: true,
63 fs: 48000,
64 fps: 120,
65 duration: 4,
66};
68export const SPEED_OF_SOUND = 343; // [m/s]
69export const RHO_AIR = 1.2; // [kg/m^3]
70export const GRAVITY = 9.81; // [m/s^2]
71export const CFL = 0.9; // fraction of the stability limit used for the time step of the full simulation
72export const HIST_MARGIN = 128; // ring buffer margin of the retarded-time listener sum [samples]
73export const STICK_RELEASE = 5e-3; // the stick is withdrawn once the contact has been open this long [s]
74export const PAD_POISSON = 0.3; // Poisson ratio of the felt, which sets its shear stiffness from E
76// ---------- Grid ----------
78export interface Grid {
79 nx: number;
80 ny: number;
81 hx: number;
82 hy: number;
83 n: number;
84 /** 1 for nodes that belong to the plate, 0 otherwise; absent for a full rectangle. */
85 mask?: Uint8Array;
86}
88/** True when the point (x, y), relative to the center of the bounding box, lies in the plate. */
89export type InsideFn = (x: number, y: number) => boolean;
91/** The grid covering the bounding box Lx by Ly with about the requested spacing, restricted to the outline when one is given. */
92export function makeGrid(Lx: number, Ly: number, spacing: number, inside?: InsideFn): Grid {
93 const nx = Math.max(6, Math.round(Lx / spacing) + 1);
94 const ny = Math.max(6, Math.round(Ly / spacing) + 1);
95 const g: Grid = { nx, ny, hx: Lx / (nx - 1), hy: Ly / (ny - 1), n: nx * ny };
96 return inside ? withMask(g, inside) : g;
97}
99/**
100 * Restricts a grid to an outline. The edge is a staircase of whole cells: a
101 * cell belongs to the plate when its center lies inside the outline, and a
102 * node when it is a corner of such a cell. Deciding by cell centers keeps
103 * the covered area unbiased. A node left with no cell of its own would carry
104 * no mass, and the trailing loop removes any such node. A hole is handled by
105 * the same rule.
106 */
107export function withMask(g: Grid, inside: InsideFn): Grid {
108 return withCellMask(g, (ci, cj) => inside(nodeX(g, ci) + g.hx / 2, nodeY(g, cj) + g.hy / 2));
109}
111/** The same restriction with the membership of each cell (ci, cj) given directly. */
112export function withCellMask(g: Grid, cellInside: (ci: number, cj: number) => boolean): Grid {
113 const mask = new Uint8Array(g.n);
114 for (let cj = 0; cj < g.ny - 1; cj++) {
115 for (let ci = 0; ci < g.nx - 1; ci++) {
116 if (cellInside(ci, cj)) {
117 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
118 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
119 }
120 }
121 }
122 const gm: Grid = { ...g, mask };
123 for (let changed = true; changed; ) {
124 changed = false;
125 for (let j = 0; j < g.ny; j++) {
126 for (let i = 0; i < g.nx; i++) {
127 if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
128 mask[j * g.nx + i] = 0;
129 changed = true;
130 }
131 }
132 }
133 }
134 return gm;
135}
137/** Node coordinates relative to the center of the bounding box. */
138export function nodeX(g: Grid, i: number): number {
139 return (i - (g.nx - 1) / 2) * g.hx;
140}
141export function nodeY(g: Grid, j: number): number {
142 return (j - (g.ny - 1) / 2) * g.hy;
143}
145/** True when node (i, j) exists and belongs to the plate. */
146export function insideNode(g: Grid, i: number, j: number): boolean {
147 if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
148 return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
149}
151/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the plate. */
152export function insideCell(g: Grid, ci: number, cj: number): boolean {
153 return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
154}
156/**
157 * Quadrature and lumped-mass weight of a node: the fraction of its four
158 * surrounding cells that belong to the plate (1 inside, 1/2 on an edge, 1/4
159 * at a corner of a full rectangle). Zero outside.
160 */
161export function nodeWeight(g: Grid, i: number, j: number): number {
162 if (!insideNode(g, i, j)) return 0;
163 let c = 0;
164 if (insideCell(g, i, j)) c++;
165 if (insideCell(g, i - 1, j)) c++;
166 if (insideCell(g, i, j - 1)) c++;
167 if (insideCell(g, i - 1, j - 1)) c++;
168 return c / 4;
169}
171/** Per-node weights as a flat array. */
172export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
173 const w = new Float32Array(g.n);
174 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);
175 return w;
176}
178/** Per-cell membership (1 or 0) as a flat array of (nx - 1) * (ny - 1). */
179export function cellMask(g: Grid): Float32Array<ArrayBuffer> {
180 const c = new Float32Array((g.nx - 1) * (g.ny - 1));
181 for (let cj = 0; cj < g.ny - 1; cj++) for (let ci = 0; ci < g.nx - 1; ci++) c[cj * (g.nx - 1) + ci] = insideCell(g, ci, cj) ? 1 : 0;
182 return c;
183}
185/** Number of nodes in the plate. */
186export function countInside(g: Grid): number {
187 if (!g.mask) return g.n;
188 let c = 0;
189 for (let k = 0; k < g.n; k++) c += g.mask[k];
190 return c;
191}
193// ---------- Material, stand, losses, plate ----------
195/** Thermal constants of the metal, from which the thermoelastic damping follows (loss.ts). */
196export interface Thermal {
197 alpha: number; // linear thermal expansion coefficient [1/K]
198 cp: number; // specific heat [J/(kg K)]
199 k: number; // thermal conductivity [W/(m K)]
200}
202export interface Material {
203 E: number; // Young's modulus [Pa]
204 nu: number; // Poisson ratio
205 rho: number; // density [kg/m^3]
206 /** Loss factor of the metal's internal friction (hysteretic, frequency independent); 0 for none. */
207 eta: number;
208 /** Thermal constants, or null when the thermoelastic damping is left out. */
209 thermal: Thermal | null;
210}
212/** Temperature of the plate [K], which enters the thermoelastic damping. */
213export const ROOM_TEMPERATURE = 293;
215/** A felt pad of the stand: an annulus of felt between the plate and the stand (or the wing nut). */
216export interface Pad {
217 innerRadius: number; // [m]
218 outerRadius: number; // [m]
219 thickness: number; // [m]
220 E: number; // Young's modulus of the felt in compression [Pa]
221 lossFactor: number; // loss factor eta of the felt
222}
224/**
225 * The stand: a rod through the center hole at (x, y) with a felt pad below
226 * the plate and, when topPad is set, the same pad above it under the wing
227 * nut, which compresses the top pad by `preload` at rest.
228 */
229export interface Stand {
230 x: number;
231 y: number;
232 pad: Pad;
233 topPad: boolean;
234 preload: number; // compression of the top pad at rest [m]
235}
237/** The stand described by the parameters, or null for a free plate. */
238export function standFromParams(p: SimParams): Stand | null {
239 if (!p.standOn) return null;
240 return {
241 x: p.standX,
242 y: p.standY,
243 pad: { innerRadius: p.padInnerRadius, outerRadius: Math.max(p.padOuterRadius, p.padInnerRadius + 1e-4), thickness: p.padThickness, E: p.padE, lossFactor: p.padLossFactor },
244 topPad: p.topPad,
245 preload: p.preload,
246 };
247}
249/** The plate as evaluated on the grid. */
250export interface Plate {
251 name: string;
252 Lx: number;
253 Ly: number;
254 grid: Grid;
255 mat: Material;
256 /** Thickness at the nodes [m]; zero outside the plate. */
257 h: Float64Array;
258 /** Stress-free height at the nodes [m], mean plane removed; zero outside the plate. */
259 w0: Float64Array;
260 stand: Stand | null;
261 /** Temperature [K], for the thermoelastic damping. */
262 temperature: number;
263}
265export function bendingStiffness(m: Material, h: number): number {
266 return (m.E * h * h * h) / (12 * (1 - m.nu * m.nu));
267}
268export function membraneStiffness(m: Material, h: number): number {
269 return (m.E * h) / (1 - m.nu * m.nu);
270}
271export function shearStiffness(m: Material, h: number): number {
272 return (m.E * h) / (2 * (1 + m.nu));
273}
274/** In-plane wave speed, which sets the time step of the full simulation. */
275export function inPlaneSpeed(m: Material): number {
276 return Math.sqrt(m.E / (m.rho * (1 - m.nu * m.nu)));
277}
278/** Bending-wave constant kappa = sqrt(D / rho h), so that omega = kappa k^2. */
279export function bendingConstant(m: Material, h: number): number {
280 return Math.sqrt(bendingStiffness(m, h) / (m.rho * h));
281}
283/** Thickness of a cell: the mean of its four corners. */
284export function cellThickness(pl: Plate, ci: number, cj: number): number {
285 const nx = pl.grid.nx;
286 return 0.25 * (pl.h[cj * nx + ci] + pl.h[cj * nx + ci + 1] + pl.h[(cj + 1) * nx + ci] + pl.h[(cj + 1) * nx + ci + 1]);
287}
289/** Lumped mass of a node: rho h dA times its weight. */
290export function nodeMass(pl: Plate, i: number, j: number): number {
291 const g = pl.grid;
292 return pl.mat.rho * pl.h[j * g.nx + i] * g.hx * g.hy * nodeWeight(g, i, j);
293}
295/** Per-node masses as a flat array. */
296export function nodeMasses(pl: Plate): Float64Array {
297 const g = pl.grid;
298 const m = new Float64Array(g.n);
299 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) m[j * g.nx + i] = nodeMass(pl, i, j);
300 return m;
301}
303/** Total mass, the mass-weighted center, and the moments of inertia about the in-plane axes through it. */
304export function massProperties(pl: Plate): { mass: number; cx: number; cy: number; Ix: number; Iy: number; area: number } {
305 const g = pl.grid;
306 let mass = 0, sx = 0, sy = 0, area = 0;
307 for (let j = 0; j < g.ny; j++) {
308 for (let i = 0; i < g.nx; i++) {
309 const m = nodeMass(pl, i, j);
310 mass += m; sx += m * nodeX(g, i); sy += m * nodeY(g, j);
311 area += g.hx * g.hy * nodeWeight(g, i, j);
312 }
313 }
314 const cx = mass > 0 ? sx / mass : 0, cy = mass > 0 ? sy / mass : 0;
315 let Ix = 0, Iy = 0;
316 for (let j = 0; j < g.ny; j++) {
317 for (let i = 0; i < g.nx; i++) {
318 const m = nodeMass(pl, i, j);
319 Ix += m * (nodeY(g, j) - cy) ** 2;
320 Iy += m * (nodeX(g, i) - cx) ** 2;
321 }
322 }
323 return { mass, cx, cy, Ix, Iy, area };
324}
326/** Mean thickness over the plate (area-weighted). */
327export function meanThickness(pl: Plate): number {
328 const g = pl.grid;
329 let s = 0, a = 0;
330 for (let j = 0; j < g.ny; j++) {
331 for (let i = 0; i < g.nx; i++) {
332 const w = nodeWeight(g, i, j);
333 s += w * pl.h[j * g.nx + i]; a += w;
334 }
335 }
336 return a > 0 ? s / a : 0;
337}
339/** Smallest and largest thickness over the plate. */
340export function thicknessRange(pl: Plate): { min: number; max: number } {
341 const g = pl.grid;
342 let lo = Infinity, hi = 0;
343 for (let j = 0; j < g.ny; j++) {
344 for (let i = 0; i < g.nx; i++) {
345 if (!insideNode(g, i, j)) continue;
346 const h = pl.h[j * g.nx + i];
347 lo = Math.min(lo, h); hi = Math.max(hi, h);
348 }
349 }
350 return { min: isFinite(lo) ? lo : 0, max: hi };
351}
353// ---------- Stand as seen by the solvers ----------
355/**
356 * The pads as per-node spring rates. `kw[k]` is the transverse stiffness of
357 * the pads at node k (both pads summed), `ku[k]` the in-plane (shear)
358 * stiffness, both in N/m; `kwBottom` is the bottom pad alone, which the full
359 * simulation needs for the lift-off rule. `eta` is the loss factor, and
360 * `preloadTop`, `preloadBottom` are the pad compressions at rest.
361 */
362export interface StandRates {
363 kw: Float64Array;
364 ku: Float64Array;
365 kwBottom: Float64Array;
366 nodes: number[];
367 eta: number;
368 preloadTop: number;
369 preloadBottom: number;
370 /** Total transverse stiffness of the pads [N/m]. */
371 totalKw: number;
372 /** Bounce frequency of the plate on the pads [Hz]. */
373 bounceHz: number;
374}
376/**
377 * Per-node pad stiffness of the stand, or null for a free plate. The pads
378 * act on the plate nodes within the annulus; because a coarse grid covers
379 * the annulus with a handful of nodes, the per-node rates are scaled so
380 * that their sum equals the stiffness of the exact annulus area.
381 */
382export function standRates(pl: Plate): StandRates | null {
383 const st = pl.stand;
384 if (!st) return null;
385 const g = pl.grid;
386 const kw = new Float64Array(g.n), ku = new Float64Array(g.n), kwBottom = new Float64Array(g.n);
387 const nodes: number[] = [];
388 const pad = st.pad;
389 const kArea = pad.E / pad.thickness; // N/m^3
390 const gArea = pad.E / (2 * (1 + PAD_POISSON)) / pad.thickness;
391 const nPads = st.topPad ? 2 : 1;
392 const annulus = Math.PI * (pad.outerRadius ** 2 - pad.innerRadius ** 2);
393 let nodeArea = 0;
394 for (let j = 0; j < g.ny; j++) {
395 for (let i = 0; i < g.nx; i++) {
396 const wgt = nodeWeight(g, i, j);
397 if (wgt === 0) continue;
398 const r = Math.hypot(nodeX(g, i) - st.x, nodeY(g, j) - st.y);
399 if (r < pad.innerRadius || r > pad.outerRadius) continue;
400 nodes.push(j * g.nx + i);
401 nodeArea += g.hx * g.hy * wgt;
402 }
403 }
404 const empty = { kw, ku, kwBottom, nodes, eta: pad.lossFactor, preloadTop: 0, preloadBottom: 0, totalKw: 0, bounceHz: 0 };
405 if (nodes.length === 0 || nodeArea === 0) return empty;
406 const scale = annulus / nodeArea;
407 let totalKw = 0;
408 for (const k of nodes) {
409 const dA = scale * g.hx * g.hy * nodeWeight(g, k % g.nx, Math.floor(k / g.nx));
410 kwBottom[k] = kArea * dA;
411 kw[k] = nPads * kArea * dA;
412 ku[k] = nPads * gArea * dA;
413 totalKw += kw[k];
414 }
415 const { mass } = massProperties(pl);
416 const preloadTop = st.topPad ? st.preload : 0;
417 // The bottom pad carries the top pad's preload plus the weight of the plate.
418 const preloadBottom = preloadTop + (mass * GRAVITY) / (kArea * annulus);
419 return { kw, ku, kwBottom, nodes, eta: pad.lossFactor, preloadTop, preloadBottom, totalKw, bounceHz: Math.sqrt(totalKw / mass) / (2 * Math.PI) };
420}
422// ---------- Patches ----------
424/** A set of plate nodes with the quantities the strike needs. */
425export interface Patch {
426 nodes: number[];
427 nearest: number;
428 mass: number;
429 cx: number;
430 cy: number;
431}
433/** The plate nodes within `radius` of (x, y), always including the nearest one. */
434export function patchAround(pl: Plate, x: number, y: number, radius: number): Patch {
435 const g = pl.grid;
436 const nodes: number[] = [];
437 let nearest = -1, best = Infinity;
438 for (let j = 0; j < g.ny; j++) {
439 for (let i = 0; i < g.nx; i++) {
440 if (nodeMass(pl, i, j) === 0) continue;
441 const d = Math.hypot(nodeX(g, i) - x, nodeY(g, j) - y);
442 if (d < best) { best = d; nearest = j * g.nx + i; }
443 if (d <= radius) nodes.push(j * g.nx + i);
444 }
445 }
446 if (nearest >= 0 && !nodes.includes(nearest)) nodes.push(nearest);
447 let mass = 0, cx = 0, cy = 0;
448 for (const k of nodes) {
449 const m = nodeMass(pl, k % g.nx, Math.floor(k / g.nx));
450 mass += m; cx += m * nodeX(g, k % g.nx); cy += m * nodeY(g, Math.floor(k / g.nx));
451 }
452 return { nodes, nearest, mass, cx: mass > 0 ? cx / mass : 0, cy: mass > 0 ? cy / mass : 0 };
453}
455/** The strike point in plate coordinates from the fractions of the bounding box. */
456export function strikePoint(pl: Plate, p: SimParams): { x: number; y: number } {
457 return { x: (p.strikeFracX - 0.5) * pl.Lx, y: (p.strikeFracY - 0.5) * pl.Ly };
458}
460export function strikePatch(pl: Plate, p: SimParams): Patch {
461 const { x, y } = strikePoint(pl, p);
462 return patchAround(pl, x, y, p.contactRadius);
463}
465// ---------- Listener ----------
467/** Listener position relative to the center of the plate, z along the plate's normal. */
468export function listenerPosition(p: SimParams): { x: number; y: number; z: number } {
469 const th = (p.listenerAngle * Math.PI) / 180, ph = (p.listenerAzimuth * Math.PI) / 180;
470 return { x: p.listenerDist * Math.sin(th) * Math.cos(ph), y: p.listenerDist * Math.sin(th) * Math.sin(ph), z: p.listenerDist * Math.cos(th) };
471}
473/**
474 * Per-node terms of the Rayleigh integral for the listener: the weight
475 * rho_air dA / (2 pi r) and the delay r / c in seconds, for every node (zero
476 * outside the plate). The profile is included in the node positions.
477 */
478export function listenerTerms(pl: Plate, p: SimParams): { weight: Float64Array; delay: Float64Array } {
479 const g = pl.grid;
480 const L = listenerPosition(p);
481 const weight = new Float64Array(g.n), delay = new Float64Array(g.n);
482 for (let j = 0; j < g.ny; j++) {
483 for (let i = 0; i < g.nx; i++) {
484 const wgt = nodeWeight(g, i, j);
485 if (wgt === 0) continue;
486 const k = j * g.nx + i;
487 const r = Math.hypot(L.x - nodeX(g, i), L.y - nodeY(g, j), L.z - pl.w0[k]);
488 weight[k] = (RHO_AIR * g.hx * g.hy * wgt) / (2 * Math.PI * Math.max(r, 1e-3));
489 delay[k] = r / SPEED_OF_SOUND;
490 }
491 }
492 return { weight, delay };
493}
495// ---------- Grid dispersion ----------
497/**
498 * Frequencies at which grid dispersion lowers a bending partial by the given
499 * fraction, for a plate of thickness h: a plane wave along a grid axis has the
500 * discrete frequency ratio (2 - 2 cos kh) / (kh)^2.
501 */
502export function dispersionFrequency(m: Material, h: number, spacing: number, fraction: number): number {
503 let lo = 0, hi = Math.PI;
504 for (let it = 0; it < 60; it++) {
505 const mid = (lo + hi) / 2;
506 const ratio = (2 - 2 * Math.cos(mid)) / (mid * mid);
507 if (1 - ratio < fraction) lo = mid; else hi = mid;
508 }
509 const k = lo / spacing;
510 return (bendingConstant(m, h) * k * k) / (2 * Math.PI);
511}