Derive the losses from physical constants instead of a prescribed decay curve
design.loss.T60 is gone; a script that sets it is refused. The decay of each mode is now the sum of four mechanisms. Internal friction: a hysteretic loss factor of the metal (design.material.lossFactor), rate eta omega / 2 on the plate's share of the strain energy. Thermoelastic damping: Zener's relaxation from the thermal constants alpha, cp, k (design.material) and the temperature, in the biaxial plate form, applied per mode to the local bending energy with the local thickness. Radiation: computed per mode from the mode shape by the Rayleigh radiation resistance matrix of Elliott and Johnson, evaluated as an FFT convolution on the padded grid (two 2-D FFTs per mode, in the workers); it reproduces Wallace's simply supported efficiencies to 2 percent. The felt of the stand as before. The mode browser shows the decay time of each mode and of each mechanism alone.
The full model fits its two-parameter operator to an analytic estimate: the material rate, the area-averaged thermoelastic rate, and the Leppington average radiation efficiency of a plate of the design's area and marching-squares perimeter, times (f / f_c)^2 for free edges, a factor that the boundary terms of the far-field integral predict and that the modal values on the presets follow within a factor of two. A finding along the way, now in the theory note: a free-edged plate radiates 10 to 1000 times less than a simply supported one below coincidence, since its modes have no monopole or dipole moment, so on the presets the internal friction sets the decay below about 5 kHz and the radiation takes over only in the top octaves; the thermoelastic term dominates the lowest modes of the bronze cymbals.
Presets carry typical published loss factors and thermal constants for stainless steel, mild steel, bronze and aluminum. New tests cover Zener's peak, the perimeter, the efficiency limits, the fit, the piston closed form, the FFT against the direct double sum, Wallace's asymptotes, and the thermoelastic rate of a uniform plate. The tutorial's decay page explains the four mechanisms with the design's numbers.
24 changed files+940−283
README.mdmodified+16−10View file
@@ -1,6 +1,6 @@
11 # cymbalsim
22
3-A browser application for designing a plate or a cymbal and simulating a drumstick strike on it. The design is a short script in MATLAB syntax (material, bounding box, outline, thickness field, an optional stress-free profile, the stand, and the decay times), run in the page by [numbl](https://numbl.org), that the user edits in the page or that an AI writes for them. The sound is computed in one of two ways. The **linear model** computes the modes of the plate once per design and synthesizes the strike from them in about a second, on the CPU, in a pool of Web Workers. The **full model** solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly on the GPU through WebGPU; it takes tens of seconds per second of sound and captures the pitch glide and the wash of a hard hit. Both use the same grid, so the linear sound is the small-amplitude limit of the full one. The model and the numerical method are described in [THEORY.md](THEORY.md).
3+A browser application for designing a plate or a cymbal and simulating a drumstick strike on it. The design is a short script in MATLAB syntax (material with its loss constants, bounding box, outline, thickness field, an optional stress-free profile, and the stand), run in the page by [numbl](https://numbl.org), that the user edits in the page or that an AI writes for them. The sound is computed in one of two ways. The **linear model** computes the modes of the plate once per design and synthesizes the strike from them in about a second, on the CPU, in a pool of Web Workers. The **full model** solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly on the GPU through WebGPU; it takes tens of seconds per second of sound and captures the pitch glide and the wash of a hard hit. Both use the same grid, so the linear sound is the small-amplitude limit of the full one. The model and the numerical method are described in [THEORY.md](THEORY.md).
44
55 The focus is on free plates, sheets, discs and ellipses, flat by default, and on a realistic cymbal stand: a rod through the center hole with a felt pad of given radius below the cymbal and optionally another above it, which hold the cymbal level and let it bounce, rock, and (in the full model) lift off.
66
@@ -19,7 +19,7 @@ The sidebar holds the run controls at its top, then the design script, the strik
1919
2020 **Run mode.** The select at the top chooses the linear or the full model. In linear mode the modes of the current design are computed as soon as the script compiles (a progress bar shows the slices), and Synthesize produces the sound; after the first synthesis, a change of the strike or the listener parameters re-synthesizes automatically. In full mode Start begins the GPU run; Pause halts it, Stop ends it and keeps the audio and movie, and Continue extends a finished run after the duration has been raised.
2121
22-**Design script.** Pick a preset from the dropdown (a flat steel rectangle, a bowed rectangle, a flat bronze disc on a stand, a tapered 16 inch crash, a 20 inch ride, the same crash with a bow and a bell, a steel ellipse, a thunder sheet) and edit it freely; the script is compiled as you type (or with Ctrl+Enter) and the view shows the resulting stress-free shape. Errors are reported under the script, and a summary lists the bounding box, the material, the thickness range, the mass, the grid, the profile's peak-to-peak height and largest slope (with a warning above the shallow-shell range), the stand with its bounce frequency, and the decay times. Revert restores the selected preset.
22+**Design script.** Pick a preset from the dropdown (a flat steel rectangle, a bowed rectangle, a flat bronze disc on a stand, a tapered 16 inch crash, a 20 inch ride, the same crash with a bow and a bell, a steel ellipse, a thunder sheet) and edit it freely; the script is compiled as you type (or with Ctrl+Enter) and the view shows the resulting stress-free shape. Errors are reported under the script, and a summary lists the bounding box, the material, the thickness range, the mass, the grid, the profile's peak-to-peak height and largest slope (with a warning above the shallow-shell range), the stand with its bounce frequency, the loss constants with the coincidence frequency, and the estimated decay. Revert restores the selected preset.
2323
2424 The script fills in a struct named `design`. All lengths are meters. The outline, the thickness, and the profile are numbers or function handles `@(x, y, r, theta)` of the point relative to the center of the bounding box, its distance from the center, and its angle from the x axis; the handles are applied to whole grid matrices, so they use elementwise operators. A script is self-contained: nothing is in scope but numbl's own functions, so a script defines the shapes it needs as ordinary MATLAB functions at the bottom of the file, and the unit constants at the top.
2525
@@ -30,6 +30,10 @@ design.name = '16 inch flat crash';
3030 design.material.E = 110e9; % Pa
3131 design.material.nu = 0.34;
3232 design.material.rho = 8800; % kg/m^3
33+design.material.lossFactor = 2e-4; % internal friction of the metal (hysteretic loss factor)
34+design.material.alpha = 18e-6; % thermal expansion [1/K] } thermal constants, for the
35+design.material.cp = 370; % specific heat [J/(kg K)] } thermoelastic damping
36+design.material.k = 40; % conductivity [W/(m K)] }
3337 design.bounds.width = D; % bounding box [m]
3438 design.bounds.height = D;
3539 design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2); % omit for the whole box
@@ -42,20 +46,19 @@ design.stand.pad.E = 0.5e6; % Pa, felt in compression
4246 design.stand.pad.lossFactor = 0.3;
4347 design.stand.topPad = true; % false for a bottom pad only
4448 design.stand.preload = 1*mm; % compression of the top pad by the wing nut
45-design.loss.T60 = [100 5; 1000 3.5; 5000 1.2; 15000 0.4]; % decay times [Hz s], or a number, or @(f)
4649
4750 function v = ramp(a, b, x)
4851 v = min(max((x - a) ./ (b - a), 0), 1);
4952 end
5053 ```
5154
52-`design.stand = []` (or no stand field) gives a free plate; `design.stand.x` and `.y` move the rod. Pad fields that are omitted take the defaults shown. `design.loss.T60` is a decay time in seconds (a number), a table of `[frequency T60]` rows interpolated in log-log coordinates, or a function handle of the frequency in hertz; the default is 4 s at low frequencies falling to about 1 s at 5 kHz.
55+`design.stand = []` (or no stand field) gives a free plate; `design.stand.x` and `.y` move the rod. Pad fields that are omitted take the defaults shown. The decay of the sound is not prescribed but derived: from the loss factor of the metal's internal friction (`lossFactor`, 0 when omitted), the thermoelastic damping that follows from the three thermal constants (all three or none; without them there is no thermoelastic loss) and `design.temperature` (293 K by default), the radiation of each mode into the air, which the linear model computes from the mode's shape, and the felt of the stand. A script that still sets `design.loss` is refused with a message saying so.
5356
54-**Modes.** In linear mode the mode browser lists every mode up to the cutoff frequency with its decay time; click a row to draw the mode's shape on the plate (scaled to 1 mm), and filter the list by frequency. On a stand the first modes are the rock and the bounce of the cymbal on the felt. The modal frequencies are also drawn as ticks on the left edge of the spectrogram, in either run mode, so that the partials of a full run can be compared with the modes.
57+**Modes.** In linear mode the mode browser lists every mode up to the cutoff frequency with its decay time and, beside it, the decay time each mechanism alone would give: the air (radiation), the metal (internal friction), thermal (thermoelastic) and the felt of a stand; click a row to draw the mode's shape on the plate (scaled to 1 mm), and filter the list by frequency. On a stand the first modes are the rock and the bounce of the cymbal on the felt. The modal frequencies are also drawn as ticks on the left edge of the spectrogram, in either run mode, so that the partials of a full run can be compared with the modes.
5558
5659 **Strike and listener.** The strike point is given as fractions of the bounding box (the yellow ring in the view), with the stick speed, its effective mass, the Hertz contact stiffness and exponent, and the tip radius. The listener distance, its angle from the plate's normal, and its azimuth can be changed at any time; the microphone in the view moves with them, drawn nearer than the listener is so that both fit in the frame.
5760
58-**Numerics.** The grid spacing (8 mm by default), the cutoff frequency of the modes, the lift-off of the pads in the full model, the sample rate, and the frame rate. The Notes under the group report the grid, the frequencies at which grid dispersion lowers the partials by 3 and 10 percent, the count and time of the modes, the time step and step count of the full model, and the decay times of the fitted two-parameter loss next to the requested ones.
61+**Numerics.** The grid spacing (8 mm by default), the cutoff frequency of the modes, the lift-off of the pads in the full model, the sample rate, and the frame rate. The Notes under the group report the grid, the frequencies at which grid dispersion lowers the partials by 3 and 10 percent, the count and time of the modes, the time step and step count of the full model, and the decay times of the fitted two-parameter loss of the full model next to the estimate it was fitted to.
5962
6063 **Playback.** After a run, drag the time slider to explore the movie; Play starts the audio from the slider position and moves the movie and the spectrogram cursor with it; Design returns to the stress-free shape and, in full mode, Live to the running view. The audio can be played or downloaded as a 16-bit WAV at any time; Absolute level maps 1 Pa to full scale shifted by the gain slider (the peak level in dBFS is shown with a clipping warning), and the two normalized modes scale each run to its loudness or its peak. When the signal to be played is loud enough to hurt in headphones (a peak within 10 dB of full scale, or a short-term level of -20 dBFS RMS over some 50 ms window, which every loudness-normalized hit reaches), Play first shows a warning to turn the volume down, with a Play anyway button. Once accepted, the same signal replays without a prompt; a new recording, another level mode or a change of the gain brings the warning back.
6164
@@ -65,24 +68,27 @@ end
6568
6669 ## Tests
6770
68-`npm test` runs seven groups of tests with Node's built-in runner; `npm run test:cpu` runs the first six without a browser.
71+`npm test` runs eight groups of tests with Node's built-in runner; `npm run test:cpu` runs the first seven without a browser.
6972
70-- `tests/design.test.ts`: the design script (compilation, the checks with their messages, hoisted functions, the mask of a disc with a hole, the orientation of the sampled matrices, the removal of the mean plane, the stand's rates and preloads, the three forms of the decay times, the fit of the two-parameter loss), and that every preset evaluates on the default grid.
73+- `tests/design.test.ts`: the design script (compilation, the checks with their messages, hoisted functions, the mask of a disc with a hole, the orientation of the sampled matrices, the removal of the mean plane, the stand's rates and preloads, the loss constants of the material), and that every preset evaluates on the default grid.
7174 - `tests/reference.test.ts`: the float64 reference of the discretization (rigid modes of the stiffness, forces as the gradient of the energy, dissipativity of the losses and the pad dashpots, thickness scaling, the stand's stiffness against translation and tilt, the lift-off rule).
7275 - `tests/assemble.test.ts`: the assembled stiffness equals the linearized reference force on a flat plate and on a bowed, tapered, holed plate on a stand.
7376 - `tests/modal.test.ts`: the banded factorization and its inertia count, the tridiagonal eigensolver, all modes of a small square and a small disc against a dense solve (degenerate pairs included), the free square and disc frequencies against Gorman's and Leissa's values, and a plate on a stand.
7477 - `tests/synth.test.ts`: a single mode under the stick follows the closed-form damped sinusoid, the radiation gain of a piston equals the analytic Rayleigh sum, the attack and the tail join, the parts computed by the workers sum to the whole.
78+- `tests/loss.test.ts`: Zener's thermoelastic loss factor peaks where it should, the outline perimeter of a rectangle and a disc, the limits of the average radiation efficiency, the two-parameter fit, the radiation rate of a piston against the closed form and the FFT convolution against the direct double sum for every mode, and the thermoelastic rate of a uniform plate.
7579 - `tests/tutorial.test.ts`: every page of the tutorial renders with no design, with a plate, and with its modes, without undefined values and with balanced math delimiters; the actions run against an in-process host; the scaling pages' comparisons give the predicted ratios (2 for twice the thickness, 4 for half the size on a grid of the same cell count) to rounding.
7680 - `tests/gpu.test.ts`: serves the app with Vite and drives a headless Chromium through Playwright: the in-page self-test of the WGSL kernel against the float64 reference on five presets (`?selftest=1` shows the same numbers in the browser), the linear model's modes and synthesis, the tutorial panel with its typeset mathematics and an action that computes the modes, the loudness warning before playback, an error in a script, short full runs on a free plate and on a stand, the absence of a high-frequency noise floor in a curved design's full run, and a light hit simulated both ways, whose spectra must agree. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU.
7781
7882 ## Layout
7983
8084 - `src/params.ts`: simulation parameters, grid and mask, the plate as evaluated on the grid, material functions, the stand's rates, patches, listener geometry, grid dispersion.
81-- `src/design.ts`, `src/mruntime.ts`, `src/presets.ts`, `src/loss.ts`: running the script with numbl, the checks on the design struct, evaluation on the grid, the preset scripts, and the decay times.
85+- `src/design.ts`, `src/mruntime.ts`, `src/presets.ts`: running the script with numbl, the checks on the design struct, evaluation on the grid, and the preset scripts.
86+- `src/loss.ts`: the losses from physical constants (internal friction, Zener's thermoelastic damping, the average radiation efficiency, the outline perimeter) and the two-parameter fit of the full model.
8287 - `src/model/reference.ts`: float64 reference of the discretization (forces, losses, energy, the time-step estimate).
8388 - `src/model/assemble.ts`: the linearized stiffness in banded storage and the lumped mass.
8489 - `src/modal/banded.ts`, `lanczos.ts`, `modes.ts`: the banded factorization with its inertia count, shift-invert Lanczos on one slice, and the sliced computation of all modes through a job runner.
85-- `src/modal/synth.ts`, `worker.ts`, `pool.ts`: the modal synthesis (contact, radiation gains, attack table, tail, movie), the worker, and the pool that runs the eigenvalue jobs and splits the synthesis between the workers.
90+- `src/modal/damping.ts`: the per-mode radiation rate from the mode shape by FFT convolution of the radiation resistance, and the per-mode thermoelastic rate weighted by the local bending energy.
91+- `src/modal/synth.ts`, `worker.ts`, `pool.ts`: the modal synthesis (contact, radiation gains, attack table, tail, movie), the worker, and the pool that runs the eigenvalue jobs, the per-mode losses, and the synthesis split between the workers.
8692 - `src/sim/shaders/`, `src/sim/simulation.ts`, `src/sim/selftest.ts`: the WGSL kernels (step, listener sum, extraction), the WebGPU driver, and the in-browser comparison with the reference.
8793 - `src/render/`: 2D and 3D rendering of the plate with the markers of the strike, the stand, and the listener, with a 2D canvas fallback.
8894 - `src/audio/`: FFT, spectrogram, WAV encoder, playback, level mapping.
THEORY.mdmodified+11−5View file
@@ -41,7 +41,7 @@ which reduces to $D\nabla^4 w$ where $D$ is constant. Bending acts on the deviat
4141
4242 **Profile.** The stress-free shape is given by the script as a function of position. Only its slopes enter the strains, so its mean is immaterial, but its tilt is not: the von Karman strains are not invariant under rotation, and a tilted stress-free shape would add a spurious membrane coupling linear in the tilt. The mass-weighted least-squares plane is therefore subtracted from every profile. The treatment is that of a shallow shell, valid while the slopes are small compared with one; the interface reports the largest slope and warns above 0.1. The bell of a real crash rises 25 mm over a radius of 50 mm, a slope of about 0.5, so the curved preset uses a flatter bell, and this is the main point of tension between the model and a real cymbal (Section 7).
4343
44-**Losses.** The design gives the decay time $T_{60}(f)$, the time in which a partial at frequency $f$ decays by 60 dB, as a number, a table interpolated in log-log coordinates, or a function of frequency. The decay rate of a partial is then $\sigma(f) = \ln(1000)/T_{60}(f)$, so that its amplitude decays as $e^{-\sigma t}$. The linear model applies this rate to each mode at its own frequency (Section 4). The full model needs a loss operator rather than a rate per mode, and uses the standard two-parameter model of plate synthesis, a rate $\sigma_0 + \sigma_1 k^2$ for a bending wave of wavenumber $k$, with $\sigma_0$ and $\sigma_1$ fitted to the design's curve (Section 5). None of these values are measured; the presets are chosen to sound plausible, and the fit is reported next to the requested curve so that the difference is visible.
44+**Losses.** The decay of the sound is not prescribed. A partial of angular frequency $\omega$ decays as $e^{-\sigma t}$, and its decay time $T_{60} = \ln(1000)/\sigma$ is the time in which it falls by 60 dB; the rate $\sigma$ is the sum of four mechanisms, each with its own physics and its own constants in the script. *Internal friction* of the metal is a hysteretic loss: each cycle dissipates the fraction $2\pi\eta$ of the stored energy, so a mode decays at $\sigma = \eta\omega/2$, and the loss factor $\eta$ is a material constant on the same footing as $E$ and $\rho$, measured rather than derived, of order $10^{-4}$ for bronze, steel and aluminum. *Thermoelastic damping* comes from the heat that flows across the thickness between the compressed and the stretched face of a bent plate; Zener's analysis gives the loss factor $\eta_{te}(\omega, h) = \Delta\,\omega\tau/(1 + \omega^2\tau^2)$ with the relaxation strength $\Delta = [E\alpha^2 T/(\rho c_p)]\,(1+\nu)/(1-\nu)$, the biaxial form of the beam value, and the relaxation time $\tau = h^2\rho c_p/(\pi^2 k)$, the time heat takes to diffuse across the thickness, from the thermal expansion $\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$ (293 K unless the script says otherwise). The loss peaks at $\omega\tau = 1$, below 100 Hz for a millimeter of metal, and is small for the presets; it is what damps thin aluminum. *Radiation* into the air is the energy the plate loses as sound, which the linear model computes exactly for each mode from its shape (Section 4) and the full model takes from an analytic average (Section 5). It is smaller than one might expect. A plate with free edges radiates poorly below its coincidence frequency, the frequency at which the bending wavelength equals the wavelength of sound (about 17 kHz for 1 mm bronze): its modes are orthogonal to the rigid translation and rotations, so they have no net volume velocity and no dipole moment, and at a free edge the bending moment and the shear vanish, which removes the edge radiation that a simply supported plate has. The efficiency rises steeply toward coincidence, so on the presets the radiation overtakes the internal friction only in the top octaves, above about 5 kHz, and it is there that it makes the high partials die first; below that the decay of a free plate is set by the metal. The *felt* of a stand is the fourth mechanism (Section 3). The linear model adds the four rates for each mode, and the mode browser shows the decay time each would give alone. The full model needs a loss operator rather than a rate per mode and uses the standard two-parameter model of plate synthesis, a rate $\sigma_0 + \sigma_1 k^2$ for a bending wave of wavenumber $k$, fitted to the estimated curve (Section 5). The loss factors and thermal constants of the presets are typical published values, not measurements of a particular cymbal.
4545
4646 ## 2. Spatial discretization
4747
@@ -89,7 +89,13 @@ $$
8989 \ddot q_m + 2\sigma_m \dot q_m + \omega_m^2 q_m = \phi_m^T f(t),
9090 $$
9191
92-with $\sigma_m = \ln(1000)/T_{60}(f_m)$ from the design's decay times, plus the pad term of Section 3 on a stand. The listener signal and the movie are then sums over the modes.
92+with $\sigma_m$ the sum of the four loss rates of the mode. The internal friction gives $\eta U_m/\omega_m$ with $U_m$ the strain energy of the plate in the mode, $\omega_m^2/2$ for a mass-normalized mode less the energy in the pads; the pads give the term of Section 3; the thermoelastic loss gives $\sum_k \eta_{te}(\omega_m, h_k)\,U_{b,k}/\omega_m$, the local loss factor weighted by the bending energy $U_{b,k}$ of the mode at each node (the same discrete bending energy as Section 2, with the cell terms shared among their corners), which reduces to $\eta_{te}\omega_m/2$ on a flat uniform plate. The radiation rate is the power the mode radiates over twice its energy. A baffled plate whose nodes move with the velocities $v_k$ radiates the power $v^T R v$ from one face, with the radiation resistance matrix $R_{jk} = (\rho_{\mathrm{air}}\omega^2/4\pi c)\,\Delta A_j \Delta A_k\, \mathrm{sinc}(k_a r_{jk})$ of Elliott and Johnson, $k_a = \omega/c$; both faces radiate, the velocity of a mass-normalized mode of amplitude $A$ is $\omega A\phi$ and its energy $\omega^2 A^2/2$, so</p>
93+
94+$$
95+\sigma_{\mathrm{rad},m} = \frac{\rho_{\mathrm{air}}\,\omega_m^2}{2\pi c} \sum_{j,k} \Delta A_j\,\Delta A_k\,\phi_m(j)\,\phi_m(k)\,\frac{\sin(k_a r_{jk})}{k_a r_{jk}},
96+$$
97+
98+which for a small piston at low frequency reduces to the textbook resistance $\rho_{\mathrm{air}}\omega^2 S^2/(2\pi c)$ divided by the mass. On a regular grid the kernel depends only on the offset $j - k$, so the double sum is a convolution, and by Parseval it is a sum over the spectrum: two two-dimensional FFTs of the zero-padded grid per mode, about 5 million operations each for a 16 inch cymbal at 8 mm, split among the workers, rather than $N^2$ terms. The treatment of a node's area as a point source holds while $k_a h$ is small; near coincidence the mode shape itself is barely resolved by the grid, so nothing is lost there. The listener signal and the movie are then sums over the modes.
9399
94100 **Strike.** The stick tip is a free mass $M_s$ (30 g by default, the effective mass of a stick in a wrist stroke) arriving with speed $v_0$ at a point of the plate. The contact patch consists of the plate nodes within the tip radius (5 mm) plus the nearest one; its contact coordinate is the mass-weighted mean transverse displacement $w_c = \sum_k (m_k/M_p) w_k$ of these nodes, the penetration is $\delta = w_c - z$ with $z$ the tip position, and the force is the Hertz law $F = K\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ and $K = 10^7\ \mathrm{N/m^{1.5}}$ by default, a softer contact than the Hertz value for a hard wooden tip on metal (of order $10^8$ to $10^9$), which stands for the compliance of a real tip and gives a rounder attack; the value is a parameter of the strike. The force is spread over the patch in proportion to the node masses, which makes the pair conjugate: the modal force is $-F g_m$ and the contact coordinate is $\sum_m g_m q_m$ with the same gains $g_m = \sum_k (m_k/M_p)\, \phi_m(k)$. The stick obeys $M_s \ddot z = F$. The contact is stepped at 2 microseconds; over each step the force is held constant, and every oscillator is advanced by its exact solution under a constant force, so the step is accurate and unconditionally stable for the plate, and only the stick is integrated explicitly. Within the hit the stick may chatter, releasing and re-engaging as the plate runs ahead of it and comes back; once the contact has been open for 5 ms the stick is withdrawn, so that a plate rocking on its stand cannot strike it again tens of milliseconds later, as a drummer lifts the stick well before that. After that instant $t_1$ every mode decays freely, and its state is one complex amplitude, $q_m(t) = \mathrm{Re}\left[A_m e^{s_m (t - t_1)}\right]$ with $s_m = -\sigma_m + i\sqrt{\omega_m^2 - \sigma_m^2}$. A thin plate is a soft target: its driving-point impedance is $8\sqrt{D\rho h}$, about 60 N s/m for 1 mm bronze, so the contact lasts about a millisecond and the peak force is a few hundred newtons; the excitation is broadband, and the sound is a short bright attack followed by the ring. A hard hit in the linear model is simply a loud one; the model does not know that the plate would have stiffened.
95101
@@ -114,7 +120,7 @@ which is exact for the free decay and includes the directivity of every mode. Th
114120
115121 The full model keeps the nonlinear strains of Section 1 and steps the semi-discrete system with the symplectic Euler scheme, $v^{n+1} = v^n + \Delta t\, a(x^n, v^n)$ and $x^{n+1} = x^n + \Delta t\, v^{n+1}$, in single precision on the GPU. With the velocities interpreted at half steps this is the leapfrog scheme, second order and symplectic for the conservative part; the losses and the pad dashpots are evaluated explicitly at the old velocity. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system, estimated at start-up by power iteration on a small uniform grid with the same spacing and the largest thickness of the design; 90 percent of the limit is used. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)}$, about 3800 m/s for bronze, set this limit, and $c_p$ does not depend on the thickness: at 8 mm spacing $\Delta t \approx 1.9\ \mu$s for bronze, and 4 s of sound take about 2 million steps. The step is rounded down so that an integer number of steps makes one audio sample. The in-plane inertia terms are usually dropped in favor of an Airy stress function, since in-plane waves are much faster than bending waves; we keep them because they make the whole system explicit and local, which suits the GPU, at the cost of the time step just described.
116122
117-**Losses.** The two-parameter model gives a bending mode of angular frequency $\omega$ and wavenumber $k$ the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes the high frequencies decay faster in proportion to frequency. The two parameters are fitted by non-negative least squares to the design's rates $\ln(1000)/T_{60}(f)$ on log-spaced frequencies from 50 Hz to 10 kHz, with $\kappa$ at the mean thickness; a decay curve that is not of this form (a table with a kink, for instance) is matched only approximately, and the interface reports the fitted decay times at 100 Hz, 1 kHz and 5 kHz next to the requested ones. The $\sigma_1$ term is implemented as $2\sigma_1$ times the graph Laplacian of the transverse velocity and does not depend on the thickness; where the plate is thicker, a given frequency corresponds to a smaller $k^2$, so the frequency-dependent loss is weaker there. The in-plane motion carries an additional frequency-independent rate, taken as the design's rate at 5 kHz: the stretching modes lie at several kilohertz, couple linearly to the transverse motion once the plate is curved, and would otherwise ring for the whole low-frequency decay time.
123+**Losses.** The two-parameter model gives a bending mode of angular frequency $\omega$ and wavenumber $k$ the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes the high frequencies decay faster in proportion to frequency. Note that the internal friction, $\eta\omega/2$, is exactly of the second form with $\sigma_1 = \eta\kappa/2$. The two parameters are fitted by non-negative least squares to an estimate of the total rate on log-spaced frequencies from 50 Hz to 10 kHz, with $\kappa$ at the mean thickness. The estimate has no modes to work with, so it takes the internal friction as is, the thermoelastic loss factor averaged over the plate area with the local thickness, and the radiation from the average radiation efficiency of a baffled plate of the design's area and perimeter (the perimeter measured along the marching-squares contour of the cell mask, so that a disc is not counted as a staircase). The average of Leppington, Broadbent and Heron for a simply supported plate grows about linearly with frequency below coincidence, peaks there at a value set by the plate's size in coincidence wavelengths, and tends to one above; the modal computation of Section 4 reproduces Wallace's classical efficiencies of simply supported mode shapes to 2 percent, and for the free-edged modes of the presets it gives values below that average by the factor $(f/f_c)^2$, within a factor of two from 100 Hz to 8 kHz, which is what the boundary terms of the far-field integral predict for edges at which the moment and the shear vanish. The estimate applies that factor. The interface reports the fitted decay times at 100 Hz, 1 kHz and 5 kHz next to the estimate, and the modal computation, when it has run, shows how far the exact per-mode rates scatter about it. The $\sigma_1$ term is implemented as $2\sigma_1$ times the graph Laplacian of the transverse velocity, whose eigenvalue on this grid is exactly $\omega/\kappa$ for a grid mode, so the intended rate holds for the discrete modes too; it does not depend on the thickness, and where the plate is thicker a given frequency corresponds to a smaller $k^2$, so the frequency-dependent loss is weaker there. The in-plane motion neither radiates nor has a thermoelastic loss in this model and carries only the internal friction, as a frequency-independent rate $\eta\omega_p/2$ at the lowest in-plane resonance $\omega_p = \pi c_p/L$ of the plate's longer side, since the stretching modes lie at several kilohertz and couple linearly to the transverse motion once the plate is curved.
118124
119125 **Free plate.** A sheet hanging from a thread is treated as free, and the reaction to the stick is distributed over the plate as the inertial load of a rigid body, which is the inertia-relief equilibrium of a free body: the force of the previous step is applied as the rigid acceleration field $a_k = F\,(1/M + y_k y_c/I_x + x_k x_c/I_y)$ that balances the patch load and its moments about the center, so the plate gains no net momentum or angular momentum. The one-step lag has zero net impulse. A plate on a stand takes the reaction through the pads. This is the counterpart of dropping the rigid modes in the linear model.
120126
@@ -124,7 +130,7 @@ The full model keeps the nonlinear strains of Section 1 and steps the semi-discr
124130
125131 ## 6. Validation
126132
127-The tests check, without a GPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the full model are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered, free and on a stand; that the loss terms and the pad dashpots only remove energy; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate on a stand; that the banded factorization solves and that its inertia counts the eigenvalues; that every mode of a small square and of a small disc, degenerate pairs included, matches a dense eigensolution to $10^{-7}$ with no mode missed; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that the stand resists a translation with its total stiffness and a tilt with the pads' second moment, and that the lift-off rule unloads the right pad; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, that the attack and the tail of the synthesized signal join, and that the parts computed by the workers sum to the whole.
133+The tests check, without a GPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the full model are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered, free and on a stand; that the loss terms and the pad dashpots only remove energy; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate on a stand; that the banded factorization solves and that its inertia counts the eigenvalues; that every mode of a small square and of a small disc, degenerate pairs included, matches a dense eigensolution to $10^{-7}$ with no mode missed; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that the stand resists a translation with its total stiffness and a tilt with the pads' second moment, and that the lift-off rule unloads the right pad; that Zener's loss factor peaks at $\omega\tau = 1$ with the value $\Delta/2$, that the outline perimeter of a disc is within a few percent of $2\pi R$, that the average radiation efficiency has its textbook limits, that the radiation rate of a uniform piston at low frequency equals the closed form $\rho_{\mathrm{air}}\omega^2 S^2/(2\pi c M)$ and that the FFT convolution equals the direct double sum to $10^{-9}$ for every mode, and that the thermoelastic rate of a flat uniform plate is $\eta_{te}\omega/2$; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, that the attack and the tail of the synthesized signal join, and that the parts computed by the workers sum to the whole.
128134
129135 With a GPU, a Playwright test runs one step of the WGSL kernel from a smooth random state of about 2.5 thicknesses amplitude, with the losses and the stand with lift-off, on the flat, bowed, tapered, curved and elliptical presets, and compares the accelerations with the float64 reference: they agree to between $1.5\times10^{-7}$ and $6\times10^{-7}$ relative, which is single-precision rounding. A short strike runs to completion on a free plate and on a stand. A light hit on the curved crash is run to the end and the top octave of its tail is required to lie more than 50 dB below the peak, which catches a cancellation in the strains of the kind described in Section 5. And a light hit (0.05 m/s) on the flat disc is simulated both ways: over the interval from 50 to 250 ms the magnitude spectra of the two listener signals below 6 kHz have a correlation of 0.985 and their levels agree within 13 percent. The remaining difference comes from the fitted two-parameter loss against the per-mode decay times, the pad dashpot against the loss factor, and the modal truncation at 16 kHz.
130136
@@ -132,6 +138,6 @@ With a GPU, a Playwright test runs one step of the WGSL kernel from a smooth ran
132138
133139 The linear model is linear: it has no pitch glide, no wash, no lift-off, and a hard hit is only a loud one. It is the right tool for the tone, the modes and the decay of a design, and for a first listen, and the full model is the right tool for the sound of a hard hit; the two agree for a light one. The modal truncation at the cutoff frequency removes the highest partials, which the grid places poorly anyway; the linear interpolation of the acceleration table during the first two milliseconds slightly blurs the attack.
134140
135-The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a strike that moves a 400 mm cymbal by a few millimeters is well within its range. The Marguerre treatment of the profile is a shallow-shell approximation, adequate for a bow of a few millimeters on a sheet and for a gentle cymbal bow, but strained by the bell of a real cymbal, whose slope is several times the range in which the approximation holds; the curved preset uses a flatter bell than a real cymbal, and the interface reports the largest slope of a design. A deep bell also behaves as a stiff cup that shifts the bell modes up in ways the shallow shell does not capture. The thickness field is treated as a slowly varying coefficient in a thin-plate model, which is appropriate for a lathed taper but not for a step or a groove on the scale of the grid. The hammering of a real cymbal, which leaves residual stresses and a slightly non-axisymmetric shape, is absent unless the script adds it. The loss model is empirical and its parameters are not calibrated against a real cymbal, and the full model matches the design's decay curve only as well as the two-parameter fit allows. Grid dispersion lowers the high partials by the amounts given in Section 2, and the staircase edge places the low frequencies within a few percent and splits modes that would be degenerate on a true disc.
141+The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a strike that moves a 400 mm cymbal by a few millimeters is well within its range. The Marguerre treatment of the profile is a shallow-shell approximation, adequate for a bow of a few millimeters on a sheet and for a gentle cymbal bow, but strained by the bell of a real cymbal, whose slope is several times the range in which the approximation holds; the curved preset uses a flatter bell than a real cymbal, and the interface reports the largest slope of a design. A deep bell also behaves as a stiff cup that shifts the bell modes up in ways the shallow shell does not capture. The thickness field is treated as a slowly varying coefficient in a thin-plate model, which is appropriate for a lathed taper but not for a step or a groove on the scale of the grid. The hammering of a real cymbal, which leaves residual stresses and a slightly non-axisymmetric shape, is absent unless the script adds it. The losses are derived from physical constants, but two of those constants are themselves empirical: the loss factor of the metal and, to a lesser degree, its thermal constants are typical published values rather than measurements of a particular alloy after hammering, and the loss factor is the one knob a user has on the decay. The radiation loss is that of a plate in a rigid wall, consistent with the listener model but an overstatement at low frequency for a plate in free air, where the two faces partly cancel; the losses to the air's viscosity and to the support other than the felt are absent; and the full model carries the losses only through the two-parameter fit to an average, so its decay agrees with the modal one on average and not mode by mode. Grid dispersion lowers the high partials by the amounts given in Section 2, and the staircase edge places the low frequencies within a few percent and splits modes that would be degenerate on a true disc.
136142
137143 The stand is a linear-elastic felt with a loss factor; real felt stiffens with compression, and the value of its modulus is an estimate. The rod and its sleeve, which stop the cymbal from sliding sideways once the felt's shear gives way, are not modeled, and neither is the friction of the felt. Gravity enters only through the preload. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The plate is unbaffled in reality, so the sound from its two faces partly cancels at low frequencies, and the Rayleigh integral overestimates the low-frequency level; an accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator. Finally, single precision limits the dynamic range of the tail of the full model's sound to roughly 100 dB below the peak, and the linear model, in double precision, has no such limit but no more physics.
package.jsonmodified+1−1View file
@@ -9,7 +9,7 @@
99 "build": "tsc --noEmit && vite build",
1010 "preview": "vite preview",
1111 "test": "node --test tests/*.test.ts",
12- "test:cpu": "node --test tests/design.test.ts tests/reference.test.ts tests/assemble.test.ts tests/modal.test.ts tests/synth.test.ts tests/tutorial.test.ts",
12+ "test:cpu": "node --test tests/design.test.ts tests/reference.test.ts tests/assemble.test.ts tests/modal.test.ts tests/synth.test.ts tests/loss.test.ts tests/tutorial.test.ts",
1313 "test:gpu": "node --test tests/gpu.test.ts"
1414 },
1515 "devDependencies": {
src/design.tsmodified+29−45View file
@@ -1,6 +1,7 @@
11 // The design: a MATLAB-syntax script, run by numbl, that fills in a struct
22 // named `design` with the material, the bounding box, the outline, the
3-// thickness field, the stress-free profile, the stand, and the losses.
3+// thickness field, the stress-free profile, the stand, and the loss
4+// constants of the material (from which the decay is derived, loss.ts).
45 // The outline, the thickness, and the profile are function handles
56 // @(x, y, r, theta), vectorized over grid matrices (or plain numbers).
67 // compileDesign runs the script and checks the struct; evaluatePlate runs it
@@ -14,14 +15,14 @@
1415
1516 import { executeCode, type RuntimeValue, type RuntimeTensor, type RuntimeStruct } from "numbl";
1617 import { RUNTIME_FILES } from "./mruntime.ts";
17-import { LOSS_FREQS, constantLoss, defaultLoss, lossFromTable } from "./loss.ts";
1818 import {
1919 type Grid,
20- type Loss,
2120 type Material,
2221 type Pad,
2322 type Plate,
2423 type Stand,
24+ type Thermal,
25+ ROOM_TEMPERATURE,
2526 insideNode,
2627 makeGrid,
2728 nodeWeight,
@@ -43,8 +44,8 @@ export interface Design {
4344 /** The profile when it is a plain number (or absent); null for a function. */
4445 profileConst: number | null;
4546 stand: Stand | null;
46- /** The losses when they do not need the script (a number, a table, or the default); null for a function handle. */
47- lossConst: Loss | null;
47+ /** Temperature of the plate [K], for the thermoelastic damping. */
48+ temperature: number;
4849 source: string;
4950 }
5051
@@ -142,9 +143,27 @@ export function compileDesign(source: string): Design {
142143 const name = (nameV !== undefined ? text(nameV) : null)?.trim() || "Untitled design";
143144
144145 const m = field(d, "material");
145- if (!isStruct(m)) throw new DesignError("design.material must be a struct with fields E, nu, rho");
146- const material: Material = { E: num(field(m, "E"), "design.material.E", 1e6), nu: num(field(m, "nu"), "design.material.nu", 0), rho: num(field(m, "rho"), "design.material.rho", 1) };
146+ if (!isStruct(m)) throw new DesignError("design.material must be a struct with fields E, nu, rho, lossFactor, and the thermal constants alpha, cp, k");
147+ const has = (s: RuntimeStruct, k: string) => field(s, k) !== undefined && !isEmpty(field(s, k));
148+ const thermalKeys = ["alpha", "cp", "k"];
149+ const given = thermalKeys.filter((k) => has(m, k));
150+ let thermal: Thermal | null = null;
151+ if (given.length === 3) {
152+ thermal = { alpha: num(field(m, "alpha"), "design.material.alpha (thermal expansion, 1/K)", 0), cp: num(field(m, "cp"), "design.material.cp (specific heat, J/(kg K))", 1), k: num(field(m, "k"), "design.material.k (thermal conductivity, W/(m K))", 1e-3) };
153+ } else if (given.length > 0) {
154+ throw new DesignError("design.material needs all three thermal constants for the thermoelastic damping, alpha (1/K), cp (J/(kg K)) and k (W/(m K)), or none of them");
155+ }
156+ const material: Material = {
157+ E: num(field(m, "E"), "design.material.E", 1e6),
158+ nu: num(field(m, "nu"), "design.material.nu", 0),
159+ rho: num(field(m, "rho"), "design.material.rho", 1),
160+ eta: has(m, "lossFactor") ? num(field(m, "lossFactor"), "design.material.lossFactor", 0) : 0,
161+ thermal,
162+ };
147163 if (material.nu >= 0.5) throw new DesignError("design.material.nu must be below 0.5");
164+ if (material.eta >= 0.1) throw new DesignError("design.material.lossFactor must be small (metals have loss factors of 1e-4 to 1e-3)");
165+ if (field(d, "loss") !== undefined) throw new DesignError("design.loss is no longer used: the decay of each mode follows from design.material.lossFactor, the thermal constants alpha, cp and k, the radiation of the plate, and the felt of the stand");
166+ const temperature = has(d, "temperature") ? num(field(d, "temperature"), "design.temperature (K)", 1) : ROOM_TEMPERATURE;
148167
149168 const b = field(d, "bounds");
150169 if (!isStruct(b)) throw new DesignError("design.bounds must be a struct with fields width, height");
@@ -177,9 +196,8 @@ export function compileDesign(source: string): Design {
177196 }
178197
179198 const stand = readStand(field(d, "stand"));
180- const lossConst = readLoss(field(d, "loss"));
181199
182- return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, stand, lossConst, source };
200+ return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, stand, temperature, source };
183201 }
184202
185203 function readStand(st: RuntimeValue | undefined): Stand | null {
@@ -208,31 +226,6 @@ function readStand(st: RuntimeValue | undefined): Stand | null {
208226 };
209227 }
210228
211-/** The losses when they can be read without the script: null when T60 is a function handle. */
212-function readLoss(lo: RuntimeValue | undefined): Loss | null {
213- if (lo === undefined || isEmpty(lo)) return defaultLoss();
214- if (!isStruct(lo)) throw new DesignError("design.loss must be a struct with the field T60: a decay time in seconds, a table [f T60; ...], or a function handle @(f)");
215- const t = field(lo, "T60");
216- if (t === undefined || isEmpty(t)) return defaultLoss();
217- if (isFunction(t)) return null;
218- const s = scalar(t);
219- if (s !== null) {
220- if (!(s > 0) || !isFinite(s)) throw new DesignError("design.loss.T60 must be a positive decay time in seconds");
221- return constantLoss(s);
222- }
223- if (isTensor(t) && t.shape.length === 2 && t.shape[1] === 2 && t.shape[0] >= 1) {
224- const rows = t.shape[0];
225- const fs: number[] = [], ts: number[] = [];
226- for (let r = 0; r < rows; r++) {
227- const f = t.data[r], t60 = t.data[rows + r]; // column-major
228- if (!(f > 0) || !(t60 > 0) || !isFinite(f) || !isFinite(t60)) throw new DesignError("design.loss.T60 table rows must be [frequency T60] with positive values");
229- fs.push(f); ts.push(t60);
230- }
231- return lossFromTable(fs, ts);
232- }
233- throw new DesignError("design.loss.T60 must be a number, a two-column table [f T60; ...], or a function handle @(f)");
234-}
235-
236229 // ---------- Evaluation on a grid ----------
237230
238231 /** Subtracts the weighted least-squares plane a + b x + c y from f over the plate nodes. */
@@ -277,9 +270,6 @@ function trailer(d: Design, g: Grid): string {
277270 "cymbalsim_IN = cymbalsim_field(design.inside, 'inside', cymbalsim_XC, cymbalsim_YC, hypot(cymbalsim_XC, cymbalsim_YC), atan2(cymbalsim_YC, cymbalsim_XC));",
278271 ].join("\n")
279272 : "",
280- d.lossConst === null
281- ? [`cymbalsim_F = [${Array.from(LOSS_FREQS, fmt).join(" ")}];`, "cymbalsim_T60 = cymbalsim_t60(design.loss.T60, cymbalsim_F);"].join("\n")
282- : "",
283273 "",
284274 ].join("\n");
285275 }
@@ -296,7 +286,7 @@ function matrix(vars: Record<string, RuntimeValue>, name: string, rows: number,
296286 * Samples the design on a grid of about the requested spacing: the mask from
297287 * the outline at the cell centers, the thickness at the nodes (which must be
298288 * positive on the plate), the profile at the nodes with its mean plane
299- * removed, and the decay times when they are a function. Only the slopes of
289+ * removed. Only the slopes of
300290 * the profile enter the physics, so its mean is immaterial, but its tilt is
301291 * not: the von Karman strains are not invariant under rotation, and a tilted
302292 * stress-free shape would add a spurious membrane coupling.
@@ -313,12 +303,6 @@ export function evaluatePlate(d: Design, spacing: number): Plate {
313303 const IN = matrix(vars, "cymbalsim_IN", ny - 1, nx - 1, "inside");
314304 g = withCellMask(full, (ci, cj) => IN[ci * (ny - 1) + cj] !== 0);
315305 }
316- let loss = d.lossConst;
317- if (loss === null) {
318- const T = matrix(vars, "cymbalsim_T60", 1, LOSS_FREQS.length, "design.loss.T60");
319- for (let i = 0; i < T.length; i++) if (!(T[i] > 0) || !isFinite(T[i])) throw new DesignError(`design.loss.T60 must be positive; it is ${T[i]} at ${LOSS_FREQS[i].toFixed(0)} Hz`);
320- loss = { freqs: LOSS_FREQS, t60: Float64Array.from(T) };
321- }
322306 const h = new Float64Array(g.n), w0 = new Float64Array(g.n);
323307 let any = false;
324308 for (let j = 0; j < ny; j++) {
@@ -338,7 +322,7 @@ export function evaluatePlate(d: Design, spacing: number): Plate {
338322 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
339323 removePlane(w0, g);
340324 for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) if (!insideNode(g, i, j)) w0[j * nx + i] = 0;
341- return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, stand: d.stand, loss };
325+ return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, stand: d.stand, temperature: d.temperature };
342326 }
343327
344328 /** Whether the plate is flat (no profile), in which case bending and in-plane motion decouple in the linear model. */
src/loss.tsmodified+212−55View file
@@ -1,73 +1,219 @@
1-// Losses given as decay times. The design script gives T60, the time in
2-// which a partial decays by 60 dB, as a number, a table [f T60] or a
3-// function of frequency; the design layer samples it on LOSS_FREQS, and this
4-// module interpolates it, converts it to decay rates for the modes, and fits
5-// the two-parameter model of the full simulation to it.
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,
13+// which the full model's loss operator is fitted to. The felt of a stand is
14+// a fourth mechanism, handled with the stand (params.ts, synth.ts).
615
7-import { type Loss, type Material, bendingConstant } from "./params.ts";
8-
9-/** Frequencies at which T60 is sampled: 128 points from 10 Hz to 40 kHz, log-spaced. */
10-export const LOSS_FREQS: Float64Array = (() => {
11- const n = 128, f = new Float64Array(n);
12- for (let i = 0; i < n; i++) f[i] = 10 * Math.pow(4000, i / (n - 1));
13- return f;
14-})();
16+import { type Grid, type Material, type Plate, RHO_AIR, SPEED_OF_SOUND, bendingConstant, insideCell, massProperties, meanThickness, nodeWeight } from "./params.ts";
1517
1618 /** ln(1000): a rate sigma [1/s] and a decay time relate by sigma = LN1000 / T60. */
1719 export const LN1000 = Math.log(1000);
1820
19-/** The default losses: 4 s at low frequencies, falling to about 1 s at 5 kHz. */
20-export function defaultT60(f: number): number {
21- return 4 / (1 + f / 2000);
21+/** 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). */
22+export function t60Text(seconds: number): string {
23+ if (!isFinite(seconds) || seconds <= 0) return "∞";
24+ if (seconds < 10) return seconds.toFixed(2);
25+ if (seconds < 1000) return seconds.toFixed(0);
26+ return seconds.toExponential(0).replace("e+", "e");
27+}
28+
29+/** The decay time of a rate for display (see t60Text). */
30+export function rateText(sigma: number): string {
31+ return t60Text(sigma > 0 ? LN1000 / sigma : Infinity);
32+}
33+
34+// ---------- Internal friction ----------
35+
36+/** 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. */
37+export function materialRate(eta: number, omega: number): number {
38+ return (eta * omega) / 2;
39+}
40+
41+// ---------- Thermoelastic damping ----------
42+
43+/**
44+ * Zener's thermoelastic relaxation of a plate of thickness h: the relaxation
45+ * strength delta and the relaxation time tau. The strength carries the
46+ * factor (1 + nu) / (1 - nu) of the biaxial stress state of a bent plate
47+ * (Zener's beam value is E alpha^2 T / (rho cp)); the time is the thermal
48+ * diffusion time across the thickness, h^2 / (pi^2 chi) with the diffusivity
49+ * chi = k / (rho cp).
50+ */
51+export function thermoelasticRelaxation(mat: Material, h: number, T: number): { delta: number; tau: number } {
52+ const th = mat.thermal;
53+ if (!th) return { delta: 0, tau: 1 };
54+ const delta = ((mat.E * th.alpha * th.alpha * T) / (mat.rho * th.cp)) * ((1 + mat.nu) / (1 - mat.nu));
55+ const tau = (h * h * mat.rho * th.cp) / (Math.PI * Math.PI * th.k);
56+ return { delta, tau };
57+}
58+
59+/** 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. */
60+export function thermoelasticLossFactor(mat: Material, h: number, omega: number, T: number): number {
61+ const { delta, tau } = thermoelasticRelaxation(mat, h, T);
62+ const x = omega * tau;
63+ return (delta * x) / (1 + x * x);
64+}
65+
66+// ---------- Geometry of the outline ----------
67+
68+/**
69+ * The length of the outline of the plate, from the marching-squares contour
70+ * of the cell mask at level one half (cells outside the grid count as
71+ * outside). A staircase count would overstate the perimeter of a disc by
72+ * 4 / pi; the contour, whose segments run in eight directions, overstates
73+ * it by a few percent, which is well within the accuracy of the average
74+ * radiation efficiency it feeds.
75+ */
76+export function outlinePerimeter(g: Grid): number {
77+ const cx = g.nx - 1, cy = g.ny - 1;
78+ const inside = (ci: number, cj: number) => ci >= 0 && cj >= 0 && ci < cx && cj < cy && insideCell(g, ci, cj);
79+ const hx = g.hx, hy = g.hy;
80+ let length = 0;
81+ // Squares between the centers of cells (ci, cj), (ci+1, cj), (ci+1, cj+1), (ci, cj+1), one ring beyond the grid.
82+ for (let cj = -1; cj < cy; cj++) {
83+ for (let ci = -1; ci < cx; ci++) {
84+ 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;
85+ const code = a | (b << 1) | (c << 2) | (d << 3);
86+ if (code === 0 || code === 15) continue;
87+ // Midpoints of the square's edges: bottom (a-b), right (b-c), top (c-d), left (d-a), as offsets in cells.
88+ const bottom = [0.5, 0], right = [1, 0.5], top = [0.5, 1], left = [0, 0.5];
89+ const seg = (p: number[], q: number[]) => { length += Math.hypot((q[0] - p[0]) * hx, (q[1] - p[1]) * hy); };
90+ switch (code) {
91+ case 1: case 14: seg(left, bottom); break;
92+ case 2: case 13: seg(bottom, right); break;
93+ case 3: case 12: seg(left, right); break;
94+ case 4: case 11: seg(right, top); break;
95+ case 6: case 9: seg(bottom, top); break;
96+ case 7: case 8: seg(left, top); break;
97+ case 5: seg(left, bottom); seg(right, top); break; // saddle: two corners
98+ case 10: seg(bottom, right); seg(left, top); break;
99+ }
100+ }
101+ }
102+ return length;
103+}
104+
105+// ---------- Radiation ----------
106+
107+/** 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). */
108+export function coincidenceFrequency(mat: Material, h: number): number {
109+ return (SPEED_OF_SOUND * SPEED_OF_SOUND) / (2 * Math.PI * bendingConstant(mat, h));
22110 }
23111
24-export function lossFromFunction(fn: (f: number) => number): Loss {
25- const t60 = new Float64Array(LOSS_FREQS.length);
26- for (let i = 0; i < t60.length; i++) t60[i] = fn(LOSS_FREQS[i]);
27- return { freqs: LOSS_FREQS, t60 };
112+/**
113+ * The radiation efficiency of a baffled plate averaged over its modes near
114+ * the frequency f. The simply supported case is that of Leppington,
115+ * Broadbent and Heron (1982), as given by Fahy and Gardonio: below
116+ * coincidence the edges radiate and the efficiency grows about linearly with
117+ * frequency, at coincidence it peaks at a value set by the size of the plate
118+ * in coincidence wavelengths, and above it tends to one. A plate with free
119+ * edges radiates much less below coincidence: its modes are orthogonal to
120+ * the rigid translation and rotations, so they have no monopole or dipole
121+ * moment, and at a free edge the bending moment and the shear vanish, so
122+ * the boundary terms that make a simply supported edge radiate are smaller
123+ * by the factor (k / k_b)^2 in amplitude, k the acoustic and k_b the bending
124+ * wavenumber, that is by (f / f_c)^2 in power. The per-mode computation of
125+ * modal/damping.ts on the presets follows the simply supported curve times
126+ * this factor within a factor of two from 100 Hz to 8 kHz, and that is what
127+ * `freeEdges` applies. `a` and `b` are the shorter and the longer side.
128+ */
129+export function averageRadiationEfficiency(f: number, fc: number, area: number, perimeter: number, a: number, b: number, freeEdges = true): number {
130+ const c = SPEED_OF_SOUND;
131+ const lambdaC = c / fc;
132+ const cap = Math.max(1, (0.5 - 0.15 * (a / b)) * Math.sqrt(a / lambdaC));
133+ if (f >= fc) return Math.min(cap, 1 / Math.sqrt(Math.max(1e-12, 1 - fc / f)));
134+ const mu = Math.sqrt(fc / f);
135+ 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));
136+ return Math.min(cap, s) * (freeEdges ? (f / fc) ** 2 : 1);
28137 }
29138
30-export function defaultLoss(): Loss {
31- return lossFromFunction(defaultT60);
139+/** 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). */
140+export function radiationRateFromEfficiency(efficiency: number, mat: Material, h: number): number {
141+ return (RHO_AIR * SPEED_OF_SOUND * efficiency) / (mat.rho * h);
32142 }
33143
34-export function constantLoss(t60: number): Loss {
35- return lossFromFunction(() => t60);
144+// ---------- The estimate for a whole plate ----------
145+
146+/** Geometric quantities of a plate that the loss estimate needs. */
147+export interface PlateGeometry {
148+ area: number;
149+ perimeter: number;
150+ hMean: number;
151+ kappa: number;
152+ fc: number;
153+ /** Shorter and longer side of the bounding box [m]. */
154+ a: number;
155+ b: number;
156+ /** Node areas over the grid [m^2] (zero outside the plate), for the area-weighted thermoelastic loss. */
157+ nodeArea: Float64Array;
36158 }
37159
38-/** Log-log interpolation of a table of (frequency, T60) rows, clamped at the ends. */
39-export function lossFromTable(freqs: ArrayLike<number>, t60s: ArrayLike<number>): Loss {
40- const n = freqs.length;
41- const idx = Array.from({ length: n }, (_, i) => i).sort((a, b) => freqs[a] - freqs[b]);
42- const fs = idx.map((i) => freqs[i]), ts = idx.map((i) => t60s[i]);
43- return lossFromFunction((f) => interpLogLog(fs, ts, f));
160+export function plateGeometry(pl: Plate): PlateGeometry {
161+ const g = pl.grid;
162+ const { area } = massProperties(pl);
163+ const hMean = meanThickness(pl);
164+ const nodeArea = new Float64Array(g.n);
165+ 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);
166+ 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 };
44167 }
45168
46-function interpLogLog(fs: number[], ts: number[], f: number): number {
47- const n = fs.length;
48- if (n === 1 || f <= fs[0]) return ts[0];
49- if (f >= fs[n - 1]) return ts[n - 1];
50- let i = 0;
51- while (i < n - 2 && fs[i + 1] < f) i++;
52- const t = (Math.log(f) - Math.log(fs[i])) / (Math.log(fs[i + 1]) - Math.log(fs[i]));
53- return Math.exp(Math.log(ts[i]) * (1 - t) + Math.log(ts[i + 1]) * t);
169+/** The three contributions to the decay rate of a bending partial at frequency f [1/s], and their sum. */
170+export interface LossEstimate {
171+ material: number;
172+ thermoelastic: number;
173+ radiation: number;
174+ total: number;
54175 }
55176
56-/** T60 at a frequency, by log-log interpolation of the sampled curve, clamped at its ends. */
57-export function t60At(loss: Loss, f: number): number {
58- return interpLogLog(Array.from(loss.freqs), Array.from(loss.t60), f);
177+/**
178+ * The decay rate of a bending partial at frequency f estimated without the
179+ * modes: the material rate, the thermoelastic rate averaged over the plate
180+ * area with the local thickness, and the radiation rate from the average
181+ * efficiency at the mean thickness. The modal computation replaces the last
182+ * two by exact per-mode values; this estimate serves the full model's loss
183+ * operator and the summary of a design before its modes exist.
184+ */
185+export function estimatedRate(pl: Plate, geom: PlateGeometry, f: number): LossEstimate {
186+ const omega = 2 * Math.PI * f;
187+ const material = materialRate(pl.mat.eta, omega);
188+ let te = 0;
189+ if (pl.mat.thermal) {
190+ let num = 0, den = 0;
191+ for (let k = 0; k < geom.nodeArea.length; k++) {
192+ const a = geom.nodeArea[k];
193+ if (a === 0) continue;
194+ num += a * thermoelasticLossFactor(pl.mat, pl.h[k], omega, pl.temperature);
195+ den += a;
196+ }
197+ te = den > 0 ? ((num / den) * omega) / 2 : 0;
198+ }
199+ const radiation = radiationRateFromEfficiency(averageRadiationEfficiency(f, geom.fc, geom.area, geom.perimeter, geom.a, geom.b), pl.mat, geom.hMean);
200+ return { material, thermoelastic: te, radiation, total: material + te + radiation };
59201 }
60202
61-/** Decay rate sigma [1/s] of a partial at frequency f: its amplitude decays as exp(-sigma t). */
62-export function decayRate(loss: Loss, f: number): number {
63- return LN1000 / t60At(loss, f);
203+/** 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). */
204+export function inPlaneRate(pl: Plate): number {
205+ const cp = Math.sqrt(pl.mat.E / (pl.mat.rho * (1 - pl.mat.nu * pl.mat.nu)));
206+ const fIp = cp / (2 * Math.max(pl.Lx, pl.Ly));
207+ return materialRate(pl.mat.eta, 2 * Math.PI * fIp);
64208 }
65209
210+// ---------- The two-parameter operator of the full model ----------
211+
66212 /**
67213 * The two-parameter loss model of the full simulation. A bending mode of
68214 * angular frequency omega decays at the rate sigma0 + sigma1 k^2 with
69215 * k^2 = omega / kappa; the fit is a non-negative least-squares fit of that
70- * line to the design's decay rates on log-spaced frequencies from 50 Hz to
216+ * line to the estimated rates on log-spaced frequencies from 50 Hz to
71217 * 10 kHz, at the mean thickness of the plate.
72218 */
73219 export interface SigmaFit {
@@ -76,24 +222,35 @@ export interface SigmaFit {
76222 kappa: number; // [m^2/s] bending constant used for the fit
77223 }
78224
79-export function fitSigma(loss: Loss, mat: Material, hMean: number, fLo = 50, fHi = 10000, n = 40): SigmaFit {
80- const kappa = bendingConstant(mat, hMean);
225+export function fitSigma(pl: Plate, geom: PlateGeometry = plateGeometry(pl), fLo = 50, fHi = 10000, n = 40): SigmaFit {
226+ const kappa = geom.kappa;
81227 const x: number[] = [], y: number[] = [];
82228 for (let i = 0; i < n; i++) {
83229 const f = fLo * Math.pow(fHi / fLo, i / (n - 1));
84230 x.push((2 * Math.PI * f) / kappa); // k^2
85- y.push(decayRate(loss, f));
231+ y.push(estimatedRate(pl, geom, f).total);
86232 }
87- // Unconstrained least squares, then clamp: with two parameters the
88- // non-negative solution is the unconstrained one, or one parameter at
89- // zero with the other refit.
233+ return fitLine(x, y, kappa);
234+}
235+
236+/** 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. */
237+export function fitLine(x: number[], y: number[], kappa: number): SigmaFit {
238+ const n = x.length;
90239 let sx = 0, sy = 0, sxx = 0, sxy = 0;
91240 for (let i = 0; i < n; i++) { sx += x[i]; sy += y[i]; sxx += x[i] * x[i]; sxy += x[i] * y[i]; }
92241 const det = n * sxx - sx * sx;
93- let s1 = (n * sxy - sx * sy) / det, s0 = (sy - s1 * sx) / n;
94- if (s1 < 0) { s1 = 0; s0 = sy / n; }
95- if (s0 < 0) { s0 = 0; s1 = sxy / sxx; }
96- return { sigma0: s0, sigma1: s1, kappa };
242+ const candidates: [number, number][] = [[0, 0]];
243+ if (det > 0) { const s1 = (n * sxy - sx * sy) / det; candidates.push([(sy - s1 * sx) / n, s1]); }
244+ candidates.push([sy / n, 0]);
245+ if (sxx > 0) candidates.push([0, sxy / sxx]);
246+ let best: [number, number] = [0, 0], bestRes = Infinity;
247+ for (const [s0, s1] of candidates) {
248+ if (!(s0 >= 0 && s1 >= 0)) continue;
249+ let res = 0;
250+ for (let i = 0; i < n; i++) res += (y[i] - s0 - s1 * x[i]) ** 2;
251+ if (res < bestRes) { bestRes = res; best = [s0, s1]; }
252+ }
253+ return { sigma0: best[0], sigma1: best[1], kappa };
97254 }
98255
99256 /** Decay rate of the fitted model at a frequency. */
src/main.tsmodified+32−15View file
@@ -20,12 +20,13 @@ import {
2020 inPlaneSpeed,
2121 CFL,
2222 } from "./params.ts";
23-import { fitSigma, fittedRate, t60At, LN1000, decayRate } from "./loss.ts";
23+import { fitSigma, fittedRate, LN1000, plateGeometry, estimatedRate, inPlaneRate, rateText } from "./loss.ts";
24+import { dampingPlate, computeDamping } from "./modal/damping.ts";
2425 import { assemble, type LinearSystem } from "./model/assemble.ts";
2526 import { maxAngularFrequency, fullLoss, type FullLoss } from "./model/reference.ts";
2627 import { type ModeSet } from "./modal/modes.ts";
2728 import { WorkerPool } from "./modal/pool.ts";
28-import { modalModel, modeShapeOnGrid, type ModalModel, type ContactResult } from "./modal/synth.ts";
29+import { modalModel, modalRates, modeShapeOnGrid, type ModalModel, type ContactResult } from "./modal/synth.ts";
2930 import { Simulation } from "./sim/simulation.ts";
3031 import { runSelfTest } from "./sim/selftest.ts";
3132 import { PlateRenderer, type ViewMode, type ColorMode } from "./render/renderer.ts";
@@ -481,17 +482,18 @@ class App {
481482 const hit = key ? this.modesCache.get(key) : undefined;
482483 if (hit && this.plate) {
483484 this.installModes(hit.sys, hit.modes, hit.model, null);
484- if (this.pool) void this.pool.setModel(hit.model);
485+ if (this.pool) void this.pool.setModel(hit.model, null);
485486 } else this.modeBrowser.clear("Press Synthesize to compute the modes of this design.");
486487 }
487488
488489 /** Shows a mode set in the browser, the notes and the spectrogram ticks. */
489490 private installModes(sys: LinearSystem, modes: ModeSet, model: ModalModel, secs: number | null): void {
490- const pl = this.plate!;
491491 const p = this.form.get();
492492 const fMax = Math.min(p.modeMaxFreq, p.fs / 2);
493493 this.sys = sys; this.modes = modes; this.model = model;
494- const rows = Array.from(modes.freqs).map((f, m) => ({ freq: f, t60: LN1000 / Math.min(decayRate(pl.loss, f) + model.padSigma[m], 0.99 * 2 * Math.PI * f) }));
494+ const sigma = modalRates(model);
495+ const t60 = (s: number) => (s > 0 ? LN1000 / s : Infinity);
496+ const rows = Array.from(modes.freqs).map((f, m) => ({ freq: f, t60: t60(sigma[m]), radiation: t60(model.sigmaRad[m]), material: t60(model.sigmaMat[m]), thermoelastic: t60(model.sigmaTe[m]), felt: t60(model.padSigma[m]) }));
495497 const kind = sys.fields === 1 ? "bending modes" : "modes (bending coupled to stretching)";
496498 const when = secs === null ? "cached" : `in ${secs.toFixed(1)} s`;
497499 this.modeBrowser.setModes(rows, `${modes.count} ${kind} up to ${(fMax / 1e3).toFixed(1)} kHz, ${when} (${sys.n} unknowns${modes.missing ? `, ${modes.missing} missing` : ""}). Click a mode to see its shape.`);
@@ -521,7 +523,12 @@ class App {
521523 const pad = pl.stand.pad;
522524 lines.push(`Stand: ${pl.stand.topPad ? "two pads" : "bottom pad only"}, ${(pad.innerRadius * 1e3).toFixed(0)} to ${(pad.outerRadius * 1e3).toFixed(0)} mm, ${(pad.thickness * 1e3).toFixed(0)} mm felt at ${(pad.E / 1e6).toFixed(2)} MPa; ${(sr.totalKw / 1e3).toFixed(0)} kN/m, bounce ${sr.bounceHz.toFixed(0)} Hz, ${sr.nodes.length} pad nodes, preload ${(sr.preloadTop * 1e3).toFixed(2)} / ${(sr.preloadBottom * 1e3).toFixed(2)} mm`);
523525 } else lines.push("Free plate (no stand)");
524- lines.push(`Decay T60: ${t60At(pl.loss, 100).toFixed(2)} s at 100 Hz, ${t60At(pl.loss, 1000).toFixed(2)} s at 1 kHz, ${t60At(pl.loss, 5000).toFixed(2)} s at 5 kHz`);
526+ const geom = plateGeometry(pl);
527+ const est = (f: number) => estimatedRate(pl, geom, f);
528+ const t60 = rateText;
529+ const th = pl.mat.thermal;
530+ lines.push(`Losses: loss factor ${pl.mat.eta.toExponential(1)}${th ? `; thermal constants alpha ${(th.alpha * 1e6).toFixed(1)}e-6 /K, cp ${th.cp.toFixed(0)} J/(kg K), k ${th.k.toFixed(0)} W/(m K)` : "; no thermal constants, so no thermoelastic damping"}; coincidence at ${(geom.fc / 1e3).toFixed(1)} kHz, perimeter ${(geom.perimeter * 1e3).toFixed(0)} mm`);
531+ lines.push(`Estimated decay T60: ${t60(est(100).total)} / ${t60(est(1000).total)} / ${t60(est(5000).total)} s at 100 Hz / 1 kHz / 5 kHz, of which radiation alone ${t60(est(100).radiation)} / ${t60(est(1000).radiation)} / ${t60(est(5000).radiation)} s; the modes give the exact values.`);
525532 void p;
526533 $("designInfo").innerHTML = lines.join("\n");
527534 }
@@ -537,12 +544,14 @@ class App {
537544 const dtMax = (CFL * 2) / wmax;
538545 const sps = Math.ceil(1 / (p.fs * dtMax));
539546 const dt = 1 / (p.fs * sps);
540- const fit = fitSigma(pl.loss, pl.mat, hm);
547+ const geom = plateGeometry(pl);
548+ const fit = fitSigma(pl, geom);
541549 const t60f = (f: number) => (LN1000 / fittedRate(fit, f)).toFixed(2);
550+ const t60e = (f: number) => (LN1000 / estimatedRate(pl, geom, f).total).toFixed(2);
542551 const lines = [
543552 `Grid ${g.nx} x ${g.ny}, ${countInside(g)} plate nodes. Grid dispersion lowers the bending partials by 3% near ${(f3 / 1e3).toFixed(1)} kHz and by 10% near ${(f10 / 1e3).toFixed(1)} kHz (at the mean thickness).`,
544553 `Linear model: modes up to ${(p.modeMaxFreq / 1e3).toFixed(1)} kHz${this.modes ? `, ${this.modes.count} found in ${this.modes.seconds.toFixed(1)} s` : ""}${this.modes && this.modes.missing ? ` <span class="warn">(${this.modes.missing} missing)</span>` : ""}.`,
545- `Full model: time step ${(dt * 1e6).toFixed(2)} us (${sps} per sample, in-plane wave speed ${inPlaneSpeed(pl.mat).toFixed(0)} m/s), ${((p.duration * p.fs * sps) / 1e6).toFixed(1)} million steps for ${p.duration} s. Losses fitted as sigma0 = ${fit.sigma0.toFixed(2)} 1/s, sigma1 = ${fit.sigma1.toExponential(2)} m^2/s: T60 ${t60f(100)} / ${t60f(1000)} / ${t60f(5000)} s at 100 Hz / 1 kHz / 5 kHz (requested ${t60At(pl.loss, 100).toFixed(2)} / ${t60At(pl.loss, 1000).toFixed(2)} / ${t60At(pl.loss, 5000).toFixed(2)}).`,
554+ `Full model: time step ${(dt * 1e6).toFixed(2)} us (${sps} per sample, in-plane wave speed ${inPlaneSpeed(pl.mat).toFixed(0)} m/s), ${((p.duration * p.fs * sps) / 1e6).toFixed(1)} million steps for ${p.duration} s. Losses fitted to the estimate as sigma0 = ${fit.sigma0.toFixed(2)} 1/s, sigma1 = ${fit.sigma1.toExponential(2)} m^2/s: T60 ${t60f(100)} / ${t60f(1000)} / ${t60f(5000)} s at 100 Hz / 1 kHz / 5 kHz (estimate ${t60e(100)} / ${t60e(1000)} / ${t60e(5000)}); in-plane rate ${inPlaneRate(pl).toFixed(2)} 1/s.`,
546555 ];
547556 this.form.setNotes("Grid and numerics", lines.join("\n"));
548557 const patch = strikePatch(pl, p);
@@ -574,7 +583,7 @@ class App {
574583 if (!pl || !key) return Promise.resolve();
575584 const hit = this.modesCache.get(key);
576585 if (hit) {
577- if (this.model !== hit.model) { this.installModes(hit.sys, hit.modes, hit.model, null); if (this.pool) return this.pool.setModel(hit.model); }
586+ if (this.model !== hit.model) { this.installModes(hit.sys, hit.modes, hit.model, null); if (this.pool) return this.pool.setModel(hit.model, null); }
578587 return Promise.resolve();
579588 }
580589 const p = this.form.get();
@@ -604,10 +613,19 @@ class App {
604613 throw e;
605614 }
606615 const model = modalModel(pl, sys, modes);
616+ // The radiation and the thermoelastic loss of every mode, from its shape.
617+ this.setProgress(1, "Computing the radiation of the modes...");
618+ const dp = dampingPlate(pl);
619+ try {
620+ if (this.pool) { await this.pool.setModel(model, dp); await this.pool.rates(model); }
621+ else computeDamping(model, dp);
622+ } catch (e) {
623+ this.setProgress(null, `Loss computation failed: ${e instanceof Error ? e.message : String(e)}`);
624+ throw e;
625+ }
607626 this.modesCache.set(key, { sys, modes, model });
608627 while (this.modesCache.size > 4) this.modesCache.delete(this.modesCache.keys().next().value!);
609628 if (gen !== this.designGeneration) { this.setProgress(null); return; }
610- if (this.pool) await this.pool.setModel(model);
611629 this.setProgress(null, "");
612630 this.installModes(sys, modes, model, (performance.now() - t0) / 1000);
613631 });
@@ -661,7 +679,7 @@ class App {
661679 try {
662680 let audio: Float32Array, frames: Float32Array, nFrames: number, contact: ContactResult;
663681 if (this.pool) {
664- const r = await this.pool.synthesize(model, pl, pl.loss, p, (a, c) => {
682+ const r = await this.pool.synthesize(model, pl, p, (a, c) => {
665683 if (gen !== this.designGeneration) return;
666684 this.installRecording(a, [], p, "linear");
667685 this.lastContact = c;
@@ -670,7 +688,7 @@ class App {
670688 audio = r.audio; frames = r.frames; nFrames = r.nFrames; contact = r.contact;
671689 } else {
672690 const { synthesize } = await import("./modal/synth.ts");
673- const r = synthesize(model, pl, pl.loss, p, true);
691+ const r = synthesize(model, pl, p, true);
674692 audio = r.audio; frames = r.frames; nFrames = r.nFrames; contact = r.contact;
675693 }
676694 if (gen !== this.designGeneration) return;
@@ -699,9 +717,8 @@ class App {
699717 if (!pl || !this.device || this.running) return;
700718 const p = this.form.get();
701719 this.sim?.destroy();
702- const hm = meanThickness(pl);
703- const fit = fitSigma(pl.loss, pl.mat, hm);
704- const loss: FullLoss = fullLoss(fit, standRates(pl), decayRate(pl.loss, 5000));
720+ const fit = fitSigma(pl);
721+ const loss: FullLoss = fullLoss(fit, standRates(pl), inPlaneRate(pl));
705722 this.sim = new Simulation(this.device, pl, p, loss);
706723 this.rec = { audio: new Float32Array(0), samples: 0, fs: p.fs, frames: [new Float32Array(pl.grid.n)], fps: p.fps, source: "full" };
707724 this.resetLoudness();
src/modal/damping.tsadded+239−0View file
@@ -0,0 +1,239 @@
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.
25+
26+import { type Material, RHO_AIR, SPEED_OF_SOUND, bendingStiffness, type Plate, nodeWeight, insideCell, insideNode } from "../params.ts";
27+import { thermoelasticLossFactor } from "../loss.ts";
28+import { type ModalModel } from "./synth.ts";
29+
30+/** What the per-mode rates need of the plate, in transferable arrays. */
31+export 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+}
43+
44+export 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+}
50+
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). */
52+export 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+}
56+
57+// ---------- Radiation ----------
58+
59+const twiddleCache = new Map<number, { cos: Float64Array; sin: Float64Array; rev: Uint32Array }>();
60+
61+function 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+}
72+
73+/** In-place radix-2 FFT of the strided sequence re[off + i * stride], im[...], i < n (n a power of two). */
74+function 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+}
97+
98+/** In-place 2-D FFT of a Px by Py array stored row-major (index y * Px + x). */
99+function 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+}
103+
104+function pow2AtLeast(n: number): number {
105+ let p = 1;
106+ while (p < n) p <<= 1;
107+ return p;
108+}
109+
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+ */
114+export 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+}
156+
157+/** The same rate by the direct double sum, for the tests and small plates. */
158+export 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+}
176+
177+// ---------- Thermoelastic ----------
178+
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+ */
184+export 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+}
215+
216+function 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+}
221+
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. */
223+export 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+}
src/modal/pool.tsmodified+26−6View file
@@ -5,7 +5,8 @@
55 // modes, and the main thread sums them.
66
77 import { type LinearSystem } from "../model/assemble.ts";
8-import { type Loss, type Plate, type SimParams } from "../params.ts";
8+import { type Plate, type SimParams } from "../params.ts";
9+import { type DampingPlate } from "./damping.ts";
910 import { type Job, type JobResult, type JobRunner, type ModeSet, computeModes } from "./modes.ts";
1011 import {
1112 type ContactResult,
@@ -116,8 +117,8 @@ export class WorkerPool {
116117 return computeModes(sys, fMax, this.runner(), { onProgress });
117118 }
118119
119- /** Splits the modes of the model between the workers. */
120- async setModel(model: ModalModel): Promise<void> {
120+ /** Splits the modes of the model between the workers, with the plate when the per-mode losses are still to be computed. */
121+ async setModel(model: ModalModel, plate: DampingPlate | null): Promise<void> {
121122 const n = this.workers.length;
122123 this.modelSplit = [];
123124 const per = Math.ceil(model.count / n);
@@ -126,18 +127,37 @@ export class WorkerPool {
126127 const from = Math.min(model.count, i * per), to = Math.min(model.count, (i + 1) * per);
127128 this.modelSplit.push({ from, to });
128129 const slice = sliceModel(model, from, to);
129- return w.send({ type: "model", model: slice });
130+ return w.send({ type: "model", model: slice, plate });
130131 }),
131132 );
132133 }
133134
135+ /** The radiation and thermoelastic rates of the model's modes, computed by the workers on their slices and written into the model. */
136+ async rates(model: ModalModel): Promise<void> {
137+ const parts = await Promise.all(
138+ this.workers.map((w, i) => {
139+ const { from, to } = this.modelSplit[i];
140+ if (to <= from) return Promise.resolve(null);
141+ return w.send({ type: "rates", id: this.nextId++ });
142+ }),
143+ );
144+ parts.forEach((r, i) => {
145+ if (!r || r.type !== "rates") return;
146+ const { from } = this.modelSplit[i];
147+ model.sigmaRad.set(r.rad, from);
148+ model.sigmaTe.set(r.te, from);
149+ });
150+ // The workers' slices were sent before the rates existed; refresh them so the synthesis sees the same model.
151+ await this.setModel(model, null);
152+ }
153+
134154 /**
135155 * Synthesizes a strike with the model last set. The audio comes first
136156 * (through `onAudio`), then the movie frames; the returned promise
137157 * resolves when both are done.
138158 */
139- async synthesize(model: ModalModel, pl: Plate, loss: Loss, p: SimParams, onAudio?: (audio: Float32Array, contact: ContactResult) => void, withMovie = true): Promise<Synthesis> {
140- const sigma = modalRates(model, loss);
159+ async synthesize(model: ModalModel, pl: Plate, p: SimParams, onAudio?: (audio: Float32Array, contact: ContactResult) => void, withMovie = true): Promise<Synthesis> {
160+ const sigma = modalRates(model);
141161 const g = strikeGains(model, pl, p);
142162 const contact = simulateContact(model, sigma, g, p, p.fps);
143163 const { weight, delay, maxDelay } = listenerWeights(pl, p);
src/modal/synth.tsmodified+31−19View file
@@ -1,6 +1,7 @@
11 // Modal synthesis of a strike. Each mode m of the linearized plate is a
2-// damped oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma from the
3-// design's decay times (and the pads' loss factor on a stand). The stick is
2+// damped oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma the sum of
3+// the metal's internal friction, the thermoelastic loss, the radiation of
4+// the mode (damping.ts) and the pads' loss factor on a stand. The stick is
45 // the Hertz contact of the full simulation, stepped against the modal
56 // system with an exact update of every oscillator under a piecewise
67 // constant force. Once the stick has left, every mode is a damped sinusoid
@@ -14,8 +15,7 @@
1415
1516 import { type LinearSystem } from "../model/assemble.ts";
1617 import { type ModeSet } from "./modes.ts";
17-import { type Loss, type Plate, type SimParams, STICK_RELEASE, listenerTerms, standRates, strikePatch, nodeMass } from "../params.ts";
18-import { decayRate } from "../loss.ts";
18+import { type Plate, type SimParams, STICK_RELEASE, listenerTerms, standRates, strikePatch, nodeMass } from "../params.ts";
1919
2020 /** What the synthesis needs of a plate and its modes; plain arrays so that it can be posted to a worker. */
2121 export interface ModalModel {
@@ -30,6 +30,12 @@ export interface ModalModel {
3030 wRow: Int32Array;
3131 /** Grid node of each row. */
3232 rowNode: Int32Array;
33+ /** Decay rate of each mode from the metal's internal friction [1/s]. */
34+ sigmaMat: Float64Array;
35+ /** Decay rate of each mode from thermoelastic damping [1/s]; filled by damping.ts, zero until then. */
36+ sigmaTe: Float64Array;
37+ /** Decay rate of each mode from its radiation into the air [1/s]; filled by damping.ts, zero until then. */
38+ sigmaRad: Float64Array;
3339 /** Decay rate added by the pads' loss factor to each mode [1/s]; zero on a free plate. */
3440 padSigma: Float64Array;
3541 /** Node masses on the grid [kg]. */
@@ -49,13 +55,17 @@ export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalM
4955 }
5056 const masses = new Float64Array(nGrid);
5157 for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) masses[j * pl.grid.nx + i] = nodeMass(pl, i, j);
52- // Hysteretic damping of the pads: a loss factor eta gives each mode the
53- // rate sigma = eta (phi' K_pad phi) / (2 omega).
54- const padSigma = new Float64Array(M);
58+ // Hysteretic damping: a loss factor eta on a store of energy U gives a mode
59+ // the rate sigma = eta U / omega. The pads hold the energy phi' K_pad phi / 2
60+ // and the plate the rest of the mode's strain energy, omega^2 / 2 for a
61+ // mass-normalized mode; each has its own loss factor. The thermoelastic and
62+ // the radiation rates need the shapes on the grid and are filled in by
63+ // damping.ts.
64+ const padSigma = new Float64Array(M), sigmaMat = new Float64Array(M);
5565 const sr = standRates(pl);
56- if (sr) {
57- for (let m = 0; m < M; m++) {
58- let e = 0;
66+ for (let m = 0; m < M; m++) {
67+ let e = 0;
68+ if (sr) {
5969 for (const k of sr.nodes) {
6070 const d = sys.nodeDof[k];
6171 if (d < 0) continue;
@@ -66,11 +76,13 @@ export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalM
6676 e += sr.ku[k] * (u * u + v * v);
6777 }
6878 }
69- const omega = 2 * Math.PI * modes.freqs[m];
70- padSigma[m] = omega > 0 ? (sr.eta * e) / (2 * omega) : 0;
7179 }
80+ const omega = 2 * Math.PI * modes.freqs[m];
81+ if (!(omega > 0)) continue;
82+ padSigma[m] = sr ? (sr.eta * e) / (2 * omega) : 0;
83+ sigmaMat[m] = Math.max(0, pl.mat.eta * (omega / 2 - e / (2 * omega)));
7284 }
73- return { nGrid, nRows, count: M, freqs: modes.freqs, wShapes, wRow, rowNode: Int32Array.from(rows), padSigma, masses };
85+ return { nGrid, nRows, count: M, freqs: modes.freqs, wShapes, wRow, rowNode: Int32Array.from(rows), sigmaMat, sigmaTe: new Float64Array(M), sigmaRad: new Float64Array(M), padSigma, masses };
7486 }
7587
7688 /** The model restricted to the modes [from, to). */
@@ -78,7 +90,7 @@ export function sliceModel(model: ModalModel, from: number, to: number): ModalMo
7890 const M = to - from;
7991 const wShapes = new Float32Array(model.nRows * M);
8092 for (let r = 0; r < model.nRows; r++) wShapes.set(model.wShapes.subarray(r * model.count + from, r * model.count + to), r * M);
81- return { ...model, count: M, freqs: model.freqs.slice(from, to), wShapes, padSigma: model.padSigma.slice(from, to) };
93+ return { ...model, count: M, freqs: model.freqs.slice(from, to), wShapes, sigmaMat: model.sigmaMat.slice(from, to), sigmaTe: model.sigmaTe.slice(from, to), sigmaRad: model.sigmaRad.slice(from, to), padSigma: model.padSigma.slice(from, to) };
8294 }
8395
8496 /** The transverse shape of mode m on the grid (zero outside the plate). */
@@ -88,12 +100,12 @@ export function modeShapeOnGrid(model: ModalModel, m: number): Float32Array {
88100 return w;
89101 }
90102
91-/** Per-mode decay rates from the design's decay times plus the pads, kept below the mode's frequency (underdamped). */
92-export function modalRates(model: ModalModel, loss: Loss): Float64Array {
103+/** Per-mode decay rates: internal friction, thermoelastic, radiation and the pads, kept below the mode's frequency (underdamped). */
104+export function modalRates(model: ModalModel): Float64Array {
93105 const sigma = new Float64Array(model.count);
94106 for (let m = 0; m < model.count; m++) {
95107 const w = 2 * Math.PI * model.freqs[m];
96- sigma[m] = Math.min(decayRate(loss, model.freqs[m]) + model.padSigma[m], 0.99 * w);
108+ sigma[m] = Math.min(model.sigmaMat[m] + model.sigmaTe[m] + model.sigmaRad[m] + model.padSigma[m], 0.99 * w);
97109 }
98110 return sigma;
99111 }
@@ -406,8 +418,8 @@ export interface Synthesis {
406418 }
407419
408420 /** Runs the whole synthesis on the calling thread (the tests, and browsers without workers). */
409-export function synthesize(model: ModalModel, pl: Plate, loss: Loss, p: SimParams, withMovie = true): Synthesis {
410- const sigma = modalRates(model, loss);
421+export function synthesize(model: ModalModel, pl: Plate, p: SimParams, withMovie = true): Synthesis {
422+ const sigma = modalRates(model);
411423 const g = strikeGains(model, pl, p);
412424 const contact = simulateContact(model, sigma, g, p, p.fps);
413425 const { weight, delay, maxDelay } = listenerWeights(pl, p);
src/modal/worker.tsmodified+13−3View file
@@ -1,16 +1,18 @@
11 // The Web Worker of the linear model. A worker holds the linear system for
22 // the eigenvalue jobs of modes.ts (inertia counts and slices) and, later, a
3-// slice of the modes for the synthesis parts of synth.ts. pool.ts drives a
4-// set of these.
3+// slice of the modes for the per-mode losses of damping.ts and the synthesis
4+// parts of synth.ts. pool.ts drives a set of these.
55
66 import { type LinearSystem } from "../model/assemble.ts";
77 import { runJob, type Job } from "./modes.ts";
8+import { type DampingPlate, radiationRates, thermoelasticRates } from "./damping.ts";
89 import { type ContactResult, type ModalModel, modalTail, movieFrames, nodalAccelerationTable, radiationGains } from "./synth.ts";
910
1011 export type WorkerRequest =
1112 | { type: "system"; n: number; bw: number; band: Float64Array; mass: Float64Array }
1213 | { type: "job"; id: number; job: Job }
13- | { type: "model"; model: ModalModel }
14+ | { type: "model"; model: ModalModel; plate: DampingPlate | null }
15+ | { type: "rates"; id: number }
1416 | {
1517 type: "synth";
1618 id: number;
@@ -29,11 +31,13 @@ export type WorkerRequest =
2931 export type WorkerResponse =
3032 | { type: "ready" }
3133 | { type: "job"; id: number; result: unknown }
34+ | { type: "rates"; id: number; rad: Float64Array; te: Float64Array }
3235 | { type: "synth"; id: number; table: Float32Array | null; audio: Float32Array | null; frames: Float32Array | null }
3336 | { type: "error"; id: number; message: string };
3437
3538 let sys: LinearSystem | null = null;
3639 let model: ModalModel | null = null;
40+let plate: DampingPlate | null = null;
3741
3842 self.onmessage = (ev: MessageEvent<WorkerRequest>) => {
3943 const msg = ev.data;
@@ -48,7 +52,13 @@ self.onmessage = (ev: MessageEvent<WorkerRequest>) => {
4852 post({ type: "job", id: msg.id, result }, transfer);
4953 } else if (msg.type === "model") {
5054 model = msg.model;
55+ plate = msg.plate;
5156 post({ type: "ready" });
57+ } else if (msg.type === "rates") {
58+ if (!model || !plate) throw new Error("no model or plate");
59+ const rad = radiationRates(model, plate, 0, model.count);
60+ const te = thermoelasticRates(model, plate, 0, model.count);
61+ post({ type: "rates", id: msg.id, rad, te }, [rad.buffer as ArrayBuffer, te.buffer as ArrayBuffer]);
5262 } else if (msg.type === "synth") {
5363 if (!model) throw new Error("no model");
5464 const c = msg.contact;
src/model/reference.tsmodified+2−2View file
@@ -11,6 +11,7 @@ import {
1111 type Material,
1212 type Plate,
1313 type StandRates,
14+ ROOM_TEMPERATURE,
1415 bendingStiffness,
1516 cellThickness,
1617 insideCell,
@@ -37,8 +38,7 @@ export function fields(n: number): Fields {
3738 export function uniformPlate(g: Grid, mat: Material, h: number, w0?: Float64Array): Plate {
3839 const hs = new Float64Array(g.n);
3940 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (insideNode(g, i, j)) hs[j * g.nx + i] = h;
40- const freqs = new Float64Array([10, 40000]);
41- return { name: "uniform", Lx: (g.nx - 1) * g.hx, Ly: (g.ny - 1) * g.hy, grid: g, mat, h: hs, w0: w0 ?? new Float64Array(g.n), stand: null, loss: { freqs, t60: new Float64Array([1e9, 1e9]) } };
41+ return { name: "uniform", Lx: (g.nx - 1) * g.hx, Ly: (g.ny - 1) * g.hy, grid: g, mat, h: hs, w0: w0 ?? new Float64Array(g.n), stand: null, temperature: ROOM_TEMPERATURE };
4242 }
4343
4444 const G3 = 1 / Math.sqrt(3);
src/mruntime.tsmodified+5−25View file
@@ -1,9 +1,9 @@
11 // The only .m file the simulator puts in scope of a design script. It is not
22 // a design helper: the trailer that evaluatePlate appends calls it to sample
3-// a field (a number or a function handle) on the grid matrices, or the decay
4-// time on a row of frequencies, and to report a bad return value in the
5-// script's own terms. Design scripts define whatever shape functions they
6-// need themselves, at the bottom of the script.
3+// a field (a number or a function handle) on the grid matrices, and to
4+// report a bad return value in the script's own terms. Design scripts define
5+// whatever shape functions they need themselves, at the bottom of the
6+// script.
77
88 import type { WorkspaceFile } from "numbl";
99
@@ -30,24 +30,4 @@ end
3030 end
3131 `;
3232
33-const T60 = `function v = cymbalsim_t60(f, freqs)
34-% Samples the decay time on a row of frequencies when it is a function handle.
35-try
36- v = f(freqs);
37-catch e
38- error('design.loss.T60: %s', e.message);
39-end
40-if isscalar(v)
41- v = double(v) * ones(size(freqs));
42-elseif numel(v) ~= numel(freqs)
43- error('design.loss.T60 returned %d values for %d frequencies; use elementwise operators (./, .^)', numel(v), numel(freqs));
44-else
45- v = double(v);
46-end
47-end
48-`;
49-
50-export const RUNTIME_FILES: WorkspaceFile[] = [
51- { name: "cymbalsim_field.m", source: SAMPLE },
52- { name: "cymbalsim_t60.m", source: T60 },
53-];
33+export const RUNTIME_FILES: WorkspaceFile[] = [{ name: "cymbalsim_field.m", source: SAMPLE }];
src/params.tsmodified+17−11View file
@@ -1,6 +1,6 @@
11 // Simulation parameters (strike, listener, numerics), the grid with its mask,
22 // the plate as evaluated on the grid (material, thickness field, stress-free
3-// shape, stand, losses), and the derived quantities both solvers need. All
3+// shape, stand, loss constants), and the derived quantities both solvers need. All
44 // values are SI. The design itself is produced by the script in design.ts.
55
66 export interface SimParams {
@@ -171,12 +171,26 @@ export function countInside(g: Grid): number {
171171
172172 // ---------- Material, stand, losses, plate ----------
173173
174+/** Thermal constants of the metal, from which the thermoelastic damping follows (loss.ts). */
175+export interface Thermal {
176+ alpha: number; // linear thermal expansion coefficient [1/K]
177+ cp: number; // specific heat [J/(kg K)]
178+ k: number; // thermal conductivity [W/(m K)]
179+}
180+
174181 export interface Material {
175182 E: number; // Young's modulus [Pa]
176183 nu: number; // Poisson ratio
177184 rho: number; // density [kg/m^3]
185+ /** Loss factor of the metal's internal friction (hysteretic, frequency independent); 0 for none. */
186+ eta: number;
187+ /** Thermal constants, or null when the thermoelastic damping is left out. */
188+ thermal: Thermal | null;
178189 }
179190
191+/** Temperature of the plate [K], which enters the thermoelastic damping. */
192+export const ROOM_TEMPERATURE = 293;
193+
180194 /** A felt pad of the stand: an annulus of felt between the plate and the stand (or the wing nut). */
181195 export interface Pad {
182196 innerRadius: number; // [m]
@@ -199,15 +213,6 @@ export interface Stand {
199213 preload: number; // compression of the top pad at rest [m]
200214 }
201215
202-/**
203- * Decay times of the plate as a function of frequency, sampled by the design
204- * script on a fixed log-spaced set of frequencies (see loss.ts).
205- */
206-export interface Loss {
207- freqs: Float64Array;
208- t60: Float64Array;
209-}
210-
211216 /** The plate as evaluated on the grid. */
212217 export interface Plate {
213218 name: string;
@@ -220,7 +225,8 @@ export interface Plate {
220225 /** Stress-free height at the nodes [m], mean plane removed; zero outside the plate. */
221226 w0: Float64Array;
222227 stand: Stand | null;
223- loss: Loss;
228+ /** Temperature [K], for the thermoelastic damping. */
229+ temperature: number;
224230 }
225231
226232 export function bendingStiffness(m: Material, h: number): number {
src/presets.tsmodified+42−11View file
@@ -16,7 +16,10 @@ const flatRectangle = `% A flat rectangular sheet of stainless steel, hanging fr
1616 % point, measured from the center of the bounding box; the handles are
1717 % applied to whole grid matrices, so use elementwise operators (.*, ./, .^)
1818 % inside them. Anything else a script needs it defines itself, as MATLAB
19-% functions at the bottom; the other presets show how.
19+% functions at the bottom; the other presets show how. The decay of the
20+% sound is not prescribed: it follows from the loss factor of the metal, its
21+% thermal constants (thermoelastic damping), the radiation of each mode into
22+% the air, and the felt of a stand.
2023 inch = 0.0254; % m
2124 mm = 1e-3; % m
2225
@@ -24,18 +27,22 @@ design.name = 'Flat steel rectangle';
2427 design.material.E = 170e9; % Young's modulus [Pa]
2528 design.material.nu = 0.29; % Poisson ratio
2629 design.material.rho = 7900; % density [kg/m^3]
30+design.material.lossFactor = 3e-4; % internal friction of the metal (typical published value)
31+design.material.alpha = 17e-6; % thermal expansion [1/K] } for the thermoelastic damping
32+design.material.cp = 500; % specific heat [J/(kg K)] }
33+design.material.k = 16; % conductivity [W/(m K)] }
2734 design.bounds.width = 12*inch; % bounding box [m]
2835 design.bounds.height = 18*inch;
2936 design.thickness = 0.79*mm;
3037 design.stand = []; % free plate
31-design.loss.T60 = @(f) 3 ./ (1 + f/1500); % decay time [s] as a function of frequency [Hz]
3238 `;
3339
3440 const aluminumSquare = `% A flat 12 inch square of 6061-T6 aluminum, 1.02 mm thick, with a center
3541 % hole on a cymbal stand: felt pads below and above the plate around the
3642 % rod. Aluminum is a third the density of steel and a third as stiff, so
3743 % the bending frequencies land near those of a steel plate of the same size
38-% and thickness; the material damping is low, hence the long decay.
44+% and thickness. Its high thermal expansion and conductivity give it the
45+% strongest thermoelastic damping of the three metals here.
3946 inch = 0.0254; mm = 1e-3;
4047 W = 12*inch; hole = 13*mm;
4148
@@ -43,6 +50,10 @@ design.name = 'Aluminum 6061-T6 square';
4350 design.material.E = 68.9e9; % Young's modulus [Pa]
4451 design.material.nu = 0.33; % Poisson ratio
4552 design.material.rho = 2700; % density [kg/m^3]
53+design.material.lossFactor = 1e-4; % internal friction of the metal (typical published value)
54+design.material.alpha = 23.6e-6; % thermal expansion [1/K] } for the thermoelastic damping
55+design.material.cp = 896; % specific heat [J/(kg K)] }
56+design.material.k = 167; % conductivity [W/(m K)] }
4657 design.bounds.width = W;
4758 design.bounds.height = W;
4859 design.inside = @(x, y, r, theta) r >= hole/2;
@@ -54,7 +65,6 @@ design.stand.pad.E = 0.5e6; % Pa, felt in compression
5465 design.stand.pad.lossFactor = 0.3;
5566 design.stand.topPad = true;
5667 design.stand.preload = 1*mm; % compression of the top pad by the wing nut
57-design.loss.T60 = @(f) 4 ./ (1 + f/2000);
5868 `;
5969
6070 const bowedRectangle = `% The same steel sheet with a 3 mm parabolic bow across its width, the
@@ -68,12 +78,15 @@ design.name = 'Bowed steel rectangle';
6878 design.material.E = 170e9;
6979 design.material.nu = 0.29;
7080 design.material.rho = 7900;
81+design.material.lossFactor = 3e-4; % internal friction of the metal (typical published value)
82+design.material.alpha = 17e-6; % thermal expansion [1/K] } for the thermoelastic damping
83+design.material.cp = 500; % specific heat [J/(kg K)] }
84+design.material.k = 16; % conductivity [W/(m K)] }
7185 design.bounds.width = W;
7286 design.bounds.height = H;
7387 design.thickness = 0.79*mm;
7488 design.profile = @(x, y, r, theta) bow(x, W, 3*mm, 0.3); % stress-free height [m]
7589 design.stand = [];
76-design.loss.T60 = @(f) 3 ./ (1 + f/1500);
7790
7891 function v = bow(s, span, sag, offset)
7992 % A parabola of the given sag across a span centered on 0; s runs along the
@@ -95,6 +108,10 @@ design.name = 'Flat bronze disc';
95108 design.material.E = 110e9; % bronze
96109 design.material.nu = 0.34;
97110 design.material.rho = 8800;
111+design.material.lossFactor = 2e-4; % internal friction of the metal (typical published value)
112+design.material.alpha = 18e-6; % thermal expansion [1/K] } for the thermoelastic damping
113+design.material.cp = 370; % specific heat [J/(kg K)] }
114+design.material.k = 40; % conductivity [W/(m K)] }
98115 design.bounds.width = D;
99116 design.bounds.height = D;
100117 design.inside = @(x, y, r, theta) (r <= D/2) & (r >= hole/2);
@@ -106,7 +123,6 @@ design.stand.pad.E = 0.5e6; % Pa, felt in compression
106123 design.stand.pad.lossFactor = 0.3;
107124 design.stand.topPad = true;
108125 design.stand.preload = 1*mm; % compression of the top pad by the wing nut
109-design.loss.T60 = @(f) 4 ./ (1 + f/2000);
110126 `;
111127
112128 const crash16 = `% A 16 inch flat crash: bronze, lathed to a taper from 1.2 mm at the
@@ -121,6 +137,10 @@ design.name = '16 inch flat crash';
121137 design.material.E = 110e9;
122138 design.material.nu = 0.34;
123139 design.material.rho = 8800;
140+design.material.lossFactor = 2e-4; % internal friction of the metal (typical published value)
141+design.material.alpha = 18e-6; % thermal expansion [1/K] } for the thermoelastic damping
142+design.material.cp = 370; % specific heat [J/(kg K)] }
143+design.material.k = 40; % conductivity [W/(m K)] }
124144 design.bounds.width = D;
125145 design.bounds.height = D;
126146 design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
@@ -132,7 +152,6 @@ design.stand.pad.E = 0.5e6;
132152 design.stand.pad.lossFactor = 0.3;
133153 design.stand.topPad = true;
134154 design.stand.preload = 1*mm;
135-design.loss.T60 = [100 5; 1000 3.5; 5000 1.2; 15000 0.4]; % table of [Hz seconds], interpolated
136155
137156 function v = ramp(a, b, x)
138157 % 0 at x = a, 1 at x = b, clamped to that range outside.
@@ -151,6 +170,10 @@ design.name = '20 inch flat ride';
151170 design.material.E = 110e9;
152171 design.material.nu = 0.34;
153172 design.material.rho = 8800;
173+design.material.lossFactor = 2e-4; % internal friction of the metal (typical published value)
174+design.material.alpha = 18e-6; % thermal expansion [1/K] } for the thermoelastic damping
175+design.material.cp = 370; % specific heat [J/(kg K)] }
176+design.material.k = 40; % conductivity [W/(m K)] }
154177 design.bounds.width = D;
155178 design.bounds.height = D;
156179 design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
@@ -162,7 +185,6 @@ design.stand.pad.E = 0.5e6;
162185 design.stand.pad.lossFactor = 0.3;
163186 design.stand.topPad = true;
164187 design.stand.preload = 0.3*mm;
165-design.loss.T60 = @(f) 8 ./ (1 + f/2500);
166188
167189 function v = ramp(a, b, x)
168190 v = min(max((x - a) ./ (b - a), 0), 1);
@@ -184,6 +206,10 @@ design.name = '16 inch crash with bell';
184206 design.material.E = 110e9;
185207 design.material.nu = 0.34;
186208 design.material.rho = 8800;
209+design.material.lossFactor = 2e-4; % internal friction of the metal (typical published value)
210+design.material.alpha = 18e-6; % thermal expansion [1/K] } for the thermoelastic damping
211+design.material.cp = 370; % specific heat [J/(kg K)] }
212+design.material.k = 40; % conductivity [W/(m K)] }
187213 design.bounds.width = D;
188214 design.bounds.height = D;
189215 design.inside = @(x, y, r, theta) (r <= R) & (r >= hole/2);
@@ -196,7 +222,6 @@ design.stand.pad.E = 0.5e6;
196222 design.stand.pad.lossFactor = 0.3;
197223 design.stand.topPad = true;
198224 design.stand.preload = 1*mm;
199-design.loss.T60 = [100 5; 1000 3.5; 5000 1.2; 15000 0.4];
200225
201226 function v = ramp(a, b, x)
202227 v = min(max((x - a) ./ (b - a), 0), 1);
@@ -219,12 +244,15 @@ design.name = 'Steel ellipse';
219244 design.material.E = 170e9;
220245 design.material.nu = 0.29;
221246 design.material.rho = 7900;
247+design.material.lossFactor = 3e-4; % internal friction of the metal (typical published value)
248+design.material.alpha = 17e-6; % thermal expansion [1/K] } for the thermoelastic damping
249+design.material.cp = 500; % specific heat [J/(kg K)] }
250+design.material.k = 16; % conductivity [W/(m K)] }
222251 design.bounds.width = A;
223252 design.bounds.height = B;
224253 design.inside = @(x, y, r, theta) (x/(A/2)).^2 + (y/(B/2)).^2 <= 1;
225254 design.thickness = 1.0*mm;
226255 design.stand = [];
227-design.loss.T60 = @(f) 3 ./ (1 + f/1500);
228256 `;
229257
230258 const thunderSheet = `% A thunder sheet: a large thin steel sheet hung from one corner. Its modes
@@ -236,11 +264,14 @@ design.name = 'Thunder sheet';
236264 design.material.E = 200e9;
237265 design.material.nu = 0.3;
238266 design.material.rho = 7850;
267+design.material.lossFactor = 3e-4; % internal friction of the metal (typical published value)
268+design.material.alpha = 12e-6; % thermal expansion [1/K] } for the thermoelastic damping
269+design.material.cp = 460; % specific heat [J/(kg K)] }
270+design.material.k = 50; % conductivity [W/(m K)] }
239271 design.bounds.width = 24*inch;
240272 design.bounds.height = 36*inch;
241273 design.thickness = 0.5*mm;
242274 design.stand = [];
243-design.loss.T60 = @(f) 2.5 ./ (1 + f/1000);
244275 `;
245276
246277 export const PRESETS: Preset[] = [
src/sim/selftest.tsmodified+3−3View file
@@ -4,7 +4,7 @@
44
55 import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode, standRates } from "../params.ts";
66 import { accelerations, addDamping, fields, fullLoss } from "../model/reference.ts";
7-import { fitSigma, decayRate } from "../loss.ts";
7+import { fitSigma, inPlaneRate } from "../loss.ts";
88 import { Simulation } from "./simulation.ts";
99
1010 export interface SelfTestResult {
@@ -36,9 +36,9 @@ export function smoothField(pl: Plate, amplitude: number, seed: number): Float64
3636 export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParams): Promise<SelfTestResult> {
3737 const g = pl.grid, n = g.n;
3838 const hm = meanThickness(pl);
39- const fit = fitSigma(pl.loss, pl.mat, hm);
39+ const fit = fitSigma(pl);
4040 const sr = standRates(pl);
41- const loss = fullLoss(fit, sr, decayRate(pl.loss, 5000));
41+ const loss = fullLoss(fit, sr, inPlaneRate(pl));
4242 const sim = new Simulation(device, pl, params, loss);
4343 const amp = 2.5 * hm;
4444 const f = { u: smoothField(pl, 0.02 * amp, 1), v: smoothField(pl, 0.02 * amp, 2), w: smoothField(pl, amp, 3) };
src/tutorial/steps.tsmodified+41−22View file
@@ -11,7 +11,7 @@
1111
1212 import { type Design, profileStats } from "../design.ts";
1313 import { type Plate, type SimParams, bendingConstant, bendingStiffness, inPlaneSpeed, massProperties, meanThickness, nodeX, nodeY, standRates, strikePoint, thicknessRange, dispersionFrequency, CFL, listenerPosition, GRAVITY } from "../params.ts";
14-import { LN1000, t60At } from "../loss.ts";
14+import { LN1000, plateGeometry, estimatedRate, thermoelasticRelaxation, rateText } from "../loss.ts";
1515 import { type LinearSystem } from "../model/assemble.ts";
1616 import { maxAngularFrequency } from "../model/reference.ts";
1717 import { type ModeSet, defaultSliceSize } from "../modal/modes.ts";
@@ -168,7 +168,7 @@ function modeLabels(pl: Plate, model: ModalModel, count: number): string[] {
168168
169169 /** Magnitude of the radiation gain of each of the first `count` modes for the current listener, times omega^2 (a pressure per unit modal amplitude). */
170170 function pressureGains(pl: Plate, model: ModalModel, p: SimParams, count: number): number[] {
171- const sigma = modalRates(model, pl.loss);
171+ const sigma = modalRates(model);
172172 const { weight, delay } = listenerWeights(pl, p);
173173 const n = Math.min(count, model.count);
174174 const re = new Array<number>(n).fill(0), im = new Array<number>(n).fill(0);
@@ -199,9 +199,10 @@ function compareTable(a: MemoModes, b: MemoModes, n: number, expected?: number):
199199 // ---------- Scripts used by the actions ----------
200200
201201 /** The flat steel rectangle of the preset with some of its numbers changed. */
202-export function rectangleScript(o: { thicknessMm?: number; widthIn?: number; heightIn?: number; t60?: string } = {}): string {
202+export function rectangleScript(o: { thicknessMm?: number; widthIn?: number; heightIn?: number; eta?: number; thermal?: boolean } = {}): string {
203203 const t = o.thicknessMm ?? 0.79, w = o.widthIn ?? 12, h = o.heightIn ?? 18;
204- const t60 = o.t60 ?? "@(f) 3 ./ (1 + f/1500)";
204+ const eta = o.eta ?? 3e-4;
205+ const thermal = o.thermal ?? true;
205206 return `% The flat steel rectangle of the tutorial (a variant of the preset).
206207 inch = 0.0254; mm = 1e-3;
207208
@@ -211,9 +212,9 @@ design.material.nu = 0.29;
211212 design.material.rho = 7900;
212213 design.bounds.width = ${w}*inch;
213214 design.bounds.height = ${h}*inch;
214-design.thickness = ${t}*mm;
215+design.material.lossFactor = ${eta};
216+${thermal ? "design.material.alpha = 17e-6;\ndesign.material.cp = 500;\ndesign.material.k = 16;\n" : ""}design.thickness = ${t}*mm;
215217 design.stand = [];
216-design.loss.T60 = ${t60};
217218 `;
218219 }
219220
@@ -423,27 +424,45 @@ ${liveText ? live(liveText) : need("Compute the modes of a design to see the siz
423424 body: (c) => {
424425 let tab = "";
425426 if (c.plate && c.model) {
426- const sigma = modalRates(c.model, c.plate.loss);
427+ const sigma = modalRates(c.model);
428+ const t60 = rateText;
427429 const rows: string[][] = [];
428- for (let m = 0; m < Math.min(6, c.model.count); m++) {
429- const f = c.model.freqs[m], s = sigma[m], t60 = LN1000 / s;
430- rows.push([String(m + 1), f.toFixed(1), t60.toFixed(2), s.toFixed(2), ((Math.PI * f) / s).toFixed(0), (f * t60).toFixed(0)]);
430+ for (let m = 0; m < Math.min(8, c.model.count); m++) {
431+ rows.push([String(m + 1), c.model.freqs[m].toFixed(1), t60(sigma[m]), t60(c.model.sigmaRad[m]), t60(c.model.sigmaMat[m]), t60(c.model.sigmaTe[m]), t60(c.model.padSigma[m])]);
431432 }
432- tab = table(["Mode", "f (Hz)", "T60 (s)", "σ (1/s)", "Q", "Cycles in T60"], rows);
433+ tab = table(["Mode", "f (Hz)", "T60 (s)", "air alone", "metal alone", "thermal alone", "felt alone"], rows);
434+ }
435+ let est = "";
436+ if (c.plate) {
437+ const pl = c.plate, geom = plateGeometry(pl);
438+ const th = pl.mat.thermal;
439+ const rel = thermoelasticRelaxation(pl.mat, geom.hMean, pl.temperature);
440+ const t60 = rateText;
441+ const e = [100, 1000, 5000].map((f) => estimatedRate(pl, geom, f));
442+ est = `This design: loss factor $\\eta$ = ${sig(pl.mat.eta, 2)}, so internal friction alone gives $T_{60}$ = ${t60(e[0].material)} / ${t60(e[1].material)} / ${t60(e[2].material)} s at 100 Hz / 1 kHz / 5 kHz.<br>` +
443+ (th ? `Thermoelastic: $\\Delta$ = ${sig(rel.delta, 2)}, $\\tau$ = ${sig(rel.tau * 1e3, 2)} ms at the mean thickness, so the loss peaks near ${(1 / (2 * Math.PI * rel.tau)).toFixed(0)} Hz at $\\eta_{te}$ = ${sig(rel.delta / 2, 2)}; alone it gives $T_{60}$ = ${t60(e[0].thermoelastic)} / ${t60(e[1].thermoelastic)} / ${t60(e[2].thermoelastic)} s.<br>` : "No thermal constants in the script, so no thermoelastic damping.<br>") +
444+ `Radiation (average efficiency, coincidence at ${(geom.fc / 1e3).toFixed(1)} kHz): alone $T_{60}$ = ${t60(e[0].radiation)} / ${t60(e[1].radiation)} / ${t60(e[2].radiation)} s. All together: ${t60(e[0].total)} / ${t60(e[1].total)} / ${t60(e[2].total)} s.` +
445+ (tab ? tab : "<br>Compute the modes to see the exact values per mode.");
433446 }
434- const loss = c.plate ? c.plate.loss : null;
435447 return `
436-<p>A real plate loses energy: to the air, to internal friction in the metal, to its support. In the linear model each mode is a damped oscillator,</p>
448+<p>A real plate loses energy. In the linear model each mode is a damped oscillator,</p>
437449 $$ \\ddot q + 2\\sigma\\dot q + \\omega^2 q = 0, \\qquad q(t) = e^{-\\sigma t}\\cos(\\omega_d t), \\quad \\omega_d = \\sqrt{\\omega^2 - \\sigma^2}, $$
438-<p>whose amplitude falls by the factor $e$ in the time $1/\\sigma$. Musicians and acousticians state a decay by the time $T_{60}$ over which the level drops by 60 decibels, a factor of 1000 in amplitude: $e^{-\\sigma T_{60}} = 10^{-3}$, so</p>
439-$$ \\sigma = \\frac{\\ln 1000}{T_{60}} = \\frac{6.91}{T_{60}}. $$
440-<p>The design script gives $T_{60}$ as a function of frequency, and the model applies it to every mode at the mode's own frequency; the shorter decay of the high partials is what turns the bright attack of a struck plate into the low ring. The quality factor $Q = \\omega/(2\\sigma) = \\pi f T_{60}/\\ln 1000$ counts, roughly, the number of cycles before the amplitude has fallen by a factor of 20; the modes of a plate have $Q$ in the hundreds to thousands, which is why each partial is a sharp line in the spectrogram. A decay much slower than the period, $\\sigma \\ll \\omega$, barely shifts the frequency, and the model caps $\\sigma$ below $\\omega$ so that no mode is overdamped. The values are estimates: nothing here is measured.</p>
441-${loss ? live(`This design: $T_{60}$ = ${t60At(loss, 100).toFixed(2)} s at 100 Hz, ${t60At(loss, 1000).toFixed(2)} s at 1 kHz, ${t60At(loss, 5000).toFixed(2)} s at 5 kHz.${tab || "<br>Compute the modes to see the values per mode."}`) : need("Load a design to see its decay times.")}
442-<p>Try it: give the plate a constant decay of 1 s at every frequency, synthesize, and listen; the high partials, which the preset lets die in a fraction of a second, now ring as long as the low ones.</p>`;
450+<p>whose amplitude falls by the factor $e$ in the time $1/\\sigma$. Musicians and acousticians state a decay by the time $T_{60}$ over which the level drops by 60 decibels, a factor of 1000 in amplitude: $e^{-\\sigma T_{60}} = 10^{-3}$, so $\\sigma = \\ln(1000)/T_{60}$. The quality factor $Q = \\omega/(2\\sigma)$ counts, roughly, the cycles before the amplitude has fallen by a factor of 20; the modes of a plate have $Q$ in the hundreds to thousands, which is why each partial is a sharp line in the spectrogram.</p>
451+<p>Nothing in the script prescribes the decay. The rate of each mode is the sum of four mechanisms, each with its own physics.</p>
452+<p><i>Internal friction of the metal.</i> A hysteretic loss: each cycle dissipates the fraction $2\\pi\\eta$ of the stored energy, so a mode decays at $\\sigma = \\eta\\omega/2$. The loss factor $\\eta$ is a material constant like $E$ and $\\rho$, measured rather than derived, of order $10^{-4}$ for bronze, steel and aluminum; hammering and alloying change it. Alone it would let the low modes ring for minutes.</p>
453+<p><i>Thermoelastic damping.</i> Bending warms the compressed face and cools the stretched one, and the heat that flows across the thickness is lost. Zener's analysis gives the loss factor</p>
454+$$ \\eta_{te} = \\Delta\\,\\frac{\\omega\\tau}{1 + (\\omega\\tau)^2}, \\qquad \\Delta = \\frac{E\\alpha^2 T}{\\rho c_p}\\,\\frac{1+\\nu}{1-\\nu}, \\qquad \\tau = \\frac{h^2 \\rho c_p}{\\pi^2 k}, $$
455+<p>from the thermal expansion $\\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$; $\\tau$ is the time heat takes to cross the thickness, and the loss peaks at $\\omega\\tau = 1$. The model applies it to the bending energy of each mode with the local thickness. For a millimeter of metal the peak lies below 100 Hz and the effect is small; it matters for thin aluminum.</p>
456+<p><i>Radiation.</i> The sound itself carries energy away. The power a mode radiates follows from its shape through the same Rayleigh integral that gives the listener's pressure, summed over all directions (the page <i>From the plate to the ear</i>), and the mode decays at $\\sigma = P/(2E)$. A plate with free edges radiates surprisingly little below its coincidence frequency, where the bending wavelength equals the wavelength of sound: its modes have no net volume velocity, since they are orthogonal to the rigid motions, and a free edge, where the moment and the shear vanish, lacks the uncancelled strip that makes a simply supported edge radiate. The efficiency rises steeply toward coincidence, so the radiation overtakes the metal only in the top octaves, and it is there that the high partials die first. The column <i>air alone</i> below shows how weak it is for the low modes.</p>
457+<p><i>The felt of a stand</i> adds its own loss to the modes that move the center (the page on the stand).</p>
458+${est ? live(est) : need("Load a design to see its losses.")}
459+<p>Try it: raise the loss factor to $10^{-3}$, the value of a heavily worked alloy, and hear the ring shorten; remove it and hear the low modes ring on with only the air to stop them.</p>`;
443460 },
444461 actions: [
445- { label: "Constant decay of 1 s", run: async (h) => { await h.loadScript(rectangleScript({ t60: "1" })); await h.synthesize(); } },
446- { label: "The preset's decay", run: async (h) => { await h.loadPreset("flat-rectangle"); await h.synthesize(); } },
462+ { label: "Compute the modes", run: (h) => h.computeModes() },
463+ { label: "Loss factor 1e-3", run: async (h) => { await h.loadScript(rectangleScript({ eta: 1e-3 })); await h.synthesize(); } },
464+ { label: "No internal friction", run: async (h) => { await h.loadScript(rectangleScript({ eta: 0 })); await h.synthesize(); } },
465+ { label: "The preset", run: async (h) => { await h.loadPreset("flat-rectangle"); await h.synthesize(); } },
447466 { label: "Play", run: (h) => h.play() },
448467 ],
449468 },
@@ -553,7 +572,7 @@ $$ p(t) = \\frac{\\rho_{\\mathrm{air}}}{2\\pi}\\sum_k \\frac{\\Delta A_k}{r_k}\\
553572 $$ G_m = \\frac{\\rho_{\\mathrm{air}}}{2\\pi}\\sum_k \\frac{\\Delta A_k}{r_k}\\,\\phi_m(k)\\, e^{-s_m r_k/c}, $$
554573 <p>and the pressure is</p>
555574 $$ p(t) = \\sum_m \\mathrm{Re}\\left[G_m s_m^2 A_m e^{s_m t}\\right], $$
556-<p>the same recurrence as before with different amplitudes. The gain contains the directivity of the mode: its positive and negative regions partly cancel, most completely on the axis in front of the plate for a mode with equal lobes, and less so from an angle, where the delays from the lobes differ. Moving the listener recomputes only the gains, which is why the sound follows a move at once.</p>
575+<p>the same recurrence as before with different amplitudes. The gain contains the directivity of the mode: its positive and negative regions partly cancel, most completely on the axis in front of the plate for a mode with equal lobes, and less so from an angle, where the delays from the lobes differ. Moving the listener recomputes only the gains, which is why the sound follows a move at once. The same sum taken over every direction gives the power the mode radiates, and with it the radiation part of the mode's decay on the page <i>Decay</i>.</p>
557576 <p>This is the one place where the model steps outside the plate, and it is approximate: the real plate has no wall, so at low frequencies the sound from its back partly cancels that from its front, and the Rayleigh integral overstates the bass.</p>
558577 ${tab ? live(tab) : need("Compute the modes to see the gains of the first modes at the listener.")}
559578 <p>Try it: on axis, the modes with equal lobes nearly vanish; at a grazing angle they return and the pattern changes.</p>`;
@@ -608,7 +627,7 @@ ${live(`Window of ${N} samples at ${fs / 1000} kHz: $\\Delta f$ = ${(fs / N).toF
608627 liveText = `${np} pad${np > 1 ? "s" : ""} of ${mm(pad.thickness, 0)} felt at $E_p$ = ${sig(pad.E / 1e6)} MPa between ${mm(pad.innerRadius, 0)} and ${mm(pad.outerRadius, 0)}; mass $M$ = ${sig(mass * 1e3)} g.<br>$K$ = ${sig(sr.totalKw / 1e3)} kN/m, so $f_b$ = ${sr.bounceHz.toFixed(1)} Hz for a rigid cymbal; $K_\\theta$ = ${sig(Ktheta)} N m, $I$ = ${sig(I)} kg m², so $f_r$ = ${fRock.toFixed(1)} Hz. Static compression of the bottom pad under the weight alone: ${mm(mass * GRAVITY / (sr.totalKw / np), 3)}.<br>First umbrella mode of a uniform free disc of this radius and mean thickness: ${fUmbrella.toFixed(1)} Hz.`;
609628 if (c.model) {
610629 const kinds = modeKinds(pl, c.model, 8);
611- const sigma = modalRates(c.model, pl.loss);
630+ const sigma = modalRates(c.model);
612631 const rows = kinds.map((k, m) => [String(m + 1), c.model!.freqs[m].toFixed(1), k.label, k.share.toFixed(2), (LN1000 / sigma[m]).toFixed(2)]);
613632 liveText += table(["Mode", "f (Hz)", "Kind", "Rigid share", "T60 (s)"], rows);
614633 } else liveText += "<br>Compute the modes to find the rock and the bounce among them.";
src/ui/modes.tsmodified+12−3View file
@@ -2,11 +2,20 @@
22 // decay time), with a filter by frequency; clicking a row selects the mode,
33 // whose shape the host then draws on the plate.
44
5+import { t60Text } from "../loss.ts";
6+
57 export interface ModeRow {
68 freq: number; // [Hz]
7- t60: number; // [s]
9+ /** Decay time of the mode [s], and the decay time each mechanism alone would give (Infinity when it contributes nothing). */
10+ t60: number;
11+ radiation: number;
12+ material: number;
13+ thermoelastic: number;
14+ felt: number;
815 }
916
17+const fmtT60 = t60Text;
18+
1019 export class ModeBrowser {
1120 private info: HTMLDivElement;
1221 private tableWrap: HTMLDivElement;
@@ -43,7 +52,7 @@ export class ModeBrowser {
4352 this.tableWrap.className = "modes-table";
4453 const table = document.createElement("table");
4554 const thead = document.createElement("thead");
46- thead.innerHTML = "<tr><th>#</th><th>f (Hz)</th><th>T60 (s)</th></tr>";
55+ thead.innerHTML = "<tr><th>#</th><th>f (Hz)</th><th title=\"Decay time of the mode\">T60 (s)</th><th title=\"Decay time from radiation into the air alone\">air</th><th title=\"Decay time from the metal's internal friction alone\">metal</th><th title=\"Decay time from thermoelastic damping alone\">thermal</th><th title=\"Decay time from the felt of the stand alone\">felt</th></tr>";
4756 this.tbody = document.createElement("tbody");
4857 table.append(thead, this.tbody);
4958 this.tableWrap.append(table);
@@ -98,7 +107,7 @@ export class ModeBrowser {
98107 tr.dataset.m = String(m);
99108 if (m === this.sel) tr.classList.add("selected");
100109 const f = r.freq < 1000 ? r.freq.toFixed(1) : r.freq.toFixed(0);
101- tr.innerHTML = `<td>${m + 1}</td><td>${f}</td><td>${r.t60.toFixed(2)}</td>`;
110+ tr.innerHTML = `<td>${m + 1}</td><td>${f}</td><td>${fmtT60(r.t60)}</td><td>${fmtT60(r.radiation)}</td><td>${fmtT60(r.material)}</td><td>${fmtT60(r.thermoelastic)}</td><td>${fmtT60(r.felt)}</td>`;
102111 tr.addEventListener("click", () => this.select(this.sel === m ? null : m));
103112 frag.append(tr);
104113 });
tests/assemble.test.tsmodified+1−1View file
@@ -6,7 +6,7 @@ import { makeGrid, withMask, type Plate, nodeMasses, nodeX, nodeY, standRates }
66 import { uniformPlate, accelerations, fields } from "../src/model/reference.ts";
77 import { assemble, stiffnessMatvec, toGrid, fromGrid, isFlatPlate } from "../src/model/assemble.ts";
88
9-const bronze = { E: 110e9, nu: 0.34, rho: 8800 };
9+const bronze = { E: 110e9, nu: 0.34, rho: 8800, eta: 0, thermal: null };
1010
1111 function checkAgainstReference(pl: Plate, fieldsN: 1 | 3, seed = 3): number {
1212 const g = pl.grid;
tests/design.test.tsmodified+17−34View file
@@ -3,7 +3,7 @@ import assert from "node:assert/strict";
33 import { compileDesign, evaluatePlate, DesignError, removePlane, isFlat } from "../src/design.ts";
44 import { PRESETS } from "../src/presets.ts";
55 import { countInside, insideNode, nodeX, nodeY, nodeWeight, standRates, massProperties } from "../src/params.ts";
6-import { t60At, decayRate, fitSigma, fittedRate, LN1000, constantLoss, lossFromFunction, lossFromTable, LOSS_FREQS } from "../src/loss.ts";
6+import { ROOM_TEMPERATURE } from "../src/params.ts";
77
88 const base = `
99 mm = 1e-3;
@@ -20,8 +20,9 @@ test("a minimal script compiles with defaults", () => {
2020 assert.equal(d.profileConst, 0);
2121 assert.equal(d.stand, null);
2222 assert.equal(d.hasInside, false);
23- assert.ok(d.lossConst);
24- assert.ok(Math.abs(t60At(d.lossConst!, 100) / (4 / 1.05) - 1) < 1e-3);
23+ assert.equal(d.material.eta, 0, "no internal friction unless the script gives a loss factor");
24+ assert.equal(d.material.thermal, null, "no thermoelastic damping without the thermal constants");
25+ assert.equal(d.temperature, ROOM_TEMPERATURE);
2526 });
2627
2728 test("missing fields and bad values are reported in the design's terms", () => {
@@ -29,7 +30,7 @@ test("missing fields and bad values are reported in the design's terms", () => {
2930 assert.throws(() => compileDesign(base.replace("design.thickness = 1*mm;", "")), /design.thickness is required/);
3031 assert.throws(() => compileDesign(base.replace("design.thickness = 1*mm;", "design.thickness = -1;")), /positive/);
3132 assert.throws(() => compileDesign(base + "design.inside = 3;"), /design.inside must be a function handle/);
32- assert.throws(() => compileDesign(base + "design.loss.T60 = -2;"), /positive decay time/);
33+ assert.throws(() => compileDesign(base + "design.loss.T60 = 2;"), /design.loss is no longer used/);
3334 assert.throws(() => compileDesign(base + "design.material.nu = 0.6;"), /below 0.5/);
3435 });
3536
@@ -123,36 +124,18 @@ design.stand.pad.E = 1e6; design.stand.pad.lossFactor = 0.2; design.stand.preloa
123124 assert.ok(sr2.totalKw > 0);
124125 });
125126
126-test("losses: constant, table, and handle forms", () => {
127- const c = compileDesign(base + "design.loss.T60 = 2.5;");
128- assert.ok(Math.abs(t60At(c.lossConst!, 3000) - 2.5) < 1e-12);
129- const t = compileDesign(base + "design.loss.T60 = [100 4; 10000 1];");
130- assert.ok(Math.abs(t60At(t.lossConst!, 100) / 4 - 1) < 1e-2);
131- assert.ok(Math.abs(t60At(t.lossConst!, 1000) / 2 - 1) < 1e-3, "log-log midpoint");
132- assert.ok(Math.abs(t60At(t.lossConst!, 20) / 4 - 1) < 1e-9, "clamped below");
133- const h = compileDesign(base + "design.loss.T60 = @(f) 3 ./ (1 + f/1000);");
134- assert.equal(h.lossConst, null);
135- const pl = evaluatePlate(h, 0.02);
136- assert.ok(Math.abs(t60At(pl.loss, 1000) / 1.5 - 1) < 1e-3);
137- assert.equal(pl.loss.freqs.length, LOSS_FREQS.length);
138- assert.ok(Math.abs(decayRate(pl.loss, 1000) / (LN1000 / 1.5) - 1) < 1e-3);
139- const bad = compileDesign(base + "design.loss.T60 = @(f) f - 5000;");
140- assert.throws(() => evaluatePlate(bad, 0.02), /must be positive/);
141-});
142-
143-test("the sigma fit reproduces a pure sigma0 + sigma1 k^2 curve", () => {
144- const mat = { E: 200e9, nu: 0.3, rho: 7850 }, h = 1e-3;
145- const kappa = Math.sqrt((mat.E * h ** 3) / (12 * (1 - mat.nu ** 2)) / (mat.rho * h));
146- const s0 = 1.2, s1 = 4e-4;
147- const loss = lossFromFunction((f) => LN1000 / (s0 + (s1 * 2 * Math.PI * f) / kappa));
148- const fit = fitSigma(loss, mat, h);
149- assert.ok(Math.abs(fit.sigma0 / s0 - 1) < 1e-3 && Math.abs(fit.sigma1 / s1 - 1) < 1e-3, `${fit.sigma0} ${fit.sigma1}`);
150- assert.ok(Math.abs(fittedRate(fit, 2000) / (s0 + (s1 * 2 * Math.PI * 2000) / kappa) - 1) < 1e-3);
151- const flat = fitSigma(constantLoss(2), mat, h);
152- assert.ok(flat.sigma1 >= 0 && flat.sigma0 >= 0);
153- const tab = lossFromTable([1000, 100], [1, 2]);
154- // The kink of the table at 100 Hz is smoothed a little by the sampling on LOSS_FREQS.
155- assert.ok(Math.abs(t60At(tab, 100) / 2 - 1) < 1e-2, "table sorted by frequency");
127+test("the loss constants of the material are read and checked", () => {
128+ const d = compileDesign(base + "design.material.lossFactor = 2e-4; design.material.alpha = 12e-6; design.material.cp = 460; design.material.k = 50; design.temperature = 300;");
129+ assert.equal(d.material.eta, 2e-4);
130+ assert.deepEqual(d.material.thermal, { alpha: 12e-6, cp: 460, k: 50 });
131+ assert.equal(d.temperature, 300);
132+ const pl = evaluatePlate(d, 0.02);
133+ assert.equal(pl.temperature, 300);
134+ assert.equal(pl.mat.eta, 2e-4);
135+ assert.throws(() => compileDesign(base + "design.material.alpha = 12e-6;"), /all three thermal constants/);
136+ assert.throws(() => compileDesign(base + "design.material.lossFactor = -1;"), /at least 0/);
137+ assert.throws(() => compileDesign(base + "design.material.lossFactor = 0.5;"), /must be small/);
138+ assert.throws(() => compileDesign(base + "design.temperature = 0;"), /at least 1/);
156139 });
157140
158141 test("every preset compiles and evaluates on the default grid", () => {
tests/loss.test.tsadded+174−0View file
@@ -0,0 +1,174 @@
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.
7+
8+import { test } from "node:test";
9+import assert from "node:assert/strict";
10+import { makeGrid, withMask, nodeWeight, RHO_AIR, SPEED_OF_SOUND, ROOM_TEMPERATURE, type Material } from "../src/params.ts";
11+import { uniformPlate } from "../src/model/reference.ts";
12+import { assemble } from "../src/model/assemble.ts";
13+import { computeModes, localRunner } from "../src/modal/modes.ts";
14+import { modalModel, modalRates, sliceModel } from "../src/modal/synth.ts";
15+import { dampingPlate, radiationRates, radiationRateDirect, thermoelasticRates, computeDamping, bendingEnergyPerNode } from "../src/modal/damping.ts";
16+import { thermoelasticRelaxation, thermoelasticLossFactor, outlinePerimeter, averageRadiationEfficiency, coincidenceFrequency, fitLine, fitSigma, fittedRate, plateGeometry, estimatedRate, materialRate, LN1000 } from "../src/loss.ts";
17+
18+const steel: Material = { E: 200e9, nu: 0.3, rho: 7850, eta: 3e-4, thermal: { alpha: 12e-6, cp: 460, k: 50 } };
19+const bronze: Material = { E: 110e9, nu: 0.34, rho: 8800, eta: 2e-4, thermal: { alpha: 18e-6, cp: 370, k: 40 } };
20+
21+test("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+});
37+
38+test("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+});
50+
51+test("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+});
67+
68+test("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+});
85+
86+/** A small free steel plate with its modes and per-mode losses. */
87+async 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+}
95+
96+test("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+});
110+
111+test("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+});
122+
123+test("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+ }
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), padSigma: new Float64Array(1), masses: new Float64Array(0) };
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+});
148+
149+test("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+});
tests/modal.test.tsmodified+1−1View file
@@ -10,7 +10,7 @@ import { factorize, solve, inertia } from "../src/modal/banded.ts";
1010 import { tridiagonalEigen, sliceModes } from "../src/modal/lanczos.ts";
1111 import { computeModes, localRunner } from "../src/modal/modes.ts";
1212
13-const steel = { E: 200e9, nu: 0.3, rho: 7850 };
13+const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
1414
1515 /** All eigenvalues of M^{-1/2} K M^{-1/2} by cyclic Jacobi (small systems only). */
1616 function denseEigenvalues(sys: LinearSystem): number[] {
tests/reference.test.tsmodified+1−1View file
@@ -8,7 +8,7 @@ import { uniformPlate, accelerations, addDamping, potentialEnergy, fields, maxAn
88 import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
99 import { fullLoss } from "../src/model/reference.ts";
1010
11-const steel = { E: 200e9, nu: 0.3, rho: 7850 };
11+const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
1212
1313 function rnd(seed: number): () => number {
1414 let s = seed >>> 0;
tests/synth.test.tsmodified+8−7View file
@@ -8,17 +8,16 @@ import { makeGrid, defaultParams, type SimParams, listenerTerms, RHO_AIR, SPEED_
88 import { uniformPlate } from "../src/model/reference.ts";
99 import { assemble } from "../src/model/assemble.ts";
1010 import { computeModes, localRunner } from "../src/modal/modes.ts";
11-import { constantLoss, LN1000 } from "../src/loss.ts";
11+import { LN1000 } from "../src/loss.ts";
1212 import {
1313 modalModel, modalRates, strikeGains, simulateContact, listenerWeights, radiationGains, tailStart, nodalAccelerationTable, attackPressure, modalTail, movieFrames, synthesize, sliceModel, sliceContact, type ModalModel,
1414 } from "../src/modal/synth.ts";
1515
16-const steel = { E: 200e9, nu: 0.3, rho: 7850 };
16+const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
1717
1818 async function smallModel(): Promise<{ model: ModalModel; pl: ReturnType<typeof uniformPlate> }> {
1919 const g = makeGrid(0.2, 0.15, 0.0125);
2020 const pl = uniformPlate(g, steel, 1e-3);
21- pl.loss = constantLoss(1.0);
2221 const sys = assemble(pl);
2322 const modes = await computeModes(sys, 4000, localRunner(sys), { sliceSize: 20 });
2423 return { model: modalModel(pl, sys, modes), pl };
@@ -29,7 +28,9 @@ test("a single mode under an impulse follows the closed-form damped sinusoid", a
2928 const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.1, fps: 100 };
3029 // Keep only mode 0 by slicing the model.
3130 const one = sliceModel(model, 0, 1);
32- const sigma = modalRates(one, pl.loss);
31+ // Give the mode a decay time of exactly 1 s through the material term alone.
32+ one.sigmaMat[0] = LN1000; one.sigmaTe[0] = 0; one.sigmaRad[0] = 0; one.padSigma[0] = 0;
33+ const sigma = modalRates(one);
3334 const g = strikeGains(one, pl, p);
3435 const c = simulateContact(one, sigma, g, p, p.fps);
3536 assert.ok(c.peakForce > 0 && c.contactTime > 0);
@@ -86,7 +87,7 @@ test("the radiation gain of a uniform piston equals the analytic Rayleigh sum",
8687 test("the attack and the tail join, and the sliced parts sum to the whole", async () => {
8788 const { model, pl } = await smallModel();
8889 const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.2, fps: 120, strikeVelocity: 2 };
89- const whole = synthesize(model, pl, pl.loss, p, true);
90+ const whole = synthesize(model, pl, p, true);
9091 const nSamples = Math.round(p.duration * p.fs);
9192 assert.equal(whole.audio.length, nSamples);
9293 // No discontinuity at the tail start beyond the local sample-to-sample variation.
@@ -99,7 +100,7 @@ test("the attack and the tail join, and the sliced parts sum to the whole", asyn
99100 const rms = (a: Float32Array) => Math.sqrt(a.reduce((s, v) => s + v * v, 0) / a.length);
100101 assert.ok(rms(whole.audio.subarray(nSamples - 960)) < rms(whole.audio.subarray(0, 960)));
101102 // Split the modes into three slices and sum the parts.
102- const sigma = modalRates(model, pl.loss);
103+ const sigma = modalRates(model);
103104 const g = strikeGains(model, pl, p);
104105 const c = simulateContact(model, sigma, g, p, p.fps);
105106 const tStart = tailStart(c, maxDelay);
@@ -130,7 +131,7 @@ test("a lighter hit gives a proportionally smaller signal in the linear model",
130131 const { model, pl } = await smallModel();
131132 const p1: SimParams = { ...defaultParams, fs: 48000, duration: 0.05, strikeVelocity: 0.2 };
132133 const p2: SimParams = { ...p1, strikeVelocity: 0.1 };
133- const a1 = synthesize(model, pl, pl.loss, p1, false).audio, a2 = synthesize(model, pl, pl.loss, p2, false).audio;
134+ const a1 = synthesize(model, pl, p1, false).audio, a2 = synthesize(model, pl, p2, false).audio;
134135 const peak = (a: Float32Array) => a.reduce((m, v) => Math.max(m, Math.abs(v)), 0);
135136 // The Hertz contact is nonlinear (a faster hit is shorter and brighter), so the ratio exceeds 2 somewhat.
136137 const ratio = peak(a1) / peak(a2);
tests/tutorial.test.tsmodified+6−3View file
@@ -12,6 +12,7 @@ import { defaultParams, type SimParams } from "../src/params.ts";
1212 import { assemble } from "../src/model/assemble.ts";
1313 import { computeModes, localRunner } from "../src/modal/modes.ts";
1414 import { modalModel, synthesize } from "../src/modal/synth.ts";
15+import { computeDamping, dampingPlate } from "../src/modal/damping.ts";
1516 import { STEPS, chapters, emptyContext, rectangleScript, stepIndex, type TutorialContext, type TutorialHost } from "../src/tutorial/steps.ts";
1617
1718 function check(html: string, where: string): void {
@@ -45,13 +46,15 @@ class TestHost implements TutorialHost {
4546 if (!pl || this.c.modes) return;
4647 const sys = assemble(pl);
4748 const modes = await computeModes(sys, this.params.modeMaxFreq, localRunner(sys));
48- this.c = { ...this.c, sys, modes, model: modalModel(pl, sys, modes) };
49+ const model = modalModel(pl, sys, modes);
50+ computeDamping(model, dampingPlate(pl));
51+ this.c = { ...this.c, sys, modes, model };
4952 }
5053 showMode(m: number | null): void { this.shown.push(m); }
5154 async synthesize(): Promise<void> {
5255 await this.computeModes();
5356 const pl = this.c.plate!, model = this.c.model!;
54- const r = synthesize(model, pl, pl.loss, { ...this.params, duration: 0.05 }, false);
57+ const r = synthesize(model, pl, { ...this.params, duration: 0.05 }, false);
5558 this.c = { ...this.c, lastContact: r.contact, recording: { source: "linear", seconds: 0.05 } };
5659 }
5760 play(): void { this.played++; }
@@ -129,7 +132,7 @@ test("the scaling actions remember mode lists that compare as predicted", async
129132 });
130133
131134 test("the rectangle script compiles with each option", () => {
132- for (const o of [{}, { thicknessMm: 1.58 }, { widthIn: 6, heightIn: 9 }, { t60: "1" }, { t60: "[100 5; 1000 2]" }]) {
135+ for (const o of [{}, { thicknessMm: 1.58 }, { widthIn: 6, heightIn: 9 }, { eta: 1e-3 }, { thermal: false }]) {
133136 const d = compileDesign(rectangleScript(o));
134137 assert.ok(d.Lx > 0);
135138 }