1// Per-mode decay rates of the linear model from the mode shapes: the
2// radiation into the air and the thermoelastic loss. Both are exact for the
3// discrete model given the shape, so that the decay of each mode is a result
4// of the geometry and the material rather than a prescribed curve.
5//
6// Radiation. A baffled plate whose nodes move with the velocities v_k
7// radiates the power P = v^T R v with the radiation resistance matrix
8// R_jk = (rho_air omega^2 / (4 pi c)) dA_j dA_k sinc(k r_jk) (Elliott and
9// Johnson 1993), one face; both faces double it. For a mass-normalized mode
10// of amplitude A the velocity is omega A phi, the energy is omega^2 A^2 / 2,
11// and the amplitude decays at sigma = P / (2 E), which gives
12// sigma_rad = (rho_air omega^2 / (2 pi c)) sum_jk a_j a_k sinc(k r_jk),
13// with a_k = dA_k phi(k). On a regular grid the kernel depends only on the
14// offset j - k, so the double sum is a_j times a convolution, and by
15// Parseval it is a sum over the spectrum: two 2-D FFTs on a zero-padded
16// grid per mode instead of N^2 terms. The point-source treatment of a node's
17// area holds while k h is small; near coincidence the shape itself is barely
18// resolved, so nothing is lost there.
19//
20// Thermoelastic. Zener's loss factor applies to the bending energy and
21// depends on the local thickness, so a mode's rate is the sum over the nodes
22// of the local loss factor times the local bending energy, divided by omega:
23// for a mass-normalized mode the total strain energy is omega^2 / 2, and a
24// flat uniform plate recovers eta_te omega / 2.
26import { type Material, RHO_AIR, SPEED_OF_SOUND, bendingStiffness, type Plate, nodeWeight, insideCell, insideNode } from "../params.ts";
27import { thermoelasticLossFactor } from "../loss.ts";
28import { type ModalModel } from "./synth.ts";
30/** What the per-mode rates need of the plate, in transferable arrays. */
31export interface DampingPlate {
32 nx: number;
33 ny: number;
34 hx: number;
35 hy: number;
36 /** Node areas [m^2], zero outside the plate. */
37 area: Float64Array;
38 /** Thickness at the nodes [m], zero outside the plate. */
39 h: Float64Array;
40 mat: Material;
41 temperature: number;
42}
44export function dampingPlate(pl: Plate): DampingPlate {
45 const g = pl.grid;
46 const area = new Float64Array(g.n);
47 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) area[j * g.nx + i] = g.hx * g.hy * nodeWeight(g, i, j);
48 return { nx: g.nx, ny: g.ny, hx: g.hx, hy: g.hy, area, h: pl.h, mat: pl.mat, temperature: pl.temperature };
49}
51/** Fills the radiation and thermoelastic rates of every mode of the model (the single-thread path; the pool splits the modes among the workers). */
52export function computeDamping(model: ModalModel, dp: DampingPlate): void {
53 model.sigmaRad.set(radiationRates(model, dp, 0, model.count));
54 model.sigmaTe.set(thermoelasticRates(model, dp, 0, model.count));
55}
57// ---------- Radiation ----------
59const twiddleCache = new Map<number, { cos: Float64Array; sin: Float64Array; rev: Uint32Array }>();
61function twiddles(n: number) {
62 let t = twiddleCache.get(n);
63 if (t) return t;
64 const half = n >> 1, cos = new Float64Array(half), sin = new Float64Array(half), rev = new Uint32Array(n);
65 for (let k = 0; k < half; k++) { cos[k] = Math.cos((-2 * Math.PI * k) / n); sin[k] = Math.sin((-2 * Math.PI * k) / n); }
66 const bits = 31 - Math.clz32(n);
67 for (let i = 0; i < n; i++) { let r = 0; for (let b = 0; b < bits; b++) r |= ((i >> b) & 1) << (bits - 1 - b); rev[i] = r; }
68 t = { cos, sin, rev };
69 twiddleCache.set(n, t);
70 return t;
71}
73/** In-place radix-2 FFT of the strided sequence re[off + i * stride], im[...], i < n (n a power of two). */
74function fft1(re: Float64Array, im: Float64Array, off: number, stride: number, n: number): void {
75 const { cos, sin, rev } = twiddles(n);
76 for (let i = 0; i < n; i++) {
77 const r = rev[i];
78 if (r > i) {
79 const a = off + i * stride, b = off + r * stride;
80 let t = re[a]; re[a] = re[b]; re[b] = t;
81 t = im[a]; im[a] = im[b]; im[b] = t;
82 }
83 }
84 for (let size = 2; size <= n; size <<= 1) {
85 const half = size >> 1, step = n / size;
86 for (let start = 0; start < n; start += size) {
87 for (let k = 0; k < half; k++) {
88 const c = cos[k * step], s = sin[k * step];
89 const a = off + (start + k) * stride, b = off + (start + k + half) * stride;
90 const tr = re[b] * c - im[b] * s, ti = re[b] * s + im[b] * c;
91 re[b] = re[a] - tr; im[b] = im[a] - ti;
92 re[a] += tr; im[a] += ti;
93 }
94 }
95 }
96}
98/** In-place 2-D FFT of a Px by Py array stored row-major (index y * Px + x). */
99function fft2(re: Float64Array, im: Float64Array, Px: number, Py: number): void {
100 for (let y = 0; y < Py; y++) fft1(re, im, y * Px, 1, Px);
101 for (let x = 0; x < Px; x++) fft1(re, im, x, Px, Py);
102}
104function pow2AtLeast(n: number): number {
105 let p = 1;
106 while (p < n) p <<= 1;
107 return p;
108}
110/**
111 * The radiation decay rates of the modes [from, to) of the model by the FFT
112 * convolution described at the top of the file.
113 */
114export function radiationRates(model: ModalModel, dp: DampingPlate, from: number, to: number): Float64Array {
115 const { nx, ny, hx, hy } = dp;
116 const Px = pow2AtLeast(2 * nx - 1), Py = pow2AtLeast(2 * ny - 1);
117 const size = Px * Py;
118 const aRe = new Float64Array(size), aIm = new Float64Array(size);
119 const sRe = new Float64Array(size), sIm = new Float64Array(size);
120 const out = new Float64Array(to - from);
121 const c = SPEED_OF_SOUND;
122 for (let m = from; m < to; m++) {
123 const omega = 2 * Math.PI * model.freqs[m];
124 const kAc = omega / c;
125 // The weighted shape a_k = dA_k phi(k) on the padded grid.
126 aRe.fill(0); aIm.fill(0);
127 for (let r = 0; r < model.nRows; r++) {
128 const k = model.rowNode[r];
129 aRe[Math.floor(k / nx) * Px + (k % nx)] = dp.area[k] * model.wShapes[r * model.count + m];
130 }
131 fft2(aRe, aIm, Px, Py);
132 // The kernel sinc(k r) over the offsets, wrapped onto the padded grid; it is even, so its transform is real.
133 sIm.fill(0);
134 for (let dy = -(ny - 1); dy <= ny - 1; dy++) {
135 const yy = (dy + Py) % Py;
136 for (let dx = -(nx - 1); dx <= nx - 1; dx++) {
137 const x = kAc * Math.hypot(dx * hx, dy * hy);
138 sRe[yy * Px + ((dx + Px) % Px)] = x < 1e-9 ? 1 : Math.sin(x) / x;
139 }
140 }
141 // Offsets beyond the plate's extent are never reached; zero them.
142 for (let y = 0; y < Py; y++) {
143 const dy = y <= Py / 2 ? y : y - Py;
144 for (let x = 0; x < Px; x++) {
145 const dx = x <= Px / 2 ? x : x - Px;
146 if (Math.abs(dx) > nx - 1 || Math.abs(dy) > ny - 1) sRe[y * Px + x] = 0;
147 }
148 }
149 fft2(sRe, sIm, Px, Py);
150 let sum = 0;
151 for (let i = 0; i < size; i++) sum += sRe[i] * (aRe[i] * aRe[i] + aIm[i] * aIm[i]);
152 out[m - from] = ((RHO_AIR * omega * omega) / (2 * Math.PI * c)) * (sum / size);
153 }
154 return out;
155}
157/** The same rate by the direct double sum, for the tests and small plates. */
158export function radiationRateDirect(model: ModalModel, dp: DampingPlate, m: number): number {
159 const { nx, hx, hy } = dp;
160 const omega = 2 * Math.PI * model.freqs[m], c = SPEED_OF_SOUND, kAc = omega / c;
161 const a = new Float64Array(model.nRows), xs = new Float64Array(model.nRows), ys = new Float64Array(model.nRows);
162 for (let r = 0; r < model.nRows; r++) {
163 const k = model.rowNode[r];
164 a[r] = dp.area[k] * model.wShapes[r * model.count + m];
165 xs[r] = (k % nx) * hx; ys[r] = Math.floor(k / nx) * hy;
166 }
167 let sum = 0;
168 for (let i = 0; i < model.nRows; i++) {
169 for (let j = 0; j < model.nRows; j++) {
170 const x = kAc * Math.hypot(xs[i] - xs[j], ys[i] - ys[j]);
171 sum += a[i] * a[j] * (x < 1e-9 ? 1 : Math.sin(x) / x);
172 }
173 }
174 return ((RHO_AIR * omega * omega) / (2 * Math.PI * c)) * sum;
175}
177// ---------- Thermoelastic ----------
179/**
180 * The bending energy of a transverse shape w on the grid, split over the
181 * nodes (the mixed term of each cell is shared by its four corners), as in
182 * the discrete energy of the reference model.
183 */
184export function bendingEnergyPerNode(dp: DampingPlate, w: ArrayLike<number>, out?: Float64Array): Float64Array {
185 const { nx, ny, hx, hy, mat } = dp;
186 const g = { nx, ny, hx, hy, n: nx * ny, mask: maskOf(dp) };
187 const U = out ?? new Float64Array(nx * ny);
188 U.fill(0);
189 const id = (i: number, j: number) => j * nx + i;
190 const D = (h: number) => bendingStiffness(mat, h);
191 const nu = mat.nu;
192 for (let j = 0; j < ny; j++) {
193 for (let i = 0; i < nx; i++) {
194 const k = id(i, j);
195 const wgt = dp.area[k] / (hx * hy);
196 if (wgt === 0) continue;
197 let a = 0, b = 0;
198 if (insideNode(g, i - 1, j) && insideNode(g, i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[k] + w[id(i - 1, j)]) / (hx * hx);
199 if (insideNode(g, i, j - 1) && insideNode(g, i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[k] + w[id(i, j - 1)]) / (hy * hy);
200 U[k] += 0.5 * hx * hy * wgt * D(dp.h[k]) * (a * a + b * b + 2 * nu * a * b);
201 }
202 }
203 for (let cj = 0; cj < ny - 1; cj++) {
204 for (let ci = 0; ci < nx - 1; ci++) {
205 if (!insideCell(g, ci, cj)) continue;
206 const k00 = id(ci, cj), k10 = id(ci + 1, cj), k01 = id(ci, cj + 1), k11 = id(ci + 1, cj + 1);
207 const cc = (w[k11] - w[k10] - w[k01] + w[k00]) / (hx * hy);
208 const hc = 0.25 * (dp.h[k00] + dp.h[k10] + dp.h[k01] + dp.h[k11]);
209 const Uc = 0.5 * hx * hy * 2 * (1 - nu) * D(hc) * cc * cc;
210 U[k00] += Uc / 4; U[k10] += Uc / 4; U[k01] += Uc / 4; U[k11] += Uc / 4;
211 }
212 }
213 return U;
214}
216function maskOf(dp: DampingPlate): Uint8Array {
217 const m = new Uint8Array(dp.nx * dp.ny);
218 for (let k = 0; k < m.length; k++) m[k] = dp.area[k] > 0 ? 1 : 0;
219 return m;
220}
222/** The thermoelastic decay rates of the modes [from, to): sum_k eta_te(omega, h_k) U_k / omega with U_k the bending energy at node k. */
223export function thermoelasticRates(model: ModalModel, dp: DampingPlate, from: number, to: number): Float64Array {
224 const out = new Float64Array(to - from);
225 if (!dp.mat.thermal) return out;
226 const n = dp.nx * dp.ny;
227 const w = new Float64Array(n), U = new Float64Array(n);
228 for (let m = from; m < to; m++) {
229 const omega = 2 * Math.PI * model.freqs[m];
230 if (!(omega > 0)) continue;
231 w.fill(0);
232 for (let r = 0; r < model.nRows; r++) w[model.rowNode[r]] = model.wShapes[r * model.count + m];
233 bendingEnergyPerNode(dp, w, U);
234 let s = 0;
235 for (let k = 0; k < n; k++) if (U[k] !== 0) s += thermoelasticLossFactor(dp.mat, dp.h[k], omega, dp.temperature) * U[k];
236 out[m - from] = s / omega;
237 }
238 return out;
239}