982850fWoodblock strike simulator: Python reference model and interactive web appJeremy Magland 1// Time-domain RK4 integration of the coupled mallet / body / cavity system.
2// Port of simulate() in woodblock.py. State: [q(N), qdot(N), xi, xidot, y, ydot].
4import type { WoodblockParams } from '../params';
5import type { ModalBasisData } from './basis';
7export interface SolveResult {
8 Qdd: Float64Array[]; // N arrays of length T: modal accelerations
9 Xdd: Float64Array; // slit-plug acceleration
10 force: Float64Array; // contact force f_c(t)
11 fsHi: number;
12 helmholtzHz: number;
13 peakForceN: number;
14 contactMs: number;
15}
17export function simulate(
18 cfg: WoodblockParams,
19 basis: ModalBasisData,
20 oversample: number,
21 onProgress?: (frac: number) => void,
22): SolveResult {
23 const N = basis.n;
24 if (N === 0) throw new Error('No modes below max_mode_hz -- check geometry/material.');
26 const fsHi = cfg.sample_rate * oversample;
27 const dt = 1.0 / fsHi;
28 const nSteps = Math.floor(cfg.duration_s * fsHi);
30 const om2 = new Float64Array(N);
31 const twoZetaOm = new Float64Array(N);
32 for (let i = 0; i < N; i++) {
33 om2[i] = basis.omega[i] * basis.omega[i];
34 twoZetaOm[i] = 2.0 * basis.zeta[i] * basis.omega[i];
35 }
36 const phiC = basis.phiC;
37 const A = basis.A;
39 const rho0 = cfg.air_density;
40 const c = cfg.speed_of_sound;
41 const V = cfg.cavity_volume_m3;
42 const S = cfg.slit_area_m2;
43 // End correction on the slit: ~0.85 * equivalent radius per open end.
44 const Leff = cfg.slit_length_m + 2.0 * 0.85 * Math.sqrt(S / Math.PI);
45 const omegaH = c * Math.sqrt(S / (V * Leff));
46 const twoZetaHOm = 2.0 * cfg.cavity_damping * omegaH;
47 const omegaH2 = omegaH * omegaH;
48 const pScale = (rho0 * c * c) / V;
49 const cavDrive = (c * c) / (V * Leff);
51 const K = cfg.mallet_stiffness;
52 const alpha = cfg.contact_exponent;
53 const mu = cfg.contact_damping;
54 const invMh = 1.0 / cfg.mallet_mass;
56 const dim = 2 * N + 4;
57 const IX = 2 * N;
58 const IXD = 2 * N + 1;
59 const IY = 2 * N + 2;
60 const IYD = 2 * N + 3;
62 const state = new Float64Array(dim);
63 state[IYD] = cfg.mallet_velocity;
65 let lastFc = 0.0;
66 let lastXidd = 0.0;
68 const qddScratch = new Float64Array(N);
70 function deriv(s: Float64Array, d: Float64Array): void {
71 // Mallet contact
72 let w = 0.0;
73 let wd = 0.0;
74 for (let i = 0; i < N; i++) {
75 w += phiC[i] * s[i];
76 wd += phiC[i] * s[N + i];
77 }
78 const y = s[IY];
79 const yd = s[IYD];
80 const delta = y - w;
81 let fc = 0.0;
82 if (delta > 0.0) {
83 // Hunt-Crossley: f = K d^a (1 + lambda * ddot); lambda is in s/m so the
84 // correction stays dimensionless and O(1) over a sub-ms collision.
85 fc = K * Math.pow(delta, alpha) * (1.0 + mu * (yd - wd));
86 if (fc < 0.0) fc = 0.0;
87 }
88 lastFc = fc;
90 // Cavity pressure from combined volume change of walls and slit plug
91 let Aq = 0.0;
92 for (let i = 0; i < N; i++) Aq += A[i] * s[i];
93 const pc = pScale * (S * s[IX] + Aq);
95 for (let i = 0; i < N; i++) {
96 const qdd = -twoZetaOm[i] * s[N + i] - om2[i] * s[i] + fc * phiC[i] - A[i] * pc;
97 d[i] = s[N + i];
98 d[N + i] = qdd;
99 qddScratch[i] = qdd;
100 }
101 const xidd = -twoZetaHOm * s[IXD] - omegaH2 * s[IX] - cavDrive * Aq;
102 d[IX] = s[IXD];
103 d[IXD] = xidd;
104 d[IY] = yd;
105 d[IYD] = -fc * invMh;
106 lastXidd = xidd;
107 }
109 const Qdd: Float64Array[] = [];
110 for (let i = 0; i < N; i++) Qdd.push(new Float64Array(nSteps));
111 const Xdd = new Float64Array(nSteps);
112 const force = new Float64Array(nSteps);
114 const k1 = new Float64Array(dim);
115 const k2 = new Float64Array(dim);
116 const k3 = new Float64Array(dim);
117 const k4 = new Float64Array(dim);
118 const tmp = new Float64Array(dim);
120 const progressEvery = Math.max(1, Math.floor(nSteps / 50));
122 for (let k = 0; k < nSteps; k++) {
123 deriv(state, k1);
124 for (let i = 0; i < N; i++) Qdd[i][k] = qddScratch[i];
125 Xdd[k] = lastXidd;
126 force[k] = lastFc;
128 for (let i = 0; i < dim; i++) tmp[i] = state[i] + 0.5 * dt * k1[i];
129 deriv(tmp, k2);
130 for (let i = 0; i < dim; i++) tmp[i] = state[i] + 0.5 * dt * k2[i];
131 deriv(tmp, k3);
132 for (let i = 0; i < dim; i++) tmp[i] = state[i] + dt * k3[i];
133 deriv(tmp, k4);
134 for (let i = 0; i < dim; i++) {
135 state[i] += (dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
136 }
138 if (onProgress && k % progressEvery === 0) onProgress(k / nSteps);
139 }
141 let peakForce = 0.0;
142 let contactSamples = 0;
143 for (let k = 0; k < nSteps; k++) {
144 if (force[k] > peakForce) peakForce = force[k];
145 if (force[k] > 0.0) contactSamples++;
146 }
148 return {
149 Qdd,
150 Xdd,
151 force,
152 fsHi,
153 helmholtzHz: omegaH / (2 * Math.PI),
154 peakForceN: peakForce,
155 contactMs: (contactSamples * dt) * 1e3,
156 };
157}