d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 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.
8cb683aInitial shape as a list of features, with random configurationsJeremy Magland 5import type { ShapeFeature } from "./shape.ts";
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 7export type PlateShape = "rectangle" | "ellipse";
8export type Excitation = "push" | "strike";
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 9export type HitSequence = "single" | "steady" | "crescendo" | "decrescendo";
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 10
11export interface PhysParams {
12 // Sheet
13 shape: PlateShape; // outline of the sheet within its bounding box
14 Lx: number; // width along x [m] (bounding box for an ellipse)
15 Ly: number; // height along y [m]
16 h: number; // thickness [m]
17 E: number; // Young's modulus [Pa]
18 nu: number; // Poisson ratio
19 rho: number; // density [kg/m^3]
20 holeInset: number; // distance of the thread hole from the edges near the top-left corner [m]
8cb683aInitial shape as a list of features, with random configurationsJeremy Magland 21 // Stress-free initial shape: the sum of these features (see shape.ts)
22 features: ShapeFeature[];
24 sigma0: number; // frequency-independent loss of all motion [1/s]
25 sigma1: number; // loss proportional to the squared wavenumber, transverse motion [m^2/s]
26 sigmaIn: number; // additional frequency-independent loss of the in-plane motion [1/s]
27 // Excitation
28 excitation: Excitation;
29 pushAmount: number; // corner displacement [m]; negative pushes against the bow
30 fingerRadius: number; // radius of the pushed patch [m]
31 rampTime: number; // duration of the push ramp [s]
32 holdMax: number; // maximum settle time before the release [s]
33 releaseTime: number; // time over which the finger force falls to zero [s]
34 settleSigma0: number; // artificial damping during push and settle [1/s]
35 settleKE: number; // non-rigid kinetic energy below which the sheet counts as settled [J]
36 strikeFracX: number; // strike point as a fraction of the width from the left edge
37 strikeFracY: number; // strike point as a fraction of the height from the bottom edge
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 38 strikeVelocity: number; // stick tip speed at impact [m/s]; the loudest hit of a sequence
39 hitSequence: HitSequence; // one hit, or a series of hits every hitInterval seconds
40 hitCount: number; // number of hits in a series
41 hitInterval: number; // time between the hits of a series [s]
42 strikeVelocityMin: number; // softest stick speed of a crescendo or decrescendo [m/s]
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 43 stickMass: number; // effective mass of the stick [kg]
44 contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
45 contactExponent: number; // Hertz exponent alpha (1.5 for a rounded tip)
46 contactRadius: number; // radius of the contact patch [m]
47 // Listener
48 listenerDist: number; // [m]
49 listenerAzimuth: number; // [deg], rotation about the sheet's vertical axis (0 = on the normal)
50 listenerElevation: number; // [deg]
51 // Numerics and output
52 spacing: number; // requested grid spacing [m]
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 53 inPlaneMass: number; // factor by which the in-plane inertia is scaled; 1 is physical, larger lengthens the time step
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 54 rigidProjection: boolean; // remove the six rigid-body modes every RB_INTERVAL steps while ringing
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 55 fs: number; // audio sample rate [Hz]
56 fps: number; // movie frame rate, also the batch rate [Hz]
57 duration: number; // recorded duration after the release or the impact [s]
58}
60export const defaultParams: PhysParams = {
61 shape: "rectangle",
62 Lx: 12 * 0.0254,
63 Ly: 18 * 0.0254,
64 h: 0.79e-3,
67 rho: 7900,
68 holeInset: 0.012,
8cb683aInitial shape as a list of features, with random configurationsJeremy Magland 69 features: [{ type: "bow", axis: "x", sag: 3e-3, offset: 0 }],
71 sigma1: 5e-4,
72 sigmaIn: 20,
75 fingerRadius: 0.01,
76 rampTime: 0.1,
77 holdMax: 0.5,
78 releaseTime: 0.02,
79 settleSigma0: 300,
80 settleKE: 1e-8,
81 strikeFracX: 0.7,
82 strikeFracY: 0.35,
83 strikeVelocity: 3,
85 hitCount: 6,
86 hitInterval: 0.5,
87 strikeVelocityMin: 0.5,
89 contactStiffness: 5e8,
90 contactExponent: 1.5,
91 contactRadius: 5e-3,
92 listenerDist: 1.0,
93 listenerAzimuth: 0,
94 listenerElevation: 0,
95 spacing: 0.008,
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 97 rigidProjection: false,
99 fps: 120,
100 duration: 4,
101};
103export const SPEED_OF_SOUND = 343; // [m/s]
104export const RHO_AIR = 1.2; // [kg/m^3]
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 105export const CFL = 0.9; // fraction of the stability limit used for the time step
106// The ring buffer of the retarded-time listener sum holds one batch of samples
107// plus this margin, which bounds the spread of the delays across the sheet.
108export const HIST_MARGIN = 128; // [samples]
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 109export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 110export const RB_INTERVAL_HELD = 16; // while the finger constraint is active
111// The held sheet counts as settled when its non-rigid kinetic energy drops
112// below settleKE, or when it has stopped falling: it has not halved over the
113// last SETTLE_WINDOW seconds (or 3 / settleSigma0, whichever is longer).
114export const SETTLE_WINDOW = 0.05; // [s]
115export const SETTLE_RATIO = 0.5;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 116
117// ---------- Grid ----------
119export interface Grid {
120 nx: number;
121 ny: number;
122 hx: number;
123 hy: number;
124 n: number;
125 /** 1 for nodes that belong to the sheet, 0 otherwise; absent for a full rectangle. */
126 mask?: Uint8Array;
127}
129export function makeGrid(p: PhysParams): Grid {
130 const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
131 const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
132 const g: Grid = { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
133 return p.shape === "ellipse" ? withEllipseMask(g) : g;
134}
136/**
137 * Restricts a grid to the ellipse inscribed in its bounding box. The edge is a
138 * staircase of whole cells: a cell belongs to the sheet when its center lies
139 * inside the ellipse, and a node when it is a corner of such a cell. Deciding
140 * by cell centers rather than by nodes keeps the covered area unbiased; a
141 * node-based mask loses about half a cell all around the edge, which raises
142 * every frequency by a few percent. A node with no cell of its own would carry
143 * no mass, and the trailing loop removes any such node.
144 */
145export function withEllipseMask(g: Grid): Grid {
146 const mask = new Uint8Array(g.n);
147 const a = ((g.nx - 1) * g.hx) / 2, b = ((g.ny - 1) * g.hy) / 2;
148 for (let cj = 0; cj < g.ny - 1; cj++) {
149 for (let ci = 0; ci < g.nx - 1; ci++) {
150 const X = (nodeX(g, ci) + g.hx / 2) / a, Y = (nodeY(g, cj) + g.hy / 2) / b;
151 if (X * X + Y * Y <= 1) {
152 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
153 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
154 }
155 }
156 }
157 const gm: Grid = { ...g, mask };
158 for (let changed = true; changed; ) {
159 changed = false;
160 for (let j = 0; j < g.ny; j++) {
161 for (let i = 0; i < g.nx; i++) {
162 if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
163 mask[j * g.nx + i] = 0;
164 changed = true;
165 }
166 }
167 }
168 }
169 return gm;
170}
172/** Node coordinates relative to the center of the bounding box. */
173export function nodeX(g: Grid, i: number): number {
174 return (i - (g.nx - 1) / 2) * g.hx;
175}
176export function nodeY(g: Grid, j: number): number {
177 return (j - (g.ny - 1) / 2) * g.hy;
178}
180/** True when node (i, j) exists and belongs to the sheet. */
181export function insideNode(g: Grid, i: number, j: number): boolean {
182 if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
183 return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
184}
186/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the sheet. */
187export function insideCell(g: Grid, ci: number, cj: number): boolean {
188 return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
189}
191/**
192 * Quadrature and lumped-mass weight of a node: the fraction of its four
193 * surrounding cells that belong to the sheet. On a full rectangle this is the
194 * trapezoidal rule (1 inside, 1/2 on an edge, 1/4 at a corner). Zero outside.
195 */
196export function nodeWeight(g: Grid, i: number, j: number): number {
197 if (!insideNode(g, i, j)) return 0;
198 let c = 0;
199 if (insideCell(g, i, j)) c++;
200 if (insideCell(g, i - 1, j)) c++;
201 if (insideCell(g, i, j - 1)) c++;
202 if (insideCell(g, i - 1, j - 1)) c++;
203 return c / 4;
204}
206/** Per-node weights as a flat array. */
207export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
208 const w = new Float32Array(g.n);
209 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);
210 return w;
211}
213/** Number of nodes in 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}
221// ---------- Material ----------
223export interface Material {
224 D: number; // bending stiffness E h^3 / (12 (1 - nu^2)) [N m]
225 C: number; // membrane stiffness E h / (1 - nu^2) [N/m]
226 G: number; // membrane shear stiffness E h / (2 (1 + nu)) [N/m]
227 nu: number;
228 rhoh: number; // mass per unit area [kg/m^2]
229}
231export function makeMaterial(p: PhysParams): Material {
232 return {
233 D: (p.E * p.h ** 3) / (12 * (1 - p.nu * p.nu)),
234 C: (p.E * p.h) / (1 - p.nu * p.nu),
235 G: (p.E * p.h) / (2 * (1 + p.nu)),
236 nu: p.nu,
237 rhoh: p.rho * p.h,
238 };
239}
241/**
242 * 60 dB decay time [s] of a bending mode of angular frequency omega under the
243 * loss model sigma0 + sigma1 k^2, using the flat-plate dispersion k^2 = omega / kappa.
244 */
245export function t60(m: Material, p: PhysParams, omega: number): number {
246 const kappa = Math.sqrt(m.D / m.rhoh);
247 const sigma = p.sigma0 + p.sigma1 * (omega / kappa);
248 return (3 * Math.LN10) / sigma;
249}
251/**
252 * Frequency [Hz] at which the three-point second difference lowers a bending
253 * frequency by the relative amount `err` on a grid of spacing hmax. A plane
254 * bending wave along a grid axis with wavenumber k has the discrete frequency
255 * ratio (2 - 2 cos kh) / (kh)^2.
256 */
257export function dispersionFrequency(m: Material, g: Grid, err: number): number {
258 const h = Math.max(g.hx, g.hy);
259 let lo = 0, hi = Math.PI;
260 for (let it = 0; it < 60; it++) {
261 const x = (lo + hi) / 2;
262 if ((2 - 2 * Math.cos(x)) / (x * x) > 1 - err) lo = x; else hi = x;
263 }
264 const k = lo / h;
265 const kappa = Math.sqrt(m.D / m.rhoh);
266 return (kappa * k * k) / (2 * Math.PI);
267}
269// ---------- Geometry of the excitation ----------
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 271/** One drumstick hit: time after the first impact [s] and stick speed [m/s]. */
272export interface Hit {
273 t: number;
274 v: number;
275}
277/**
278 * Schedule of the drumstick hits. The first hit is always at t = 0. A series
279 * has hitCount hits every hitInterval seconds, all at strikeVelocity for
280 * "steady", rising from strikeVelocityMin to strikeVelocity for "crescendo",
281 * and falling for "decrescendo".
282 */
283export function hitSchedule(p: PhysParams): Hit[] {
284 if (p.hitSequence === "single") return [{ t: 0, v: p.strikeVelocity }];
285 const count = Math.max(1, Math.round(p.hitCount));
286 const dt = Math.max(0, p.hitInterval);
287 const lo = Math.min(p.strikeVelocityMin, p.strikeVelocity), hi = p.strikeVelocity;
288 const hits: Hit[] = [];
289 for (let k = 0; k < count; k++) {
290 const f = count > 1 ? k / (count - 1) : 1;
291 let v = hi;
292 if (p.hitSequence === "crescendo") v = lo + (hi - lo) * f;
293 else if (p.hitSequence === "decrescendo") v = hi - (hi - lo) * f;
294 hits.push({ t: k * dt, v });
295 }
296 return hits;
297}
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 299/** Strike point in sheet coordinates. */
300export function strikePoint(p: PhysParams): [number, number] {
301 return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
302}
304/**
305 * Center of the thread hole. Rectangle: inset from the two edges of the
306 * top-left corner. Ellipse: inset along the radius from the rim point in the
307 * top-left direction (parametric angle 135 degrees).
308 */
309export function holePoint(p: PhysParams): [number, number] {
310 if (p.shape === "ellipse") {
311 const c = Math.SQRT1_2;
312 return [-(p.Lx / 2 - p.holeInset) * c, (p.Ly / 2 - p.holeInset) * c];
313 }
314 return [-p.Lx / 2 + p.holeInset, p.Ly / 2 - p.holeInset];
315}
317/** Center of the pushed patch: the bottom-right corner, or the rim point opposite the hole. */
318export function pushPoint(p: PhysParams): [number, number] {
319 if (p.shape === "ellipse") {
320 const c = Math.SQRT1_2;
321 return [(p.Lx / 2) * c, -(p.Ly / 2) * c];
322 }
323 return [p.Lx / 2, -p.Ly / 2];
324}
326/**
327 * Nodes of the drumstick contact patch: the sheet node nearest the strike
328 * point, plus every sheet node within the contact radius. The contact force
329 * is spread over these nodes in proportion to their mass, and the contact
330 * coordinate is their mass-weighted mean displacement, so the two are
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 331 * conjugate and the contact has a potential energy. The mass-weighted center
332 * of the patch is where the inertia-relief reaction to the force is taken.
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 334export function contactPatch(p: PhysParams, g: Grid): { nodes: number[]; nearest: number; mass: number; center: [number, number]; box: [number, number, number, number] } {
336 const [sx, sy] = strikePoint(p);
337 let nearest = -1, best = Infinity;
338 for (let j = 0; j < g.ny; j++) {
339 for (let i = 0; i < g.nx; i++) {
340 if (nodeWeight(g, i, j) === 0) continue;
341 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
342 if (d2 < best) { best = d2; nearest = j * g.nx + i; }
343 }
344 }
345 const nodes: number[] = [];
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 346 let mass = 0, cx = 0, cy = 0;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 347 const box: [number, number, number, number] = [g.nx, -1, g.ny, -1]; // i0, i1, j0, j1 (inclusive)
348 for (let j = 0; j < g.ny; j++) {
349 for (let i = 0; i < g.nx; i++) {
350 const k = j * g.nx + i;
351 const wgt = nodeWeight(g, i, j);
352 if (wgt === 0) continue;
353 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
354 if (d2 <= p.contactRadius ** 2 || k === nearest) {
355 nodes.push(k);
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 356 const mk = m.rhoh * g.hx * g.hy * wgt;
357 mass += mk; cx += mk * nodeX(g, i); cy += mk * nodeY(g, j);
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 358 box[0] = Math.min(box[0], i); box[1] = Math.max(box[1], i);
359 box[2] = Math.min(box[2], j); box[3] = Math.max(box[3], j);
360 }
361 }
362 }
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 363 return { nodes, nearest, mass, center: [cx / mass, cy / mass], box };