2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 1// The losses of the plate, derived from physical constants rather than
2// prescribed. Three mechanisms drain the energy of a vibrating plate in this
3// model. Internal friction of the metal is a hysteretic loss with the
4// material's loss factor eta: a mode of angular frequency omega decays at the
5// rate eta omega / 2. Thermoelastic damping (Zener) comes from the heat that
6// flows across the thickness between the compressed and the stretched face
7// of a bent plate; it follows from the thermal constants of the metal and
8// the local thickness, with a peak at the frequency set by the thermal
9// diffusion time across the thickness. Radiation into the air is the energy
10// the plate loses as sound; the linear model computes it exactly for each
11// mode from the mode's shape (modal/damping.ts), and this module gives the
12// analytic average of Leppington for a plate of given area and perimeter,
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 13
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 14import { type Grid, type Material, type Plate, RHO_AIR, SPEED_OF_SOUND, bendingConstant, insideCell, massProperties, meanThickness, nodeWeight } from "./params.ts";
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 15
16/** ln(1000): a rate sigma [1/s] and a decay time relate by sigma = LN1000 / T60. */
17export const LN1000 = Math.log(1000);
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 19/** A decay time in seconds for display: an infinity sign for no decay, two decimals below 10 s, whole seconds below 1000 s, a power of ten beyond (a loss that is effectively absent). */
20export function t60Text(seconds: number): string {
21 if (!isFinite(seconds) || seconds <= 0) return "∞";
22 if (seconds < 10) return seconds.toFixed(2);
23 if (seconds < 1000) return seconds.toFixed(0);
24 return seconds.toExponential(0).replace("e+", "e");
25}
27/** The decay time of a rate for display (see t60Text). */
28export function rateText(sigma: number): string {
29 return t60Text(sigma > 0 ? LN1000 / sigma : Infinity);
30}
32// ---------- Internal friction ----------
34/** Decay rate of a mode of angular frequency omega from a hysteretic loss factor eta: the energy lost per cycle is 2 pi eta times the stored energy. */
35export function materialRate(eta: number, omega: number): number {
36 return (eta * omega) / 2;
37}
39// ---------- Thermoelastic damping ----------
41/**
42 * Zener's thermoelastic relaxation of a plate of thickness h: the relaxation
43 * strength delta and the relaxation time tau. The strength carries the
44 * factor (1 + nu) / (1 - nu) of the biaxial stress state of a bent plate
45 * (Zener's beam value is E alpha^2 T / (rho cp)); the time is the thermal
46 * diffusion time across the thickness, h^2 / (pi^2 chi) with the diffusivity
47 * chi = k / (rho cp).
48 */
49export function thermoelasticRelaxation(mat: Material, h: number, T: number): { delta: number; tau: number } {
50 const th = mat.thermal;
51 if (!th) return { delta: 0, tau: 1 };
52 const delta = ((mat.E * th.alpha * th.alpha * T) / (mat.rho * th.cp)) * ((1 + mat.nu) / (1 - mat.nu));
53 const tau = (h * h * mat.rho * th.cp) / (Math.PI * Math.PI * th.k);
54 return { delta, tau };
55}
57/** Thermoelastic loss factor of bending at angular frequency omega in a plate of thickness h: delta omega tau / (1 + (omega tau)^2), largest (delta / 2) at omega tau = 1. */
58export function thermoelasticLossFactor(mat: Material, h: number, omega: number, T: number): number {
59 const { delta, tau } = thermoelasticRelaxation(mat, h, T);
60 const x = omega * tau;
61 return (delta * x) / (1 + x * x);
62}
64// ---------- Geometry of the outline ----------
66/**
67 * The length of the outline of the plate, from the marching-squares contour
68 * of the cell mask at level one half (cells outside the grid count as
69 * outside). A staircase count would overstate the perimeter of a disc by
70 * 4 / pi; the contour, whose segments run in eight directions, overstates
71 * it by a few percent, which is well within the accuracy of the average
72 * radiation efficiency it feeds.
73 */
74export function outlinePerimeter(g: Grid): number {
75 const cx = g.nx - 1, cy = g.ny - 1;
76 const inside = (ci: number, cj: number) => ci >= 0 && cj >= 0 && ci < cx && cj < cy && insideCell(g, ci, cj);
77 const hx = g.hx, hy = g.hy;
78 let length = 0;
79 // Squares between the centers of cells (ci, cj), (ci+1, cj), (ci+1, cj+1), (ci, cj+1), one ring beyond the grid.
80 for (let cj = -1; cj < cy; cj++) {
81 for (let ci = -1; ci < cx; ci++) {
82 const a = inside(ci, cj) ? 1 : 0, b = inside(ci + 1, cj) ? 1 : 0, c = inside(ci + 1, cj + 1) ? 1 : 0, d = inside(ci, cj + 1) ? 1 : 0;
83 const code = a | (b << 1) | (c << 2) | (d << 3);
84 if (code === 0 || code === 15) continue;
85 // Midpoints of the square's edges: bottom (a-b), right (b-c), top (c-d), left (d-a), as offsets in cells.
86 const bottom = [0.5, 0], right = [1, 0.5], top = [0.5, 1], left = [0, 0.5];
87 const seg = (p: number[], q: number[]) => { length += Math.hypot((q[0] - p[0]) * hx, (q[1] - p[1]) * hy); };
88 switch (code) {
89 case 1: case 14: seg(left, bottom); break;
90 case 2: case 13: seg(bottom, right); break;
91 case 3: case 12: seg(left, right); break;
92 case 4: case 11: seg(right, top); break;
93 case 6: case 9: seg(bottom, top); break;
94 case 7: case 8: seg(left, top); break;
95 case 5: seg(left, bottom); seg(right, top); break; // saddle: two corners
96 case 10: seg(bottom, right); seg(left, top); break;
97 }
98 }
99 }
100 return length;
101}
103// ---------- Radiation ----------
105/** Coincidence frequency of bending waves in a plate of thickness h: where the bending wavelength equals the wavelength of sound, f_c = c^2 / (2 pi kappa). */
106export function coincidenceFrequency(mat: Material, h: number): number {
107 return (SPEED_OF_SOUND * SPEED_OF_SOUND) / (2 * Math.PI * bendingConstant(mat, h));
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 108}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 110/**
111 * The radiation efficiency of a baffled plate averaged over its modes near
112 * the frequency f. The simply supported case is that of Leppington,
113 * Broadbent and Heron (1982), as given by Fahy and Gardonio: below
114 * coincidence the edges radiate and the efficiency grows about linearly with
115 * frequency, at coincidence it peaks at a value set by the size of the plate
116 * in coincidence wavelengths, and above it tends to one. A plate with free
117 * edges radiates much less below coincidence: its modes are orthogonal to
118 * the rigid translation and rotations, so they have no monopole or dipole
119 * moment, and at a free edge the bending moment and the shear vanish, so
120 * the boundary terms that make a simply supported edge radiate are smaller
121 * by the factor (k / k_b)^2 in amplitude, k the acoustic and k_b the bending
122 * wavenumber, that is by (f / f_c)^2 in power. The per-mode computation of
123 * modal/damping.ts on the presets follows the simply supported curve times
124 * this factor within a factor of two from 100 Hz to 8 kHz, and that is what
125 * `freeEdges` applies. `a` and `b` are the shorter and the longer side.
126 */
127export function averageRadiationEfficiency(f: number, fc: number, area: number, perimeter: number, a: number, b: number, freeEdges = true): number {
128 const c = SPEED_OF_SOUND;
129 const lambdaC = c / fc;
130 const cap = Math.max(1, (0.5 - 0.15 * (a / b)) * Math.sqrt(a / lambdaC));
131 if (f >= fc) return Math.min(cap, 1 / Math.sqrt(Math.max(1e-12, 1 - fc / f)));
132 const mu = Math.sqrt(fc / f);
133 const s = ((perimeter * c) / (4 * Math.PI * Math.PI * area * fc)) * (1 / Math.sqrt(mu * mu - 1)) * (Math.log((mu + 1) / (mu - 1)) + (2 * mu) / (mu * mu - 1));
134 return Math.min(cap, s) * (freeEdges ? (f / fc) ** 2 : 1);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 135}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 137/** Decay rate from a radiation efficiency: the two faces of a baffled plate of surface mass rho h radiate the power 2 rho_air c sigma_rad <v^2> A, so sigma = rho_air c sigma_rad / (rho h). */
138export function radiationRateFromEfficiency(efficiency: number, mat: Material, h: number): number {
139 return (RHO_AIR * SPEED_OF_SOUND * efficiency) / (mat.rho * h);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 140}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 142// ---------- The estimate for a whole plate ----------
144/** Geometric quantities of a plate that the loss estimate needs. */
145export interface PlateGeometry {
146 area: number;
147 perimeter: number;
148 hMean: number;
149 kappa: number;
150 fc: number;
151 /** Shorter and longer side of the bounding box [m]. */
152 a: number;
153 b: number;
154 /** Node areas over the grid [m^2] (zero outside the plate), for the area-weighted thermoelastic loss. */
155 nodeArea: Float64Array;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 156}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 158export function plateGeometry(pl: Plate): PlateGeometry {
159 const g = pl.grid;
160 const { area } = massProperties(pl);
161 const hMean = meanThickness(pl);
162 const nodeArea = new Float64Array(g.n);
163 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) nodeArea[j * g.nx + i] = g.hx * g.hy * nodeWeight(g, i, j);
164 return { area, perimeter: outlinePerimeter(g), hMean, kappa: bendingConstant(pl.mat, hMean), fc: coincidenceFrequency(pl.mat, hMean), a: Math.min(pl.Lx, pl.Ly), b: Math.max(pl.Lx, pl.Ly), nodeArea };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 165}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 167/** The three contributions to the decay rate of a bending partial at frequency f [1/s], and their sum. */
168export interface LossEstimate {
169 material: number;
170 thermoelastic: number;
171 radiation: number;
172 total: number;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 173}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 175/**
176 * The decay rate of a bending partial at frequency f estimated without the
177 * modes: the material rate, the thermoelastic rate averaged over the plate
178 * area with the local thickness, and the radiation rate from the average
179 * efficiency at the mean thickness. The modal computation replaces the last
180 * two by exact per-mode values; this estimate serves the full model's loss
181 * operator and the summary of a design before its modes exist.
182 */
183export function estimatedRate(pl: Plate, geom: PlateGeometry, f: number): LossEstimate {
184 const omega = 2 * Math.PI * f;
185 const material = materialRate(pl.mat.eta, omega);
186 let te = 0;
187 if (pl.mat.thermal) {
188 let num = 0, den = 0;
189 for (let k = 0; k < geom.nodeArea.length; k++) {
190 const a = geom.nodeArea[k];
191 if (a === 0) continue;
192 num += a * thermoelasticLossFactor(pl.mat, pl.h[k], omega, pl.temperature);
193 den += a;
194 }
195 te = den > 0 ? ((num / den) * omega) / 2 : 0;
196 }
197 const radiation = radiationRateFromEfficiency(averageRadiationEfficiency(f, geom.fc, geom.area, geom.perimeter, geom.a, geom.b), pl.mat, geom.hMean);
198 return { material, thermoelastic: te, radiation, total: material + te + radiation };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 199}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 201/** The in-plane motion of the full model has no radiation and no thermoelastic loss in this model; it decays by internal friction at the rate of the lowest in-plane resonance, c_p / (2 L). */
202export function inPlaneRate(pl: Plate): number {
203 const cp = Math.sqrt(pl.mat.E / (pl.mat.rho * (1 - pl.mat.nu * pl.mat.nu)));
204 const fIp = cp / (2 * Math.max(pl.Lx, pl.Ly));
205 return materialRate(pl.mat.eta, 2 * Math.PI * fIp);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 206}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 208// ---------- The two-parameter operator of the full model ----------
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 210/**
211 * The two-parameter loss model of the full simulation. A bending mode of
212 * angular frequency omega decays at the rate sigma0 + sigma1 k^2 with
213 * k^2 = omega / kappa; the fit is a non-negative least-squares fit of that
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 214 * line to the estimated rates on log-spaced frequencies from 50 Hz to
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 215 * 10 kHz, at the mean thickness of the plate.
216 */
217export interface SigmaFit {
218 sigma0: number; // [1/s]
219 sigma1: number; // [m^2/s]
220 kappa: number; // [m^2/s] bending constant used for the fit
221}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 223export function fitSigma(pl: Plate, geom: PlateGeometry = plateGeometry(pl), fLo = 50, fHi = 10000, n = 40): SigmaFit {
224 const kappa = geom.kappa;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 225 const x: number[] = [], y: number[] = [];
226 for (let i = 0; i < n; i++) {
227 const f = fLo * Math.pow(fHi / fLo, i / (n - 1));
228 x.push((2 * Math.PI * f) / kappa); // k^2
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 229 y.push(estimatedRate(pl, geom, f).total);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 230 }
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 231 return fitLine(x, y, kappa);
232}
234/** Non-negative least squares of y = s0 + s1 x: with two parameters the solution is the unconstrained one when it is feasible, else the best of the fits with one parameter at zero, else zero. */
235export function fitLine(x: number[], y: number[], kappa: number): SigmaFit {
236 const n = x.length;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 237 let sx = 0, sy = 0, sxx = 0, sxy = 0;
238 for (let i = 0; i < n; i++) { sx += x[i]; sy += y[i]; sxx += x[i] * x[i]; sxy += x[i] * y[i]; }
239 const det = n * sxx - sx * sx;
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 240 const candidates: [number, number][] = [[0, 0]];
241 if (det > 0) { const s1 = (n * sxy - sx * sy) / det; candidates.push([(sy - s1 * sx) / n, s1]); }
242 candidates.push([sy / n, 0]);
243 if (sxx > 0) candidates.push([0, sxy / sxx]);
244 let best: [number, number] = [0, 0], bestRes = Infinity;
245 for (const [s0, s1] of candidates) {
246 if (!(s0 >= 0 && s1 >= 0)) continue;
247 let res = 0;
248 for (let i = 0; i < n; i++) res += (y[i] - s0 - s1 * x[i]) ** 2;
249 if (res < bestRes) { bestRes = res; best = [s0, s1]; }
250 }
251 return { sigma0: best[0], sigma1: best[1], kappa };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 252}
254/** Decay rate of the fitted model at a frequency. */
255export function fittedRate(fit: SigmaFit, f: number): number {
256 return fit.sigma0 + (fit.sigma1 * 2 * Math.PI * f) / fit.kappa;
257}