Cymbal: free plate with rigid-body motion removed (no support); lighter default damping - proof-of-concept/toy-instruments

proof-of-concept / toy-instruments

Cymbal: free plate with rigid-body motion removed (no support); lighter default damping

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012zsyegMfQnzJEzgMiFLzew

Jeremy Magland <jeremy.magland@gmail.com> committed 12 days ago commit 591834d83045 parent 15fbb00 Browse files

8 changed files+100−42

README.mdmodified+1−1View file

@@ -39,7 +39,7 @@ String tests (`test/string.test.ts`): modes vs the discrete-string formula; smal

3939 - Thin bronze plate (E = 100 GPa, 8600 kg/m^3, 0.8 mm by default), round, free rim, centre hole; flat by default, optional bow (dome, y = h (1 - (r/R)^2)) and bell (raised cosine). Beads on a triangular lattice lifted to the profile; the mesh is scaled so its area equals the disk's.

4040 - In-plane: springs with k = sqrt(3)/2 E t / softening (an edge with one triangle gets half). Softening (default 30x in the fast preset) speeds the run; bending is unaffected.

4141 - Bending: dihedral angles across interior edges (`src/sim/hinge.ts`, Bridson et al. 2003), sign-normalized so a positive change folds toward the upward normal. Each triangle builds a shape operator S = sum_e phi_e l_e / (2A) m_e m_e^T and stores the plate energy A D/2 [(tr S)^2 - 2 (1 - nu) det S] (nu = 1/3) as a 3x3 matrix in its hinge angles, condensed over free edges so they carry no moment. Plain hinge springs get the twisting term wrong on a regular mesh (effective nu = -1/3), which matters for a free rim. Damping uses the same matrices (Rayleigh, sigma1 omega^2).

42-- Centre support: a small rigid hub. Beads within the hub radius (default: the first ring around the hole) are held fixed; no springs, no added damping, so the rest rings like a free plate. At least three beads are needed, or the plate can tilt freely about the support.

42+- Support: none. The plate is free; its rigid-body motion (translation and rotation, mass-weighted about the current centre of mass) is removed from the accelerations at every force evaluation and from the velocities after every step (`BeadSystem.freeRigid`). The elastic modes are mass-orthogonal to the rigid ones, so they are unchanged; nothing stiffens or damps them. Pinning only the centre bead was tried: an off-centre hit sets the plate tumbling (tilts past 89 degrees within a second).

4343 - Excitation: drumstick (Hertz, K delta^1.5) or felt mallet as a 3D point mass with contact along the surface normal, regularized Coulomb friction and a contact reach; hit angle tilts the stroke toward the rim. Also strike and force pulse.

4444 - Sound: far-field microphone via the Rayleigh integral, p = rho / (2 pi) sum A_i a_n,i(t - R_i / c) / R_i (baffled-plate assumption: no diffraction around the rim).

4545 - Tests (`test/cymbal.test.ts`, `test/hinge.test.ts`): hinge gradient vs finite differences; bending energy matches the plate for cylinder, sphere, saddle and twist to 0.1%; a free flat plate's (2,0) mode within 3.4% of Leissa's value on a 9-ring lattice; energy balance with a glancing stick, friction, felts and damping.

src/instruments/bead-cymbal/model.tsmodified+5−19View file

@@ -3,8 +3,9 @@ import { hingeAngleGrad } from '../../sim/hinge';

33 import { dampingFromT60 } from '../bead-string/model';

44

55 // A cymbal: a thin bronze plate, round, free at the rim, with a hole in the

6-// centre; the beads around the hole are held fixed (a small rigid hub), so the

7-// rest of the plate rings like a free plate, with no added damping. Beads on a triangular lattice in

6+// centre. Nothing holds it: it is simulated as a free plate whose rigid-body

7+// motion is removed (see BeadSystem.freeRigid), so it rings with no support

8+// stiffness or damping. Beads on a triangular lattice in

89 // the x-z plane, lifted to the profile height y(r) (flat by default; optional

910 // bow and bell). In-plane springs on every edge; plate bending from the

