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, 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 // Numerics and output
20 spacing: number; // requested grid spacing [m]
21 modeMaxFreq: number; // modes up to this frequency are kept in the linear model [Hz]
22 fs: number; // audio sample rate [Hz]
23 fps: number; // movie frame rate, also the batch rate of the full simulation [Hz]
24 duration: number; // recorded duration after the impact [s]
25}
27export const defaultParams: SimParams = {
28 strikeFracX: 0.9,
29 strikeFracY: 0.5,
30 strikeVelocity: 3,
31 stickMass: 0.03,
32 contactStiffness: 1e7,
33 contactExponent: 1.5,
34 contactRadius: 5e-3,
35 listenerDist: 0.5,
36 listenerAngle: 30,
37 listenerAzimuth: 0,
38 spacing: 0.008,
39 modeMaxFreq: 16000,
40 fs: 48000,
41 fps: 120,
42 duration: 4,
43};
45export const SPEED_OF_SOUND = 343; // [m/s]
46export const RHO_AIR = 1.2; // [kg/m^3]
47export const CFL = 0.9; // fraction of the stability limit used for the time step of the full simulation
48export const HIST_MARGIN = 128; // ring buffer margin of the retarded-time listener sum [samples]
49export const STICK_RELEASE = 5e-3; // the stick is withdrawn once the contact has been open this long [s]
51// ---------- Grid ----------
53export interface Grid {
54 nx: number;
55 ny: number;
56 hx: number;
57 hy: number;
58 n: number;
59 /** 1 for nodes that belong to the plate, 0 otherwise; absent for a full rectangle. */
60 mask?: Uint8Array;
61}
63/** True when the point (x, y), relative to the center of the bounding box, lies in the plate. */
64export type InsideFn = (x: number, y: number) => boolean;
66/** The grid covering the bounding box Lx by Ly with about the requested spacing, restricted to the outline when one is given. */
67export function makeGrid(Lx: number, Ly: number, spacing: number, inside?: InsideFn): Grid {
68 const nx = Math.max(6, Math.round(Lx / spacing) + 1);
69 const ny = Math.max(6, Math.round(Ly / spacing) + 1);
70 const g: Grid = { nx, ny, hx: Lx / (nx - 1), hy: Ly / (ny - 1), n: nx * ny };
71 return inside ? withMask(g, inside) : g;
72}
74/**
75 * Restricts a grid to an outline. The edge is a staircase of whole cells: a
76 * cell belongs to the plate when its center lies inside the outline, and a
77 * node when it is a corner of such a cell. Deciding by cell centers keeps
78 * the covered area unbiased. A node left with no cell of its own would carry
79 * no mass, and the trailing loop removes any such node. A hole is handled by
80 * the same rule.
81 */
82export function withMask(g: Grid, inside: InsideFn): Grid {
83 return withCellMask(g, (ci, cj) => inside(nodeX(g, ci) + g.hx / 2, nodeY(g, cj) + g.hy / 2));
84}
86/** The same restriction with the membership of each cell (ci, cj) given directly. */
87export function withCellMask(g: Grid, cellInside: (ci: number, cj: number) => boolean): Grid {
88 const mask = new Uint8Array(g.n);
89 for (let cj = 0; cj < g.ny - 1; cj++) {
90 for (let ci = 0; ci < g.nx - 1; ci++) {
91 if (cellInside(ci, cj)) {
92 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
93 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
94 }
95 }
96 }
97 const gm: Grid = { ...g, mask };
98 for (let changed = true; changed; ) {
99 changed = false;
100 for (let j = 0; j < g.ny; j++) {
101 for (let i = 0; i < g.nx; i++) {
102 if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
103 mask[j * g.nx + i] = 0;
104 changed = true;
105 }
106 }
107 }
108 }
109 return gm;
110}
112/** Node coordinates relative to the center of the bounding box. */
113export function nodeX(g: Grid, i: number): number {
114 return (i - (g.nx - 1) / 2) * g.hx;
115}
116export function nodeY(g: Grid, j: number): number {
117 return (j - (g.ny - 1) / 2) * g.hy;
118}
120/** True when node (i, j) exists and belongs to the plate. */
121export function insideNode(g: Grid, i: number, j: number): boolean {
122 if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
123 return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
124}
126/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the plate. */
127export function insideCell(g: Grid, ci: number, cj: number): boolean {
128 return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
129}
131/**
132 * Quadrature and lumped-mass weight of a node: the fraction of its four
133 * surrounding cells that belong to the plate (1 inside, 1/2 on an edge, 1/4
134 * at a corner of a full rectangle). Zero outside.
135 */
136export function nodeWeight(g: Grid, i: number, j: number): number {
137 if (!insideNode(g, i, j)) return 0;
138 let c = 0;
139 if (insideCell(g, i, j)) c++;
140 if (insideCell(g, i - 1, j)) c++;
141 if (insideCell(g, i, j - 1)) c++;
142 if (insideCell(g, i - 1, j - 1)) c++;
143 return c / 4;
144}
146/** Per-node weights as a flat array. */
147export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
148 const w = new Float32Array(g.n);
149 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);
150 return w;
151}
153/** Per-cell membership (1 or 0) as a flat array of (nx - 1) * (ny - 1). */
154export function cellMask(g: Grid): Float32Array<ArrayBuffer> {
155 const c = new Float32Array((g.nx - 1) * (g.ny - 1));
156 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;
157 return c;
158}
160/** Number of nodes in the plate. */
161export function countInside(g: Grid): number {
162 if (!g.mask) return g.n;
163 let c = 0;
164 for (let k = 0; k < g.n; k++) c += g.mask[k];
165 return c;
166}
168// ---------- Material, losses, plate ----------
170/** Thermal constants of the metal, from which the thermoelastic damping follows (loss.ts). */
171export interface Thermal {
172 alpha: number; // linear thermal expansion coefficient [1/K]
173 cp: number; // specific heat [J/(kg K)]
174 k: number; // thermal conductivity [W/(m K)]
175}
177export interface Material {
178 E: number; // Young's modulus [Pa]
179 nu: number; // Poisson ratio
180 rho: number; // density [kg/m^3]
181 /** Loss factor of the metal's internal friction (hysteretic, frequency independent); 0 for none. */
182 eta: number;
183 /** Thermal constants, or null when the thermoelastic damping is left out. */
184 thermal: Thermal | null;
185}
187/** Temperature of the plate [K], which enters the thermoelastic damping. */
188export const ROOM_TEMPERATURE = 293;
190/** The plate as evaluated on the grid. */
191export interface Plate {
192 name: string;
193 Lx: number;
194 Ly: number;
195 grid: Grid;
196 mat: Material;
197 /** Thickness at the nodes [m]; zero outside the plate. */
198 h: Float64Array;
199 /** Stress-free height at the nodes [m], mean plane removed; zero outside the plate. */
200 w0: Float64Array;
201 /** Temperature [K], for the thermoelastic damping. */
202 temperature: number;
203}
205export function bendingStiffness(m: Material, h: number): number {
206 return (m.E * h * h * h) / (12 * (1 - m.nu * m.nu));
207}
208export function membraneStiffness(m: Material, h: number): number {
209 return (m.E * h) / (1 - m.nu * m.nu);
210}
211export function shearStiffness(m: Material, h: number): number {
212 return (m.E * h) / (2 * (1 + m.nu));
213}
214/** In-plane wave speed, which sets the time step of the full simulation. */
215export function inPlaneSpeed(m: Material): number {
216 return Math.sqrt(m.E / (m.rho * (1 - m.nu * m.nu)));
217}
218/** Bending-wave constant kappa = sqrt(D / rho h), so that omega = kappa k^2. */
219export function bendingConstant(m: Material, h: number): number {
220 return Math.sqrt(bendingStiffness(m, h) / (m.rho * h));
221}
223/** Thickness of a cell: the mean of its four corners. */
224export function cellThickness(pl: Plate, ci: number, cj: number): number {
225 const nx = pl.grid.nx;
226 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]);
227}
229/** Lumped mass of a node: rho h dA times its weight. */
230export function nodeMass(pl: Plate, i: number, j: number): number {
231 const g = pl.grid;
232 return pl.mat.rho * pl.h[j * g.nx + i] * g.hx * g.hy * nodeWeight(g, i, j);
233}
235/** Per-node masses as a flat array. */
236export function nodeMasses(pl: Plate): Float64Array {
237 const g = pl.grid;
238 const m = new Float64Array(g.n);
239 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);
240 return m;
241}
243/** Total mass, the mass-weighted center, and the moments of inertia about the in-plane axes through it. */
244export function massProperties(pl: Plate): { mass: number; cx: number; cy: number; Ix: number; Iy: number; area: number } {
245 const g = pl.grid;
246 let mass = 0, sx = 0, sy = 0, area = 0;
247 for (let j = 0; j < g.ny; j++) {
248 for (let i = 0; i < g.nx; i++) {
249 const m = nodeMass(pl, i, j);
250 mass += m; sx += m * nodeX(g, i); sy += m * nodeY(g, j);
251 area += g.hx * g.hy * nodeWeight(g, i, j);
252 }
253 }
254 const cx = mass > 0 ? sx / mass : 0, cy = mass > 0 ? sy / mass : 0;
255 let Ix = 0, Iy = 0;
256 for (let j = 0; j < g.ny; j++) {
257 for (let i = 0; i < g.nx; i++) {
258 const m = nodeMass(pl, i, j);
259 Ix += m * (nodeY(g, j) - cy) ** 2;
260 Iy += m * (nodeX(g, i) - cx) ** 2;
261 }
262 }
263 return { mass, cx, cy, Ix, Iy, area };
264}
266/** Mean thickness over the plate (area-weighted). */
267export function meanThickness(pl: Plate): number {
268 const g = pl.grid;
269 let s = 0, a = 0;
270 for (let j = 0; j < g.ny; j++) {
271 for (let i = 0; i < g.nx; i++) {
272 const w = nodeWeight(g, i, j);
273 s += w * pl.h[j * g.nx + i]; a += w;
274 }
275 }
276 return a > 0 ? s / a : 0;
277}
279/** Smallest and largest thickness over the plate. */
280export function thicknessRange(pl: Plate): { min: number; max: number } {
281 const g = pl.grid;
282 let lo = Infinity, hi = 0;
283 for (let j = 0; j < g.ny; j++) {
284 for (let i = 0; i < g.nx; i++) {
285 if (!insideNode(g, i, j)) continue;
286 const h = pl.h[j * g.nx + i];
287 lo = Math.min(lo, h); hi = Math.max(hi, h);
288 }
289 }
290 return { min: isFinite(lo) ? lo : 0, max: hi };
291}
294// ---------- Patches ----------
296/** A set of plate nodes with the quantities the strike needs. */
297export interface Patch {
298 nodes: number[];
299 nearest: number;
300 mass: number;
301 cx: number;
302 cy: number;
303}
305/** The plate nodes within `radius` of (x, y), always including the nearest one. */
306export function patchAround(pl: Plate, x: number, y: number, radius: number): Patch {
307 const g = pl.grid;
308 const nodes: number[] = [];
309 let nearest = -1, best = Infinity;
310 for (let j = 0; j < g.ny; j++) {
311 for (let i = 0; i < g.nx; i++) {
312 if (nodeMass(pl, i, j) === 0) continue;
313 const d = Math.hypot(nodeX(g, i) - x, nodeY(g, j) - y);
314 if (d < best) { best = d; nearest = j * g.nx + i; }
315 if (d <= radius) nodes.push(j * g.nx + i);
316 }
317 }
318 if (nearest >= 0 && !nodes.includes(nearest)) nodes.push(nearest);
319 let mass = 0, cx = 0, cy = 0;
320 for (const k of nodes) {
321 const m = nodeMass(pl, k % g.nx, Math.floor(k / g.nx));
322 mass += m; cx += m * nodeX(g, k % g.nx); cy += m * nodeY(g, Math.floor(k / g.nx));
323 }
324 return { nodes, nearest, mass, cx: mass > 0 ? cx / mass : 0, cy: mass > 0 ? cy / mass : 0 };
325}
327/** The strike point in plate coordinates from the fractions of the bounding box. */
328export function strikePoint(pl: Plate, p: SimParams): { x: number; y: number } {
329 return { x: (p.strikeFracX - 0.5) * pl.Lx, y: (p.strikeFracY - 0.5) * pl.Ly };
330}
332export function strikePatch(pl: Plate, p: SimParams): Patch {
333 const { x, y } = strikePoint(pl, p);
334 return patchAround(pl, x, y, p.contactRadius);
335}
337// ---------- Listener ----------
339/** Listener position relative to the center of the plate, z along the plate's normal. */
340export function listenerPosition(p: SimParams): { x: number; y: number; z: number } {
341 const th = (p.listenerAngle * Math.PI) / 180, ph = (p.listenerAzimuth * Math.PI) / 180;
342 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) };
343}
345/**
346 * Per-node terms of the Rayleigh integral for the listener: the weight
347 * rho_air dA / (2 pi r) and the delay r / c in seconds, for every node (zero
348 * outside the plate). The profile is included in the node positions.
349 */
350export function listenerTerms(pl: Plate, p: SimParams): { weight: Float64Array; delay: Float64Array } {
351 const g = pl.grid;
352 const L = listenerPosition(p);
353 const weight = new Float64Array(g.n), delay = new Float64Array(g.n);
354 for (let j = 0; j < g.ny; j++) {
355 for (let i = 0; i < g.nx; i++) {
356 const wgt = nodeWeight(g, i, j);
357 if (wgt === 0) continue;
358 const k = j * g.nx + i;
359 const r = Math.hypot(L.x - nodeX(g, i), L.y - nodeY(g, j), L.z - pl.w0[k]);
360 weight[k] = (RHO_AIR * g.hx * g.hy * wgt) / (2 * Math.PI * Math.max(r, 1e-3));
361 delay[k] = r / SPEED_OF_SOUND;
362 }
363 }
364 return { weight, delay };
365}
367// ---------- Grid dispersion ----------
369/**
370 * Frequencies at which grid dispersion lowers a bending partial by the given
371 * fraction, for a plate of thickness h: a plane wave along a grid axis has the
372 * discrete frequency ratio (2 - 2 cos kh) / (kh)^2.
373 */
374export function dispersionFrequency(m: Material, h: number, spacing: number, fraction: number): number {
375 let lo = 0, hi = Math.PI;
376 for (let it = 0; it < 60; it++) {
377 const mid = (lo + hi) / 2;
378 const ratio = (2 - 2 * Math.cos(mid)) / (mid * mid);
379 if (1 - ratio < fraction) lo = mid; else hi = mid;
380 }
381 const k = lo / spacing;
382 return (bendingConstant(m, h) * k * k) / (2 * Math.PI);
383}