2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 1// The losses derived from physical constants: Zener's thermoelastic loss
2// factor has its peak where it should; the outline perimeter of a disc and
3// a rectangle; the limits of the average radiation efficiency; the fit of
4// the two-parameter operator; the per-mode radiation rate against the
5// closed form for a baffled piston and against the direct double sum; the
6// per-mode thermoelastic rate of a uniform plate against eta_te omega / 2.
8import { test } from "node:test";
9import assert from "node:assert/strict";
10import { makeGrid, withMask, nodeWeight, RHO_AIR, SPEED_OF_SOUND, ROOM_TEMPERATURE, type Material } from "../src/params.ts";
11import { uniformPlate } from "../src/model/reference.ts";
12import { assemble } from "../src/model/assemble.ts";
13import { computeModes, localRunner } from "../src/modal/modes.ts";
14import { modalModel, modalRates, sliceModel } from "../src/modal/synth.ts";
15import { dampingPlate, radiationRates, radiationRateDirect, thermoelasticRates, computeDamping, bendingEnergyPerNode } from "../src/modal/damping.ts";
16import { thermoelasticRelaxation, thermoelasticLossFactor, outlinePerimeter, averageRadiationEfficiency, coincidenceFrequency, fitLine, fitSigma, fittedRate, plateGeometry, estimatedRate, materialRate, LN1000 } from "../src/loss.ts";
18const steel: Material = { E: 200e9, nu: 0.3, rho: 7850, eta: 3e-4, thermal: { alpha: 12e-6, cp: 460, k: 50 } };
19const bronze: Material = { E: 110e9, nu: 0.34, rho: 8800, eta: 2e-4, thermal: { alpha: 18e-6, cp: 370, k: 40 } };
21test("Zener's thermoelastic loss peaks at omega tau = 1 with the value delta / 2", () => {
22 const h = 1e-3;
23 const { delta, tau } = thermoelasticRelaxation(bronze, h, ROOM_TEMPERATURE);
24 const th = bronze.thermal!;
25 assert.ok(Math.abs(tau - (h * h * bronze.rho * th.cp) / (Math.PI * Math.PI * th.k)) < 1e-15);
26 assert.ok(Math.abs(delta - ((bronze.E * th.alpha ** 2 * ROOM_TEMPERATURE) / (bronze.rho * th.cp)) * ((1 + bronze.nu) / (1 - bronze.nu))) < 1e-15);
27 const peak = thermoelasticLossFactor(bronze, h, 1 / tau, ROOM_TEMPERATURE);
28 assert.ok(Math.abs(peak - delta / 2) < 1e-15);
29 assert.ok(thermoelasticLossFactor(bronze, h, 0.1 / tau, ROOM_TEMPERATURE) < peak);
30 assert.ok(thermoelasticLossFactor(bronze, h, 10 / tau, ROOM_TEMPERATURE) < peak);
31 // For a millimeter of bronze the peak lies below 100 Hz and the factor is of order 1e-3.
32 assert.ok(1 / (2 * Math.PI * tau) < 100 && delta / 2 > 1e-4 && delta / 2 < 1e-2);
33 // Twice the thickness quarters the peak frequency; without thermal constants there is no loss.
34 assert.ok(Math.abs(thermoelasticRelaxation(bronze, 2 * h, ROOM_TEMPERATURE).tau / tau - 4) < 1e-12);
35 assert.equal(thermoelasticLossFactor({ ...bronze, thermal: null }, h, 1000, ROOM_TEMPERATURE), 0);
36});
38test("the outline perimeter of a rectangle and of a disc", () => {
39 const rect = makeGrid(0.3, 0.2, 0.01);
40 // The contour cuts each corner diagonally, which costs about 0.3 cells per corner.
41 assert.ok(Math.abs(outlinePerimeter(rect) / (2 * (0.3 + 0.2)) - 1) < 0.02, `rectangle perimeter ${outlinePerimeter(rect)}`);
42 const R = 0.2;
43 const disc = withMask(makeGrid(2 * R, 2 * R, 0.008), (x, y) => x * x + y * y <= R * R);
44 // The contour's segments run in eight directions only, which overstates the length of a curve by a few percent (a staircase would overstate it by 27 percent).
45 assert.ok(Math.abs(outlinePerimeter(disc) / (2 * Math.PI * R) - 1) < 0.08, `disc perimeter ${outlinePerimeter(disc)}`);
46 assert.ok(outlinePerimeter(disc) > 2 * Math.PI * R, "an overestimate, never an underestimate");
47 const ring = withMask(makeGrid(2 * R, 2 * R, 0.008), (x, y) => x * x + y * y <= R * R && x * x + y * y >= 0.05 ** 2);
48 assert.ok(Math.abs(outlinePerimeter(ring) / (2 * Math.PI * (R + 0.05)) - 1) < 0.1, "a hole adds its own perimeter");
49});
51test("the average radiation efficiency has the textbook limits", () => {
52 const fc = coincidenceFrequency(bronze, 1e-3);
53 assert.ok(fc > 15000 && fc < 20000, `coincidence ${fc} Hz for 1 mm bronze`);
54 const A = 0.13, U = 1.28, a = 0.4, b = 0.4;
55 // Simply supported: about linear in frequency far below coincidence.
56 const s1 = averageRadiationEfficiency(100, fc, A, U, a, b, false), s2 = averageRadiationEfficiency(200, fc, A, U, a, b, false);
57 assert.ok(s1 > 0 && s1 < 0.1);
58 assert.ok(Math.abs(s2 / s1 - 2) < 0.1, "grows about linearly with frequency far below coincidence");
59 assert.ok(averageRadiationEfficiency(fc / 2, fc, A, U, a, b, false) > s2);
60 assert.ok(Math.abs(averageRadiationEfficiency(4 * fc, fc, A, U, a, b) - 1 / Math.sqrt(0.75)) < 1e-9, "above coincidence");
61 assert.ok(averageRadiationEfficiency(fc * 1.0001, fc, A, U, a, b) <= (0.5 - 0.15) * Math.sqrt(a / (SPEED_OF_SOUND / fc)) + 1e-9, "capped at the coincidence value");
62 // Free edges: the same times (f / fc)^2 below coincidence, continuous at coincidence.
63 const f1 = averageRadiationEfficiency(100, fc, A, U, a, b), f2 = averageRadiationEfficiency(200, fc, A, U, a, b);
64 assert.ok(Math.abs(f1 / s1 - (100 / fc) ** 2) < 1e-12 && Math.abs(f2 / f1 - 8) < 0.4, "free edges: about cubic in frequency far below coincidence");
65 assert.ok(Math.abs(averageRadiationEfficiency(fc * 0.9999, fc, A, U, a, b) / averageRadiationEfficiency(fc * 0.9999, fc, A, U, a, b, false) - 1) < 1e-3);
66});
68test("the two-parameter fit reproduces a line and bounds the estimate", () => {
69 const kappa = 1.1;
70 const x = [1, 2, 3, 5, 8], y = x.map((v) => 0.7 + 2e-4 * v);
71 const fit = fitLine(x, y, kappa);
72 assert.ok(Math.abs(fit.sigma0 - 0.7) < 1e-9 && Math.abs(fit.sigma1 - 2e-4) < 1e-12);
73 assert.ok(fitLine(x, x.map((v) => 3 - v), kappa).sigma1 === 0, "a falling curve gives sigma1 = 0");
74 const g = makeGrid(0.4, 0.4, 0.02);
75 const pl = uniformPlate(g, bronze, 1e-3);
76 const geom = plateGeometry(pl);
77 const f2 = fitSigma(pl, geom);
78 assert.ok(f2.sigma0 >= 0 && f2.sigma1 >= 0);
79 const e100 = estimatedRate(pl, geom, 100), e5k = estimatedRate(pl, geom, 5000);
80 assert.ok(e5k.total > e100.total, "the estimate rises with frequency");
81 assert.ok(Math.abs(e100.material - materialRate(bronze.eta, 2 * Math.PI * 100)) < 1e-12);
82 assert.ok(e100.radiation > 0 && e100.thermoelastic > 0);
83 assert.ok(fittedRate(f2, 1000) > 0.3 * e100.total && fittedRate(f2, 1000) < 3 * e5k.total);
84});
86/** A small free steel plate with its modes and per-mode losses. */
87async function plateWithModes(mat: Material, spacing = 0.0125) {
88 const g = makeGrid(0.2, 0.15, spacing);
89 const pl = uniformPlate(g, mat, 1e-3);
90 const sys = assemble(pl);
91 const modes = await computeModes(sys, 4000, localRunner(sys), { sliceSize: 20 });
92 const model = modalModel(pl, sys, modes);
93 return { pl, model, dp: dampingPlate(pl) };
94}
96test("the radiation rate of a uniform piston at low frequency equals the closed form", async () => {
97 const { model, dp } = await plateWithModes(steel);
98 // One fake mode of uniform shape 1 / sqrt(M) at 20 Hz, where the plate is small compared with the wavelength.
99 const one = sliceModel(model, 0, 1);
100 let M = 0, S = 0;
101 for (let k = 0; k < one.nGrid; k++) { M += one.masses[k]; S += dp.area[k]; }
102 for (let r = 0; r < one.nRows; r++) one.wShapes[r] = 1 / Math.sqrt(M);
103 one.freqs = new Float64Array([20]);
104 const omega = 2 * Math.PI * 20;
105 const expected = (RHO_AIR * omega * omega * S * S) / (2 * Math.PI * SPEED_OF_SOUND * M); // both faces of a baffled piston
106 const fftRate = radiationRates(one, dp, 0, 1)[0];
107 assert.ok(Math.abs(fftRate / expected - 1) < 2e-3, `fft ${fftRate} vs ${expected}`);
108 assert.ok(Math.abs(radiationRateDirect(one, dp, 0) / expected - 1) < 2e-3);
109});
111test("the FFT convolution equals the direct double sum for every mode", async () => {
112 const { model, dp } = await plateWithModes(steel);
113 const rates = radiationRates(model, dp, 0, model.count);
114 for (let m = 0; m < model.count; m += Math.max(1, Math.floor(model.count / 12))) {
115 const direct = radiationRateDirect(model, dp, m);
116 assert.ok(Math.abs(rates[m] - direct) < 1e-9 * Math.max(1, direct), `mode ${m}: ${rates[m]} vs ${direct}`);
117 }
118 // The high modes radiate faster than the low ones in absolute terms, and the rates are positive.
119 for (let m = 0; m < model.count; m++) assert.ok(rates[m] >= 0);
120 assert.ok(rates[model.count - 1] > rates[0]);
121});
123test("simply supported shapes reproduce Wallace's low-frequency radiation efficiencies", async () => {
124 // sin(m pi x / a) sin(n pi y / b) on a baffled rectangle at k a << 1 (Wallace 1972):
125 // odd-odd 32 k^2 ab / (pi^5 m^2 n^2), odd-even 8 k^4 a b^3 / (3 pi^5 m^2 n^2), even-even 2 k^6 a^3 b^3 / (15 pi^5 m^2 n^2).
126 const a = 0.3, b = 0.45, h = 1e-3, f = 30, c = SPEED_OF_SOUND, k = (2 * Math.PI * f) / c;
127 const g = makeGrid(a, b, 0.01), pl = uniformPlate(g, steel, h), dp = dampingPlate(pl);
128 const rhoH = steel.rho * h;
129 for (const [m, n] of [[1, 1], [1, 2], [2, 1], [2, 2], [3, 3]]) {
130 const rows: number[] = [], vals: number[] = [];
131 let norm = 0;
132 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
133 const kk = j * g.nx + i, w = nodeWeight(g, i, j);
134 if (!w) continue;
135 const v = Math.sin((m * Math.PI * i * g.hx) / a) * Math.sin((n * Math.PI * j * g.hy) / b);
136 rows.push(kk); vals.push(v); norm += rhoH * g.hx * g.hy * w * v * v;
137 }
2481a63Remove the standJeremy Magland 138 const one = { nGrid: g.n, nRows: rows.length, count: 1, freqs: new Float64Array([f]), wShapes: Float32Array.from(vals.map((v) => v / Math.sqrt(norm))), wRow: new Int32Array(0), rowNode: Int32Array.from(rows), sigmaMat: new Float64Array(1), sigmaTe: new Float64Array(1), sigmaRad: new Float64Array(1), masses: new Float64Array(0) };
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 139 const efficiency = (radiationRateDirect(one, dp, 0) * rhoH) / (RHO_AIR * c); // sigma = rho_air c sigma_bar / (rho h)
140 const odd = (q: number) => q % 2 === 1;
141 const wallace = odd(m) && odd(n) ? (32 * k * k * a * b) / (Math.PI ** 5 * m * m * n * n)
142 : !odd(m) && !odd(n) ? (2 * k ** 6 * a ** 3 * b ** 3) / (15 * Math.PI ** 5 * m * m * n * n)
143 : odd(m) ? (8 * k ** 4 * a * b ** 3) / (3 * Math.PI ** 5 * m * m * n * n) : (8 * k ** 4 * a ** 3 * b) / (3 * Math.PI ** 5 * m * m * n * n);
144 // The 1 cm grid samples the mean of the (3,3) shape about 3 percent low; the 8 mm grid is within 2 percent.
145 assert.ok(Math.abs(efficiency / wallace - 1) < 0.04, `(${m},${n}): ${efficiency} vs ${wallace}`);
146 }
147});
149test("the thermoelastic rate of a flat uniform plate is eta_te omega / 2, and the total rate is the sum", async () => {
150 const { pl, model, dp } = await plateWithModes(steel);
151 const rates = thermoelasticRates(model, dp, 0, model.count);
152 for (let m = 0; m < model.count; m++) {
153 const omega = 2 * Math.PI * model.freqs[m];
154 const expected = (thermoelasticLossFactor(steel, 1e-3, omega, ROOM_TEMPERATURE) * omega) / 2;
155 assert.ok(Math.abs(rates[m] / expected - 1) < 1e-4, `mode ${m}: ${rates[m]} vs ${expected}`);
156 }
157 // The bending energy of a mass-normalized mode of a flat free plate is omega^2 / 2.
158 const w = new Float64Array(dp.nx * dp.ny);
159 for (let r = 0; r < model.nRows; r++) w[model.rowNode[r]] = model.wShapes[r * model.count + 3];
160 const U = bendingEnergyPerNode(dp, w);
161 let total = 0;
162 for (const u of U) total += u;
163 const omega3 = 2 * Math.PI * model.freqs[3];
164 assert.ok(Math.abs(total / (omega3 * omega3 / 2) - 1) < 1e-4, `bending energy ${total} vs ${omega3 * omega3 / 2}`);
165 computeDamping(model, dp);
166 const sigma = modalRates(model);
167 for (let m = 0; m < model.count; m++) {
168 const omega = 2 * Math.PI * model.freqs[m];
169 assert.ok(Math.abs(sigma[m] - (materialRate(steel.eta, omega) + model.sigmaTe[m] + model.sigmaRad[m])) < 1e-12);
170 assert.ok(model.sigmaMat[m] > 0 && model.sigmaRad[m] > 0 && model.sigmaTe[m] > 0);
171 assert.ok(LN1000 / sigma[m] > 0.1 && LN1000 / sigma[m] < 1000, `T60 ${LN1000 / sigma[m]} s`);
172 }
173 void pl;
174});