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