17959acInteractive tutorial on the theory and the solverJeremy Magland 1// The contact law and what it does on impact. Left: f = K [xi]+^alpha, the
2// Hertz-type compression law of the mallet head. Right: the head thrown at a
3// small plate-like target, five modal oscillators sharing the contact point,
4// integrated through and past the contact. The same simulation renders a
5// second of the plate's response as audio, so mallet hardness can be heard,
6// not just plotted: a soft head excites only the low modes, a hard head all
7// of them.
9import { useMemo, useState } from 'react';
10import { useCanvas, drawFrame, PLOT_COLORS } from '../../components/plots/canvasUtils';
11import { Figure, FigControls, FigSlider, usePlayer } from '../figparts';
13const HEIGHT = 260;
14const M_HEAD = 0.028;
15// A plate-like ladder of point-coupled modes: frequency (Hz), modal mass
16// (kg), and amplitude decay rate (1/s).
17const MODES = [
18 { f: 180, m: 0.35, g: 2.5 },
19 { f: 476, m: 0.22, g: 3.5 },
20 { f: 894, m: 0.16, g: 5.0 },
21 { f: 1520, m: 0.12, g: 7.0 },
22 { f: 2350, m: 0.10, g: 10.0 },
23];
24const PLOT_T = 0.012;
25const AUDIO_SR = 25000;
26const DT = 1e-6; // 40 substeps per audio sample
28interface StrikeState {
29 x1: number; v1: number;
30 q: Float64Array; qd: Float64Array;
31}
33function makeState(v0: number): StrikeState {
34 return { x1: -1e-5, v1: v0, q: new Float64Array(MODES.length), qd: new Float64Array(MODES.length) };
35}
37// One symplectic-Euler substep; returns the contact force and the summed
38// modal acceleration (the collocated pickup).
39function step(s: StrikeState, K: number, alpha: number): { fc: number; acc: number } {
40 let ub = 0;
41 for (let i = 0; i < MODES.length; i++) ub += s.q[i];
42 const xi = s.x1 - ub;
43 const fc = xi > 0 ? K * Math.pow(xi, alpha) : 0;
44 s.v1 += (-fc / M_HEAD) * DT;
45 s.x1 += s.v1 * DT;
46 let acc = 0;
47 for (let i = 0; i < MODES.length; i++) {
48 const md = MODES[i];
49 const w = 2 * Math.PI * md.f;
50 const ai = (fc - md.m * w * w * s.q[i] - 2 * md.m * md.g * s.qd[i]) / md.m;
51 s.qd[i] += ai * DT;
52 s.q[i] += s.qd[i] * DT;
53 acc += ai;
54 }
55 return { fc, acc };
56}
58function simulate(K: number, alpha: number, v0: number) {
59 const s = makeState(v0);
60 const n = Math.round(PLOT_T / DT);
61 const keep = 600;
62 const stride = Math.max(1, Math.floor(n / keep));
63 const t: number[] = [];
64 const f: number[] = [];
65 let peakF = 0;
66 let peakXi = 0;
67 let contacts = 0;
68 let inC = false;
69 for (let k = 0; k < n; k++) {
70 let ub = 0;
71 for (let i = 0; i < MODES.length; i++) ub += s.q[i];
72 const xi = s.x1 - ub;
73 if (xi > peakXi) peakXi = xi;
74 const { fc } = step(s, K, alpha);
75 if (fc > 0 && !inC) contacts++;
76 inC = fc > 0;
77 if (fc > peakF) peakF = fc;
78 if (k % stride === 0) {
79 t.push(k * DT);
80 f.push(fc);
81 }
82 }
83 return { t, f, peakF, peakXi, contacts };
84}
86export function renderStrikeAudio(K: number, alpha: number, v0: number): Float32Array<ArrayBuffer> {
87 const dur = 1.2;
88 const nOut = Math.round(AUDIO_SR * dur);
89 const sub = Math.round(1 / (AUDIO_SR * DT));
90 const s = makeState(v0);
91 const y = new Float32Array(nOut);
92 for (let k = 0; k < nOut; k++) {
93 let acc = 0;
94 for (let j = 0; j < sub; j++) acc = step(s, K, alpha).acc;
95 y[k] = acc;
96 }
97 let pk = 0;
98 for (let k = 0; k < nOut; k++) pk = Math.max(pk, Math.abs(y[k]));
99 const fadeOut = Math.round(0.1 * AUDIO_SR);
100 for (let k = 0; k < nOut; k++) {
101 let g = 0.5 / (pk || 1);
102 if (k > nOut - fadeOut) g *= (nOut - k) / fadeOut;
103 y[k] *= g;
104 }
105 return y;
106}
108export default function HertzFig() {
109 const [alpha, setAlpha] = useState(2.54);
110 const [K, setK] = useState(1.6e8);
111 const [v0, setV0] = useState(5);
112 const sim = useMemo(() => simulate(K, alpha, v0), [K, alpha, v0]);
113 const play = usePlayer();
115 const ref = useCanvas(
116 (ctx, w, h) => {
117 drawFrame(ctx, w, h);
118 const pw = (w - 54) / 2;
119 const ph = h - 62;
120 const y0 = 14;
122 // ---- left: force vs compression ----
123 const xiMax = Math.max(1.5 * sim.peakXi, 1e-4);
124 const fAtXiMax = K * Math.pow(xiMax, alpha);
125 const x0 = 34;
126 ctx.strokeStyle = PLOT_COLORS.grid;
127 ctx.strokeRect(x0, y0, pw, ph);
128 ctx.strokeStyle = PLOT_COLORS.trace;
129 ctx.lineWidth = 2;
130 ctx.beginPath();
131 for (let i = 0; i <= 120; i++) {
132 const xi = (i / 120) * xiMax;
133 const fv = K * Math.pow(xi, alpha);
134 const x = x0 + (xi / xiMax) * pw;
135 const y = y0 + ph - (fv / fAtXiMax) * ph;
136 if (i === 0) ctx.moveTo(x, y);
137 else ctx.lineTo(x, y);
138 }
139 ctx.stroke();
140 ctx.lineWidth = 1;
141 // operating point
142 if (sim.peakXi > 0) {
143 ctx.fillStyle = PLOT_COLORS.accent;
144 ctx.beginPath();
145 ctx.arc(
146 x0 + (sim.peakXi / xiMax) * pw,
147 y0 + ph - (sim.peakF / fAtXiMax) * ph,
148 4, 0, 2 * Math.PI,
149 );
150 ctx.fill();
151 }
152 ctx.fillStyle = PLOT_COLORS.axis;
153 ctx.font = '11px system-ui';
154 ctx.textAlign = 'center';
155 ctx.fillText(
156 `f = K·ξ^${alpha.toFixed(2)} (compression 0…${(xiMax * 1000).toFixed(2)} mm)`,
157 x0 + pw / 2,
158 y0 + ph + 18,
159 );
160 ctx.textAlign = 'left';
161 ctx.fillText('force vs compression', x0 + 6, y0 + 14);
163 // ---- right: force vs time ----
164 const x1p = x0 + pw + 20;
165 const fMax = Math.max(1, sim.peakF) * 1.08;
166 ctx.strokeStyle = PLOT_COLORS.grid;
167 ctx.strokeRect(x1p, y0, pw, ph);
168 ctx.strokeStyle = PLOT_COLORS.accent;
169 ctx.lineWidth = 2;
170 ctx.beginPath();
171 for (let i = 0; i < sim.t.length; i++) {
172 const x = x1p + (sim.t[i] / PLOT_T) * pw;
173 const y = y0 + ph - (sim.f[i] / fMax) * ph;
174 if (i === 0) ctx.moveTo(x, y);
175 else ctx.lineTo(x, y);
176 }
177 ctx.stroke();
178 ctx.lineWidth = 1;
179 ctx.fillStyle = PLOT_COLORS.axis;
180 ctx.textAlign = 'center';
181 ctx.fillText(`force vs time (0…${(PLOT_T * 1000).toFixed(0)} ms)`, x1p + pw / 2, y0 + ph + 18);
182 ctx.textAlign = 'left';
183 ctx.fillText(
184 `peak ${sim.peakF.toFixed(0)} N, ${sim.contacts} contact${sim.contacts === 1 ? '' : 's'}`,
185 x1p + 6,
186 y0 + 14,
187 );
188 },
189 [sim, K, alpha],
190 HEIGHT,
191 );
193 return (
194 <Figure
195 caption={
196 <>
197 Left: the contact law, a power law in the compression with exponent α. Hertz
198 theory gives α = 1.5 for a hard sphere on a plate; the paper's measured felt
199 mallet is much more strongly hardening, α = 2.54, so the felt is soft at
200 grazing contact and stiffens as it packs. Right: the force the head actually
201 delivers when thrown at a small plate-like target (five modes sharing the contact
202 point). Soften the head and the pulse lengthens; stiffen it and the strike becomes
203 a brief spike; at some settings the target rebounds into the head and the pulse
204 splits, the double peak Samejima's Fig. 12 shows for the full stick model. The play
205 button renders a second of the target's response to this exact strike: a long soft
206 pulse cannot excite the upper modes (the thud of felt), a short hard one excites
207 them all (the ping of a wood tip). The excitation spectrum is the pulse's Fourier
208 transform, and you are hearing its width change.
209 </>
210 }
211 >
212 <canvas ref={ref} />
213 <FigControls>
214 <FigSlider
215 label="exponent α"
216 value={alpha}
217 min={1.2}
218 max={3.2}
219 step={0.02}
220 format={(v) => v.toFixed(2)}
221 onChange={setAlpha}
222 />
223 <FigSlider
224 label="stiffness K"
225 value={K}
226 min={1e6}
227 max={1e10}
228 log
229 format={(v) => v.toExponential(1)}
230 onChange={setK}
231 />
232 <FigSlider
233 label="impact speed"
234 value={v0}
235 min={0.5}
236 max={8}
237 step={0.1}
238 format={(v) => `${v.toFixed(1)} m/s`}
239 onChange={setV0}
240 />
241 <button onClick={() => play(renderStrikeAudio(K, alpha, v0), AUDIO_SR)}>
242 ▶ Hear this strike
243 </button>
244 </FigControls>
245 </Figure>
246 );
247}