5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// The float64 reference discretization: symmetry and rigid modes of the
2// linear stiffness, forces as the gradient of the energy, dissipativity,
3// thickness scaling, the stand, and the free-plate frequencies.
4import { test } from "node:test";
5import assert from "node:assert/strict";
6import { makeGrid, withMask, type Grid, type Plate, nodeMasses, standRates, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 7import { uniformPlate, accelerations, addDamping, addPadMemory, padState, stepPadMemory, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 8import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
9import { fullLoss } from "../src/model/reference.ts";
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 10import { feltModulus } from "../src/felt.ts";
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 11
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 12const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 13
14function rnd(seed: number): () => number {
15 let s = seed >>> 0;
16 return () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296 - 0.5; };
17}
18function randomFields(pl: Plate, amp: number, seed: number): Fields {
19 const g = pl.grid, r = rnd(seed), f = fields(g.n);
20 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
21 if (!insideNode(g, i, j)) continue;
22 const k = j * g.nx + i;
23 f.u[k] = 0.02 * amp * r(); f.v[k] = 0.02 * amp * r(); f.w[k] = amp * r();
24 }
25 return f;
26}
27function bowed(g: Grid, sag: number): Float64Array {
28 const w0 = new Float64Array(g.n);
29 const Lx = (g.nx - 1) * g.hx;
30 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = sag * (1 - (2 * nodeX(g, i) / Lx) ** 2) + 0.3 * sag * Math.sin(3 * nodeY(g, j) / Lx);
31 return w0;
32}
33function withStand(pl: Plate): Plate {
34 return { ...pl, stand: { x: 0, y: 0, pad: { innerRadius: 0.005, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
35}
37test("the linear stiffness is symmetric with exactly three rigid modes per block (uniform, tapered, disc with hole)", () => {
38 const cases: { name: string; pl: Plate; fields: 1 | 3 }[] = [];
39 const g1 = makeGrid(0.12, 0.09, 0.015);
40 cases.push({ name: "uniform rectangle", pl: uniformPlate(g1, steel, 1e-3), fields: 3 });
41 const tapered = uniformPlate(g1, steel, 1e-3);
42 for (let k = 0; k < g1.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.5 * (k % g1.nx) / g1.nx);
43 cases.push({ name: "tapered rectangle", pl: tapered, fields: 3 });
44 const gd = withMask(makeGrid(0.16, 0.16, 0.01), (x, y) => x * x + y * y <= 0.08 ** 2 && x * x + y * y >= 0.015 ** 2);
45 cases.push({ name: "disc with hole", pl: uniformPlate(gd, steel, 1e-3), fields: 1 });
46 for (const c of cases) {
47 const sys = assemble(c.pl, c.fields);
48 // Symmetry is built into the banded storage; check the rigid modes: K x = 0 for translations and rotations.
49 const g = c.pl.grid;
50 const nullVectors: Float64Array[] = [];
51 const mk = (fn: (i: number, j: number, comp: number) => number) => {
52 const x = new Float64Array(sys.n);
53 for (let d = 0; d < sys.n; d++) {
54 const k = sys.dofNode[d];
55 const comp = sys.fields === 1 ? 2 : d % 3;
56 x[d] = fn(k % g.nx, Math.floor(k / g.nx), comp);
57 }
58 return x;
59 };
60 nullVectors.push(mk((_i, _j, c) => (c === 2 ? 1 : 0)));
61 nullVectors.push(mk((i, _j, c) => (c === 2 ? nodeX(g, i) : 0)));
62 nullVectors.push(mk((_i, j, c) => (c === 2 ? nodeY(g, j) : 0)));
63 if (c.fields === 3) {
64 nullVectors.push(mk((_i, _j, c) => (c === 0 ? 1 : 0)));
65 nullVectors.push(mk((_i, _j, c) => (c === 1 ? 1 : 0)));
66 nullVectors.push(mk((i, j, c) => (c === 0 ? -nodeY(g, j) : c === 1 ? nodeX(g, i) : 0)));
67 }
68 const y = new Float64Array(sys.n);
69 for (const x of nullVectors) {
70 const { stiffnessMatvec } = require_assemble();
71 stiffnessMatvec(sys, x, y);
72 let ny = 0, nx = 0;
73 for (let d = 0; d < sys.n; d++) { ny += y[d] * y[d]; nx += x[d] * x[d]; }
74 const scale = Math.abs(stiffnessEntry(sys, 0, 0)) * Math.sqrt(nx);
75 assert.ok(Math.sqrt(ny) < 1e-9 * scale, `${c.name}: rigid mode residual ${Math.sqrt(ny) / scale}`);
76 }
77 // A non-rigid vector must not be in the null space.
78 const z = mk((i, j, c) => (c === 2 ? Math.sin(i * 0.9) * Math.cos(j * 1.1) : 0));
79 const { stiffnessMatvec } = require_assemble();
80 stiffnessMatvec(sys, z, y);
81 assert.ok(y.some((v) => Math.abs(v) > 1e-3), `${c.name}: bending is not trivial`);
82 }
83});
85import * as assembleModule from "../src/model/assemble.ts";
86function require_assemble(): typeof assembleModule { return assembleModule; }
88test("forces are the exact gradient of the discrete energy (flat and bowed, tapered, free and on a stand)", () => {
89 const g = makeGrid(0.1, 0.08, 0.01);
90 const variants: { name: string; pl: Plate }[] = [];
91 variants.push({ name: "flat free", pl: uniformPlate(g, steel, 1e-3) });
92 variants.push({ name: "bowed free", pl: uniformPlate(g, steel, 1e-3, bowed(g, 3e-3)) });
93 const tapered = uniformPlate(g, steel, 1e-3, bowed(g, 2e-3));
94 for (let k = 0; k < g.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.4 * (k % g.nx) / g.nx);
95 variants.push({ name: "bowed tapered", pl: tapered });
96 variants.push({ name: "bowed on a stand", pl: withStand(uniformPlate(g, steel, 1e-3, bowed(g, 3e-3))) });
97 for (const v of variants) {
98 const pl = v.pl;
99 const sr = standRates(pl);
100 const f = randomFields(pl, 2.5e-3, 7);
101 const masses = nodeMasses(pl);
102 const acc = fields(g.n);
103 accelerations(pl, f, acc, false, sr);
104 const dir = randomFields(pl, 1, 11);
105 // Directional derivative of U along dir versus -F . dir.
106 const eps = 1e-7;
107 const shift = (s: number): Fields => ({ u: f.u.map((x, k) => x + s * dir.u[k]), v: f.v.map((x, k) => x + s * dir.v[k]), w: f.w.map((x, k) => x + s * dir.w[k]) });
108 const dU = (potentialEnergy(pl, shift(eps), sr) - potentialEnergy(pl, shift(-eps), sr)) / (2 * eps);
109 let fd = 0;
110 for (let k = 0; k < g.n; k++) fd -= masses[k] * (acc.u[k] * dir.u[k] + acc.v[k] * dir.v[k] + acc.w[k] * dir.w[k]);
111 assert.ok(Math.abs(dU - fd) < 1e-6 * Math.abs(dU), `${v.name}: dU ${dU} vs -F.dir ${fd}`);
112 }
113});
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 115test("the loss terms only remove energy", () => {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 116 const g = makeGrid(0.1, 0.08, 0.01);
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 117 const pl = uniformPlate(g, steel, 1e-3);
118 const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, 20);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 119 const masses = nodeMasses(pl);
120 for (const seed of [1, 2, 3]) {
121 const vel = randomFields(pl, 1, seed);
122 const acc = fields(g.n);
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 123 addDamping(pl, loss, vel, acc);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 124 let power = 0;
125 for (let k = 0; k < g.n; k++) power += masses[k] * (acc.u[k] * vel.u[k] + acc.v[k] * vel.v[k] + acc.w[k] * vel.w[k]);
126 assert.ok(power < 0, `power ${power}`);
127 }
128});
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 130test("the felt's chain, stepped in time, dissipates the work of its loss stiffness per cycle", () => {
131 const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
132 const pl = withStand(uniformPlate(g, steel, 1e-3));
133 const sr = standRates(pl)!;
134 const masses = nodeMasses(pl);
135 // A rigid translation w = A sin(omega t) of the plate, stepped finely; the work done
136 // against the arms over a cycle, once the transient of the arms has passed, is
137 // pi L(omega) K A^2 with L the loss stiffness per unit of K (eta at the bounce).
138 for (const f of [sr.bounceHz, 2000]) {
139 const omega = 2 * Math.PI * f, A = 1e-4;
140 const steps = 400, dt = 1 / (f * steps);
141 const pad = padState(pl, sr);
142 const x = fields(g.n), dw = fields(g.n), acc = fields(g.n);
143 const setW = (t: number, out: Fields) => { for (let k = 0; k < g.n; k++) out.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? A * Math.sin(omega * t) : 0; };
144 let work = 0;
145 const cycles = 4;
146 for (let s = 0; s < cycles * steps; s++) {
147 const t0 = s * dt, t1 = (s + 1) * dt;
148 setW(t0, x); setW(t1, dw);
149 for (let k = 0; k < g.n; k++) dw.w[k] -= x.w[k];
150 acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
151 addPadMemory(pl, sr, x, pad, false, acc);
152 // Work of the plate against the arms over the step, with the force at mid-step (trapezoid).
153 let F0 = 0;
154 for (let k = 0; k < g.n; k++) F0 += masses[k] * acc.w[k];
155 stepPadMemory(pl, sr, dt, dw, x, pad, false);
156 setW(t1, x);
157 acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
158 addPadMemory(pl, sr, x, pad, false, acc);
159 let F1 = 0;
160 for (let k = 0; k < g.n; k++) F1 += masses[k] * acc.w[k];
161 if (s >= (cycles - 1) * steps) work -= 0.5 * (F0 + F1) * dw.w[sr.nodes[0]];
162 }
163 const L = feltModulus(sr.chain, omega).loss;
164 const expected = Math.PI * L * sr.totalKw * A * A;
165 assert.ok(work > 0);
166 assert.ok(Math.abs(work / expected - 1) < 0.02, `at ${f.toFixed(0)} Hz: work per cycle ${work} vs ${expected}`);
167 if (f === sr.bounceHz) assert.ok(Math.abs(L / sr.eta - 1) < 0.02, `loss stiffness at the bounce ${L} vs eta ${sr.eta}`);
168 }
169});
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 171test("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
172 const g = makeGrid(0.1, 0.1, 0.01);
173 const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
174 // The largest frequency is set by the in-plane waves, independent of h.
175 assert.ok(Math.abs(w1 / w2 - 1) < 1e-6, `in-plane omega_max ${w1} vs ${w2}`);
176 const { stiffnessEntry: ke } = assembleModule;
177 const s1 = assemble(uniformPlate(g, steel, 1e-3)), s2 = assemble(uniformPlate(g, steel, 2e-3));
178 const c = Math.floor(s1.n / 2);
179 // Bending stiffness scales as h^3 and the mass as h: omega^2 scales as h^2.
180 assert.ok(Math.abs(ke(s2, c, c) / ke(s1, c, c) - 8) < 1e-9);
181 assert.ok(Math.abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12);
182 assert.ok(Math.abs(bendingStiffness(steel, 2e-3) / bendingStiffness(steel, 1e-3) - 8) < 1e-12);
183});
185test("the stand resists a rigid translation with its total stiffness and a tilt with the pads' second moment", () => {
186 const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
187 const pl = withStand(uniformPlate(g, steel, 1e-3));
188 const sr = standRates(pl)!;
189 const masses = nodeMasses(pl);
190 const f = fields(g.n), acc = fields(g.n);
191 const dz = 1e-4;
192 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? dz : 0;
193 accelerations(pl, f, acc, false, sr);
194 let F = 0;
195 for (let k = 0; k < g.n; k++) F += masses[k] * acc.w[k];
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 196 // The relaxed felt holds kInf of the stiffness at the bounce; the arms of the chain hold the rest.
197 const kInf = sr.chain.kInf;
198 assert.ok(kInf > 0.2 && kInf < 0.8, `relaxed fraction ${kInf}`);
199 assert.ok(Math.abs(-F / dz / (kInf * sr.totalKw) - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${kInf * sr.totalKw}`);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 200 // Tilt about y: w = theta x; restoring moment = theta * sum_k kw_k x_k^2 (bending of a plane is
201 // zero; the von Karman membrane strain of a rotation is theta^2 / 2, so a tiny tilt keeps it negligible).
202 const th = 1e-7;
203 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? th * nodeX(g, k % g.nx) : 0;
204 accelerations(pl, f, acc, false, sr);
205 let M = 0, I2 = 0;
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 206 for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += kInf * sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 207 assert.ok(Math.abs(-M / th / I2 - 1) < 1e-6, `tilt stiffness ${-M / th} vs ${I2}`);
208 // With lift-off, a small motion within the preload is the same linear spring; a large upward motion unloads the bottom pad.
209 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 1e-5 : 0;
210 const lin = fields(g.n), lo = fields(g.n);
211 accelerations(pl, f, lin, false, sr);
212 accelerations(pl, f, lo, true, sr);
213 for (let k = 0; k < g.n; k++) assert.ok(Math.abs(lin.w[k] - lo.w[k]) < 1e-9 * Math.max(1, Math.abs(lin.w[k])));
214 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 0.01 : 0;
215 accelerations(pl, f, lo, true, sr);
216 let Flo = 0;
217 for (let k = 0; k < g.n; k++) Flo += masses[k] * lo.w[k];
218 // Bottom pad unloaded: only the top pad (compressed by preload + 10 mm) pushes down, minus the rest balance.
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 219 const kb = (kInf * sr.totalKw) / 2;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 220 const expected = -(kb * (sr.preloadTop + 0.01) - kb * sr.preloadTop) - kb * sr.preloadBottom;
221 assert.ok(Math.abs(Flo / expected - 1) < 1e-6, `lift-off force ${Flo} vs ${expected}`);
222});