982850fWoodblock strike simulator: Python reference model and interactive web appJeremy Magland 1// Modal basis: free-free Euler-Bernoulli bar, mass-normalized bending modes.
2// Port of ModalBasis in woodblock.py.
4import type { WoodblockParams } from '../params';
6// Roots of cos(bL)cosh(bL) = 1 for the free-free beam.
7const BETA_L = [4.73004074, 7.85320462, 10.99560784, 14.13716549];
9export function betaL(n: number): number {
10 if (n <= BETA_L.length) return BETA_L[n - 1];
11 return ((2 * n + 1) * Math.PI) / 2.0;
12}
14// Free-free mode shape on s in [0,1], unnormalized.
15export function modeShape(bl: number, s: Float64Array): Float64Array {
16 const sigma = (Math.cosh(bl) - Math.cos(bl)) / (Math.sinh(bl) - Math.sin(bl));
17 const out = new Float64Array(s.length);
18 for (let i = 0; i < s.length; i++) {
19 const x = bl * s[i];
20 out[i] = Math.cosh(x) + Math.cos(x) - sigma * (Math.sinh(x) + Math.sin(x));
21 }
22 return out;
23}
25function trapz(y: Float64Array, ds: number): number {
26 let sum = 0;
27 for (let i = 0; i < y.length; i++) sum += y[i];
28 sum -= 0.5 * (y[0] + y[y.length - 1]);
29 return sum * ds;
30}
32function interp1(x: number, xs: Float64Array, ys: Float64Array): number {
33 if (x <= xs[0]) return ys[0];
34 const n = xs.length;
35 if (x >= xs[n - 1]) return ys[n - 1];
36 // uniform grid
37 const dx = xs[1] - xs[0];
38 const j = Math.min(n - 2, Math.floor((x - xs[0]) / dx));
39 const t = (x - xs[j]) / dx;
40 return ys[j] * (1 - t) + ys[j + 1] * t;
41}
43export interface ModalBasisData {
44 n: number;
45 f: Float64Array; // modal frequencies (Hz)
46 omega: Float64Array;
47 zeta: Float64Array;
48 phiC: Float64Array; // mode shape at strike point (1/sqrt(kg))
49 A: Float64Array; // cavity volume coupling (m^2/sqrt(kg))
50 mass: number;
51 // Mass-normalized shapes sampled on the dense s grid, for later
52 // interpolation at arbitrary radiating-element positions.
53 sGrid: Float64Array;
54 shapes: Float64Array[]; // length n, each sampled on sGrid
55}
57const NS = 4096;
59export function buildBasis(cfg: WoodblockParams): ModalBasisData {
60 const L = cfg.length_m;
61 const W = cfg.width_m;
62 const H = cfg.wall_thickness_m;
63 const area = W * H;
64 const mass = cfg.participating_mass_kg || cfg.density * area * L;
65 const cLong = Math.sqrt(cfg.youngs_modulus / cfg.density);
67 // f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
68 const prefactor = (cLong * H) / (2.0 * Math.PI * Math.sqrt(12.0) * L * L);
70 const s = new Float64Array(NS);
71 for (let i = 0; i < NS; i++) s[i] = i / (NS - 1);
72 const ds = 1.0 / (NS - 1);
74 const freqs: number[] = [];
75 const phiC: number[] = [];
76 const A: number[] = [];
77 const shapes: Float64Array[] = [];
78 const lo = cfg.cavity_span_lo;
79 const hi = cfg.cavity_span_hi;
81 for (let n = 1; n <= cfg.n_modes; n++) {
82 const bl = betaL(n);
83 let f = prefactor * bl * bl * cfg.hollow_softening;
84 if (n > 1 && cfg.inharmonicity !== 1.0) {
85 const f0 = prefactor * betaL(1) ** 2 * cfg.hollow_softening;
86 f = f0 * (f / f0) ** cfg.inharmonicity;
87 }
88 if (f > cfg.max_mode_hz) break;
90 const raw = modeShape(bl, s);
91 const sq = new Float64Array(NS);
92 for (let i = 0; i < NS; i++) sq[i] = raw[i] * raw[i];
93 const rms = Math.sqrt(trapz(sq, ds));
94 const norm = 1.0 / (rms * Math.sqrt(mass));
95 const shape = new Float64Array(NS);
96 for (let i = 0; i < NS; i++) shape[i] = raw[i] * norm;
98 freqs.push(f);
99 phiC.push(interp1(cfg.strike_position, s, shape));
100 const masked = new Float64Array(NS);
101 for (let i = 0; i < NS; i++) masked[i] = s[i] >= lo && s[i] <= hi ? shape[i] : 0.0;
102 A.push(W * L * trapz(masked, ds));
103 shapes.push(shape);
104 }
106 const nKept = freqs.length;
107 const fArr = Float64Array.from(freqs);
108 const omega = new Float64Array(nKept);
109 const zeta = new Float64Array(nKept);
110 const AArr = new Float64Array(nKept);
111 for (let i = 0; i < nKept; i++) {
112 omega[i] = 2 * Math.PI * fArr[i];
113 zeta[i] = 0.5 * (cfg.loss_factor + cfg.loss_slope * omega[i]);
114 AArr[i] = A[i] * cfg.cavity_coupling;
115 }
117 return {
118 n: nKept,
119 f: fArr,
120 omega,
121 zeta,
122 phiC: Float64Array.from(phiC),
123 A: AArr,
124 mass,
125 sGrid: s,
126 shapes,
127 };
128}
130// Interpolate the stored mode shapes at M radiating-element centres.
131// Returns phiElem as an array of n Float64Arrays of length M, plus sElem.
132export function elementShapes(basis: ModalBasisData, M: number): {
133 phiElem: Float64Array[];
134 sElem: Float64Array;
135} {
136 const sElem = new Float64Array(M);
137 for (let j = 0; j < M; j++) sElem[j] = (j + 0.5) / M;
138 const phiElem: Float64Array[] = [];
139 for (let n = 0; n < basis.n; n++) {
140 const row = new Float64Array(M);
141 for (let j = 0; j < M; j++) row[j] = interp1(sElem[j], basis.sGrid, basis.shapes[n]);
142 phiElem.push(row);
143 }
144 return { phiElem, sElem };
145}