1011 // dihedral angles across interior edges, per triangle (see Bending in sim/types). Vibration is mostly along y.

@@ -24,8 +25,6 @@ export interface CymbalParams {

2425 youngs: number; // Pa

2526 density: number; // kg/m^3

2627 soften: number; // in-plane stiffness divided by this (1 = real bronze)

27- // centre support: beads within this radius are held fixed (a rigid hub)

28- hubRadius: number; // m

2928 // damping

3029 t60Low: number;

3130 fLow: number;

@@ -128,7 +127,6 @@ export interface CymbalDerived {

128127 mass: Float64Array;

129128 k: number; // interior spring stiffness

130129 D: number; // bending stiffness

131- support: number[]; // beads held fixed at the hub

132130 sigma0: number;

133131 sigma1: number;

134132 waveSpeedP: number; // in-plane P waves (after softening)

@@ -176,17 +174,7 @@ export function deriveCymbal(p: CymbalParams): CymbalDerived {

176174 const k = (Math.sqrt(3) / 2) * Et;

177175 const nu = 1 / 3;

178176 const D = (p.youngs * p.thickness ** 3) / (12 * (1 - nu * nu));

179- const support: number[] = [];

180- mesh.pos.forEach(([x, z], i) => {

181- if (Math.hypot(x, z) <= p.hubRadius + 1e-9) support.push(i);

182- });

183177 const warnings: string[] = [];

184- if (support.length < 3)

185- warnings.push(

186- support.length

187- ? 'Fewer than three beads are held at the hub: the plate can tilt freely about them and will swing or spin after an off-centre hit. Increase the hub radius.'

188- : 'No beads are held at the hub: the plate is free and will drift away after a hit. Increase the hub radius.',

189- );

190178 if (p.soften > 1)

191179 warnings.push(

192180 `In-plane stiffness is softened ${p.soften}×. A flat plate bends the same at small amplitude; curved shapes lose curvature stiffening and big hits exaggerate the nonlinear coupling.`,

@@ -203,7 +191,6 @@ export function deriveCymbal(p: CymbalParams): CymbalDerived {

203191 mass,

204192 k,

205193 D,

206- support,

207194 sigma0,

208195 sigma1,

209196 waveSpeedP: Math.sqrt(lambdaPlus2Mu / rhoT),

@@ -235,15 +222,14 @@ export function buildCymbal(d: CymbalDerived): BeadSystem {

235222 sL0.push(len(e.i, e.j));

236223 sc.push(0);

237224 }

238- const pinned = new Uint8Array(n);

239- for (const i of d.support) pinned[i] = 1;

240225 const bending = buildBending(d, 1 / 3);

241226 return {

242227 n,

243228 rest,

244229 mass: d.mass,

245- pinned,

230+ pinned: new Uint8Array(n),

246231 anchors: [],

232+ freeRigid: true,

247233 si: Int32Array.from(si),

248234 sj: Int32Array.from(sj),

249235 sk: Float64Array.from(sk),

src/instruments/bead-cymbal/page.tsmodified+4−16View file

@@ -38,19 +38,6 @@ const SECTIONS: Section<CymbalParams>[] = [

3838 { key: 'soften', label: 'In-plane softening', kind: 'number', unit: '×', help: 'Divide the in-plane stiffness by this for speed (bending stays real). 1 = real bronze. Flat plates are unaffected at small amplitude; curved shapes and big hits are.' },

3939 ],

4040 },

41- {

42- title: 'Centre support',

43- params: [

44- {

45- key: 'hubRadius',

46- label: 'Hub radius',

47- kind: 'number',

48- unit: 'mm',

49- scale: 1000,

50- help: 'Beads within this radius are held fixed: a small rigid hub, no springs or damping. Keep it small (just past the first ring of beads) for a nearly free plate.',

51- },

52- ],

53- },

5441 {

5542 title: 'Damping (0 = off)',

5643 params: [

@@ -128,7 +115,7 @@ export function mountBeadCymbal(root: HTMLElement) {

128115 excKind: p.excitation,

129116 pickup: mic.weights,

130117 pickupLabel: 'microphone',

131- post: { y: profileHeight(p, 0) },

118+ note: 'free plate: nothing holds it, rigid motion removed',

132119 mic: {

133120 dir: micDirection(p),

134121 from: [0, profileHeight(p, 0), 0],

@@ -148,12 +135,13 @@ export function mountBeadCymbal(root: HTMLElement) {

148135 const d = deriveCymbal(p);

149136 const R = p.diameter / 2;

150137 const rows: [string, string][] = [

151- ['Beads', `${d.n} (${d.support.length} held at the hub)`],

138+ ['Beads', `${d.n}`],

139+ ['Support', 'none: a free plate with its rigid-body motion removed'],

152140 ['Lattice spacing', `${fmtNum(d.mesh.a * 1000)} mm`],

153141 ['Mass', `${fmtNum(d.mass.reduce((s, m) => s + m, 0) * 1000)} g`],

154142 ['Bending stiffness D', `${fmtNum(d.D)} N·m`],

155143 ['In-plane wave speed', `${fmtNum(d.waveSpeedP)} m/s${p.soften > 1 ? ` (softened ${fmtNum(p.soften)}×)` : ''}`],

156- ['Free flat plate (2,0) mode', `${fmtNum((5.253 / (2 * Math.PI * R * R)) * d.bendingCoef)} Hz (ideal, no support)`],

144+ ['Free flat plate (2,0) mode', `${fmtNum((5.253 / (2 * Math.PI * R * R)) * d.bendingCoef)} Hz (ideal)`],

157145 ['Bending resolved to', `~${fmtNum(bendingLimitHz(d))} Hz`],

158146 ];

159147 const warnings = [...d.warnings];

src/instruments/bead-cymbal/presets.tsmodified+3−4View file

@@ -12,10 +12,9 @@ const toy: CymbalParams = {

1212 youngs: 1.0e11,

1313 density: 8600,

1414 soften: 30,

15- hubRadius: 0.011, // the first ring of beads around the hole

16- t60Low: 4,

15+ t60Low: 8,

1716 fLow: 200,

18- t60High: 1,

17+ t60High: 4,

1918 fHigh: 2000,

2019 excitation: 'stick',

2120 excR: 0.85,

@@ -52,7 +51,7 @@ export const CYMBAL_PRESETS: CymbalPreset[] = [

5251 {

5352 name: 'Flat 8" toy, stick (fast)',

5453 description:

55- 'Flat 0.8 mm bronze plate, 20 cm, ~370 beads, held only at a small hub (the 6 beads around the hole). In-plane stiffness softened 30× for speed (bending is real). Stick hit near the edge, spread over a 15 mm patch (a single bead would make the stick chatter). Bending is resolved to about 1.3 kHz.',

54+ 'Flat 0.8 mm bronze plate, 20 cm, ~370 beads, free (nothing holds it; its rigid motion is removed). In-plane stiffness softened 30× for speed (bending is real). Stick hit near the edge, spread over a 15 mm patch (a single bead would make the stick chatter). Bending is resolved to about 1.3 kHz.',

5655 params: toy,

5756 },

5857 {

src/sim/solver.tsmodified+56−0View file

@@ -357,8 +357,53 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog

357357 A[3 * i + 1] = F[3 * i + 1] * im;

358358 A[3 * i + 2] = F[3 * i + 2] * im;

359359 }

360+ if (sys.freeRigid) removeRigid(A);

360361 };

361362

363+ // subtract the rigid-body part (mass-weighted) of a per-bead field u;

364+ // returns the kinetic energy of that part when u is a velocity

365+ const removeRigid = (u: Float64Array): number => {

366+ let M = 0, cx = 0, cy = 0, cz = 0;

367+ for (let i = 0; i < n; i++) {

368+ M += mass[i];

369+ cx += mass[i] * X[3 * i];

370+ cy += mass[i] * X[3 * i + 1];

371+ cz += mass[i] * X[3 * i + 2];

372+ }

373+ cx /= M;

374+ cy /= M;

375+ cz /= M;

376+ let px = 0, py = 0, pz = 0, lx = 0, ly = 0, lz = 0;

377+ let ixx = 0, iyy = 0, izz = 0, ixy = 0, ixz = 0, iyz = 0;

378+ for (let i = 0; i < n; i++) {

379+ const m = mass[i];

380+ const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;

381+ const ux = u[3 * i], uy = u[3 * i + 1], uz = u[3 * i + 2];

382+ px += m * ux;

383+ py += m * uy;

384+ pz += m * uz;

385+ lx += m * (ry * uz - rz * uy);

386+ ly += m * (rz * ux - rx * uz);

387+ lz += m * (rx * uy - ry * ux);

388+ ixx += m * (ry * ry + rz * rz);

389+ iyy += m * (rx * rx + rz * rz);

390+ izz += m * (rx * rx + ry * ry);

391+ ixy -= m * rx * ry;

392+ ixz -= m * rx * rz;

393+ iyz -= m * ry * rz;

394+ }

395+ const [wx, wy, wz] = solve3(ixx, ixy, ixz, iyy, iyz, izz, lx, ly, lz);

396+ const vx = px / M, vy = py / M, vz = pz / M;

397+ for (let i = 0; i < n; i++) {

398+ const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;

399+ u[3 * i] -= vx + (wy * rz - wz * ry);

400+ u[3 * i + 1] -= vy + (wz * rx - wx * rz);

401+ u[3 * i + 2] -= vz + (wx * ry - wy * rx);

402+ }

403+ return 0.5 * (px * vx + py * vy + pz * vz) + 0.5 * (wx * lx + wy * ly + wz * lz);

404+ };

405+ if (sys.freeRigid) removeRigid(V); // an initial strike carries no rigid motion either

406+

362407 const velocityForces =

363408 sys.sigma0 !== 0 || sc.some((c) => c !== 0) || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;

364409 let maxStretch = 0;

@@ -458,6 +503,8 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog

458503 workIn += dt * pExt;

459504 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];

460505 if (mallet) for (let c = 0; c < 3; c++) qv[c] += 0.5 * dt * qa[c];

506+ // the rigid-free subspace turns with the body; drop what drifts out of it

507+ if (sys.freeRigid) dissipated += removeRigid(V);

461508 syncAnchors();

462509 }

463510 // at a sample time: recompute with the full-step velocity

@@ -481,3 +528,12 @@ export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, prog

481528 malletPos,

482529 };

483530 }

531+

532+// solve the symmetric 3x3 system [a b c; b d e; c e f] x = (r0, r1, r2)

533+function solve3(a: number, b: number, c: number, d: number, e: number, f: number, r0: number, r1: number, r2: number): [number, number, number] {

534+ const A00 = d * f - e * e, A01 = c * e - b * f, A02 = b * e - c * d;

535+ const A11 = a * f - c * c, A12 = b * c - a * e, A22 = a * d - b * b;

536+ const det = a * A00 + b * A01 + c * A02;

537+ if (!(Math.abs(det) > 0)) return [0, 0, 0];

538+ return [(A00 * r0 + A01 * r1 + A02 * r2) / det, (A01 * r0 + A11 * r1 + A12 * r2) / det, (A02 * r0 + A12 * r1 + A22 * r2) / det];

539+}

src/sim/types.tsmodified+4−0View file

@@ -21,6 +21,10 @@ export interface BeadSystem {

2121 sL0: Float64Array; // rest length (m)

2222 sc: Float64Array; // dashpot on relative velocity, all 3 components (N s/m)

2323 sigma0: number; // mass-proportional damping: force -2 sigma0 m v

24+ // remove rigid-body motion (translation and rotation, mass-weighted) from the

25+ // beads' accelerations and velocities: a free body whose rigid modes are not

26+ // simulated. The elastic modes are unchanged (they are mass-orthogonal to them).

27+ freeRigid?: boolean;

2428 bending?: Bending;

2529 }

2630

src/ui/scene3d.tsmodified+2−0View file

@@ -16,6 +16,7 @@ export interface SceneSpec {

1616 pickupLabel: string;

1717 hoop?: { radius: number }; // round rim in the x-z plane, centred on the origin (instead of walls)

1818 post?: { y: number }; // centre support: a stand rod from (0, y, 0) down to the floor

19+ note?: string; // extra legend text

1920 mic?: { dir: [number, number, number]; label: string; from: [number, number, number] }; // direction to a far-field microphone

2021 }

2122

@@ -439,6 +440,7 @@ export class BeadScene {

439440 this.pickupArrow = { arrow, bead: i };

440441 }

441442 legend.push(spec.mic ? `<span>${sw(cPick, 'arrow')}${spec.mic.label}</span>` : `<span>${sw(cPick)}pickup: ${spec.pickupLabel}</span>`);

443+ if (spec.note) legend.push(`<span>${spec.note}</span>`);

442444 this.legend.innerHTML = legend.join('');

443445 this.update();

444446 }

test/cymbal.test.tsmodified+25−2View file

@@ -39,7 +39,7 @@ describe('cymbal lattice', () => {

3939

4040 it('a free flat plate rings near the (2,0) mode of a free circular plate', () => {

4141 // no felts, no hole, no damping; initial shape r^2 cos(2 phi), small

42- const p = base({ holeRadius: 0, rings: 9, hubRadius: -1, t60Low: 0, t60High: 0, duration: 0.3, frameRate: 50 });

42+ const p = base({ holeRadius: 0, rings: 9, t60Low: 0, t60High: 0, duration: 0.3, frameRate: 50 });

4343 const d = deriveCymbal(p);

4444 const sys = buildCymbal(d);

4545 const exc = buildCymbalExcitation({ ...p, excitation: 'strike', strikeVel: 0 }, d);

@@ -61,7 +61,7 @@ describe('cymbal lattice', () => {

6161 expect(Math.abs(s.freq[best] / ideal - 1)).toBeLessThan(0.05);

6262 }, 60000);

6363

64- it('energy balances with a glancing stick hit, friction, hub and damping', () => {

64+ it('energy balances with a glancing stick hit, friction, damping and rigid motion removed', () => {

6565 const p = base({ rings: 7, hitAngle: 60, stickSpeed: 6, bowHeight: 0.01, bellHeight: 0.01, duration: 0.03, frameRate: 2000 });

6666 const d = deriveCymbal(p);

6767 const mic = buildMic(p, d);

@@ -74,4 +74,27 @@ describe('cymbal lattice', () => {

7474 expect(cf[cf.length - 1]).toBe(0);

7575 expect(r.audio.some((v) => v !== 0)).toBe(true);

7676 });

77+

78+ it('removing rigid motion leaves the plate with no net momentum and keeps the hit point vibrating', () => {

79+ const p = base({ rings: 7, duration: 0.05, frameRate: 2000 });

80+ const d = deriveCymbal(p);

81+ const sys = buildCymbal(d);

82+ const mic = buildMic(p, d);

83+ const r = simulate(sys, buildCymbalExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay });

84+ // mean vertical displacement (rigid translation) stays near zero compared with the vibration

85+ const n = sys.n;

86+ let meanMax = 0, ampMax = 0;

87+ for (let f = 0; f < r.frameTimes.length; f++) {

88+ let m = 0, mt = 0;

89+ for (let i = 0; i < n; i++) {

90+ m += sys.mass[i] * r.frames[f * 3 * n + 3 * i + 1];

91+ mt += sys.mass[i];

92+ ampMax = Math.max(ampMax, Math.abs(r.frames[f * 3 * n + 3 * i + 1]));

93+ }

94+ meanMax = Math.max(meanMax, Math.abs(m / mt));

95+ }

96+ console.log('mean y', (meanMax * 1000).toFixed(4), 'mm, max |y|', (ampMax * 1000).toFixed(3), 'mm');

97+ expect(meanMax / ampMax).toBeLessThan(1e-2);

98+ expect(Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy)).toBeLessThan(1e-2);

99+ }, 60000);

77100 });