Model the felt as a Maxwell chain so that the full model has the loss factor of the linear one
The full model damped the pads with a dashpot matched to the bounce frequency, which over-damps a mode at omega by omega / omega_b; on the aluminum square on the default stand the kilohertz modes decayed in hundredths of a second instead of seconds. The felt is now a spring in parallel with a chain of twelve Maxwell arms, two per decade from 0.5 Hz, with geometric weights that give a loss factor flat within 6 percent across the band, normalized so that the storage stiffness is one at the bounce. The full model steps the arm extensions exactly per step in a new storage buffer, with the relaxed part compression-only and the arm force scaled by the pads in contact; the linear model uses the chain's loss stiffness at each mode's frequency; the static preload uses the relaxed stiffness.
Tests: chain flatness and the exact arm step, work per cycle of the stepped chain on a plate, the self-test now checks the felt state too, and a browser test compares what the stand adds to the decay of a light hit on the aluminum square between the models with the stand on and off (within 3 dB in four bands). A dense eigensolution of the time-domain system confirms the full model and shows the linear model's first-order felt rates over-estimate the damping of felt-heavy modes by 13 to 66 percent; documented in THEORY.md.
18 changed files+497−90
README.mdmodified+3−2View file
@@ -49,7 +49,7 @@ The decay of the sound is not prescribed but derived: from the loss factor of th
4949
5050 **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.
5151
52-**Stand.** The Stand group puts the plate on the stand or takes it off (a free plate), and sets the felt annulus (inner and outer radius, thickness, Young's modulus in compression, loss factor), whether there is a top pad under the wing nut, its preload, and the position of the rod relative to the center of the bounding box. The pads act on the plate nodes within the annulus whether or not the outline has a hole there; the Notes under the group give the total stiffness, the bounce frequency of the rigid plate, the number of pad nodes and the preloads. The stand is on by default, so a preset without a hole, such as the flat rectangle, sits on it until the box is unticked.
52+**Stand.** The Stand group puts the plate on the stand or takes it off (a free plate), and sets the felt annulus (inner and outer radius, thickness, Young's modulus in compression, loss factor), whether there is a top pad under the wing nut, its preload, and the position of the rod relative to the center of the bounding box. The pads act on the plate nodes within the annulus whether or not the outline has a hole there; the Notes under the group give the total stiffness, the bounce frequency of the rigid plate, the number of pad nodes, the preloads, and the fraction of the stiffness the felt keeps at rest. The felt is viscoelastic, a chain of relaxation times that gives it the same loss factor at every frequency, so the modulus is the value at the bounce frequency and rises slowly above it; both models use the same chain (THEORY.md, Section 3). The stand is on by default, so a preset without a hole, such as the flat rectangle, sits on it until the box is unticked.
5353
5454 **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.
5555
@@ -66,11 +66,12 @@ The decay of the sound is not prescribed but derived: from the loss factor of th
6666 `npm test` runs eight groups of tests with Node's built-in runner; `npm run test:cpu` runs the first seven without a browser.
6767
6868 - `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 from the parameters with its rates and preloads, the loss constants of the material), and that every preset evaluates on the default grid.
69-- `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).
69+- `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, the work per cycle of the felt's chain, thickness scaling, the stand's stiffness against translation and tilt, the lift-off rule).
7070 - `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.
7171 - `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.
7272 - `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.
7373 - `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.
74+- `tests/felt.test.ts`: the felt's chain (flat loss factor across the band, storage of one at the reference, the exact step of an arm against a ramp and a sinusoid).
7475 - `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.
7576 - `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.
7677
THEORY.mdmodified+12−6View file
@@ -73,9 +73,15 @@ $$
7373
7474 where the sums run over the plate nodes within the annulus, $k_k = n_p (E_p/t)\, \Delta A_k$ and $g_k = n_p (G_p/t)\, \Delta A_k$, $n_p$ is 1 or 2 for one or two pads, and $\Delta A_k$ is the node's share of the area, scaled so that the shares add up to the exact area of the annulus (a coarse grid covers the annulus with a dozen nodes). A uniform translation of the plate then feels the total stiffness $K = n_p (E_p/t) \pi (r_2^2 - r_1^2)$, and a tilt feels a restoring moment set by the second moment of the annulus, so the finite radius of the felt is what holds the cymbal level. With the defaults of the presets (a felt of 12 mm at 0.5 MPa between 7 and 20 mm, two pads) a 16 inch bronze cymbal of about 1 kg has a total stiffness of about 90 kN/m, and treated as a rigid body it would bounce on the stand at about 45 Hz and rock at about 5 Hz. The rock does come out near that value. The bounce comes out well below it, near 15 Hz for the flat disc preset and 19 Hz for the thicker ride, because a thin plate is not rigid: the pads push on its center, whose compliance adds, loosely speaking, in series with that of the felt, so the bounce mode is a rigid translation blended with the umbrella-shaped plate mode (one nodal circle, near 38 Hz for the uniform disc), and the higher plate modes that move the center take a share of the felt's damping. The felt modulus is an estimate; felt is a nonlinear material whose modulus rises with compression, and the value can be changed in the Stand parameters.
7575
76-**Preload and lift-off.** The wing nut compresses the top pad by a preload $\delta_t$ at rest (1 mm by default), and the bottom pad carries that preload plus the weight of the plate: $\delta_b = \delta_t + Mg/K_1$ with $K_1$ the stiffness of one pad. Without a top pad $\delta_t = 0$ and the bottom pad carries the weight alone. Gravity enters only through this preload; the static sag of the plate under its own weight, a fraction of a millimeter, is neglected, and $w = 0$ is the rest state. In the full model each pad pushes only while compressed: the bottom pad by $\delta_b - w$, the top pad by $\delta_t + w$, and the constant part of the two forces, which balances the weight, is subtracted, so that within the preload the pads are the linear spring above and beyond it the cymbal lifts off one pad and loads the other. A hard hit near the edge lifts the cymbal off the bottom pad and starts a rock of a centimeter or two at the edge, as on a real stand. The linear model uses the linear spring throughout; it cannot represent the lift-off, which is a nonlinear effect of the stand rather than of the plate.
76+**Preload and lift-off.** The wing nut compresses the top pad by a preload $\delta_t$ at rest (1 mm by default), and the bottom pad carries that preload plus the weight of the plate: $\delta_b = \delta_t + Mg/(k_\infty K_1)$ with $K_1$ the stiffness of one pad and $k_\infty$ the relaxed fraction of the felt's stiffness defined below. Without a top pad $\delta_t = 0$ and the bottom pad carries the weight alone. Gravity enters only through this preload; the static sag of the plate under its own weight, a fraction of a millimeter, is neglected, and $w = 0$ is the rest state. In the full model each pad pushes only while compressed: the bottom pad by $\delta_b - w$, the top pad by $\delta_t + w$, and the constant part of the two forces, which balances the weight, is subtracted, so that within the preload the pads are the linear spring above and beyond it the cymbal lifts off one pad and loads the other. A hard hit near the edge lifts the cymbal off the bottom pad and starts a rock of a centimeter or two at the edge, as on a real stand. The linear model uses the linear spring throughout; it cannot represent the lift-off, which is a nonlinear effect of the stand rather than of the plate.
7777
78-**Damping of the felt.** A loss factor $\eta$ means a stiffness with a phase lag, so that a mode of angular frequency $\omega$ whose shape puts the energy $\tfrac12 \phi^T K_s \phi$ into the pads decays at the additional rate $\sigma = \eta\, \phi^T K_s \phi / (2\omega)$; the linear model applies exactly this to each mode. The full model, which has no modes, replaces the loss factor by a dashpot with the rate $c_k = \eta k_k / \omega_b$ at every pad node, where $\omega_b$ is the bounce frequency; the two agree for the bounce and the rock and differ, in the sense of a stronger dashpot loss on the high modes, above it. The high modes of the plate put next to no energy into the pads, so the difference is small.
78+**Damping of the felt.** A loss factor $\eta$ means a stiffness with a phase lag, $K^*(\omega) = K'(\omega)(1 + i\eta)$, so that a mode of angular frequency $\omega$ whose shape puts the energy $\tfrac12 \phi^T K_s \phi$ into the pads decays at the additional rate $\sigma = \eta\, \phi^T K_s \phi / (2\omega)$. A loss factor that is the same at every frequency cannot, however, come from a spring and a dashpot: a dashpot of rate $c$ gives $\eta = c\omega/K$, which grows with frequency. An earlier version of the full model used a dashpot matched to the bounce frequency $\omega_b$, and it over-damped the modes at $\omega$ by the factor $\omega/\omega_b$, which for a light plate on the default stand, whose bounce is near 100 Hz and many of whose modes move the center, shortened the decay of the kilohertz modes from seconds to a few hundredths of a second. We therefore model the felt as a viscoelastic solid, a spring $k_\infty$ in parallel with a chain of Maxwell arms (a spring $k_i$ in series with a dashpot, relaxation time $\tau_i$), whose complex stiffness per unit of the pad stiffness is
79+
80+$$
81+\hat K(\omega) = k_\infty + \sum_i k_i \frac{i\omega\tau_i}{1 + i\omega\tau_i}.
82+$$
83+
84+With the arm frequencies $1/(2\pi\tau_i)$ spaced two per decade from 0.5 Hz to 160 kHz and the stiffnesses growing geometrically, $k_{i+1}/k_i = r^{2\delta/\pi}$ with $r = \sqrt{10}$ the frequency ratio and $\delta = \arctan\eta$ the loss angle, the chain approximates the power law $\hat K \propto (i\omega)^{2\delta/\pi}$ of a material with a constant loss angle, and its loss factor $\mathrm{Im}\,\hat K/\mathrm{Re}\,\hat K$ stays within 6 percent of $\eta$ from 5 Hz to 10 kHz for $\eta \le 0.3$ (within 10 percent for $\eta = 0.6$). The chain is normalized so that its storage stiffness $\mathrm{Re}\,\hat K$ is 1 at the bounce frequency, where the felt's stiffness matters most: the modulus $E_p$ of the parameters is the modulus at the bounce, and the felt is softer below it and stiffer above, by the factor $10^{2\delta/\pi}$ per decade, about 1.5 for $\eta = 0.3$, as the Kramers-Kronig relations require of any material with a loss. At rest it holds the fraction $k_\infty$ of that stiffness, about 0.4 for $\eta = 0.3$, which is the stiffness the static preload and the weight see. The linear model uses the storage stiffness at the bounce in $K_s$ and gives each mode the rate $\sigma = \mathrm{Im}\,\hat K(\omega)\, \phi^T K_s \phi / (2\omega)$, which is the hysteretic rate above with the felt's actual loss stiffness at the mode's frequency. The full model carries the chain in time: each arm has at every pad node and for each of $u$, $v$, $w$ one state, its extension $p_i$, which obeys $\dot p_i = \dot w - p_i/\tau_i$, and the felt's force on the node is $-k_k[k_\infty w + \sum_i k_i p_i]$, with the compression-only rule of the previous paragraph applied to the relaxed part and the arm force scaled by the number of pads in contact. The update over a step in which the node moves by $\Delta w$ is exact for a linear motion within the step, $p_i \leftarrow p_i e^{-\Delta t/\tau_i} + \Delta w\,(\tau_i/\Delta t)(1 - e^{-\Delta t/\tau_i})$, so that arms far faster than the step are handled as well as slow ones; a node that has left every pad relaxes its transverse arms at once. The two solvers now describe the same felt, with one difference that remains. The rate of the linear model is first order in the felt's coupling, and a strongly damped element at one place also reshapes the modes that load it; for a mode with a large share of its stiffness in the felt the first-order rate is an over-estimate. We checked this on the aluminum square preset (12 mm grid) by a dense eigensolution of the linearized time-domain system, plate, pads and arms together: the mode at 117 Hz, which holds about a sixth of its stiffness in the felt, decays in 0.41 s against the 0.37 s of the first-order rate, one at 359 Hz in 0.21 s against 0.15 s, and one at 745 Hz in 0.59 s against 0.36 s, while modes with a small share agree within a few percent; the rock, whose stiffness is all felt, lies at 9.6 Hz against 11.3 Hz in the linear model, because the felt is softer at the rock's frequency than at the bounce where the modulus is defined. The full model, which steps the chain, reproduces the exact values. The test of Section 6 confirms that what the stand adds to the decay of a light hit agrees between the two.
7985
8086 ## 4. The linear model
8187
@@ -118,7 +124,7 @@ which is exact for the free decay and includes the directivity of every mode. Th
118124
119125 ## 5. The full model
120126
121-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.
127+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 are evaluated explicitly at the old velocity, and the felt's arms at their state from the previous step, which is then advanced with the new displacement. 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.
122128
123129 **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.
124130
@@ -130,9 +136,9 @@ The full model keeps the nonlinear strains of Section 1 and steps the semi-discr
130136
131137 ## 6. Validation
132138
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.
139+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 only remove energy; that the felt's chain has the loss factor it was built for within 6 percent across the band and a storage stiffness of 1 at the reference, that the exact step of an arm reproduces its response to a ramp for steps from a fiftieth to twice the relaxation time, and that the chain stepped in time on a plate on a stand dissipates the work $\pi\,\mathrm{Im}\,\hat K\, K A^2$ per cycle of a rigid translation of amplitude $A$; 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.
134140
135-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.
141+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 and the felt's arms at random extensions, 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, and the updated arm extensions agree to better than $10^{-6}$. 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 and the modal truncation at 16 kHz. Finally, a light hit on the aluminum square, the case that exposed the dashpot, is simulated both ways with the stand on and off, and what the stand adds to the decay of each band between 0.25 and 0.55 s (the drop with the stand on minus the drop without it, which cancels the fitted plate loss against the per-mode one) must agree between the models within 4 dB in the four bands from 50 Hz to 6 kHz; the differences are 1.4, 0.4, 0.7 and 2.8 dB. With the dashpot the same measure differed by tens of decibels above 600 Hz.
136142
137143 ## 7. Limitations
138144
@@ -140,4 +146,4 @@ The linear model is linear: it has no pitch glide, no wash, no lift-off, and a h
140146
141147 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.
142148
143-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.
149+The stand is a linear viscoelastic felt with a loss factor that is the same at every frequency; real felt stiffens with compression, its loss factor is not exactly flat, 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/loss.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/felt.test.ts tests/tutorial.test.ts",
1313 "test:gpu": "node --test tests/gpu.test.ts"
1414 },
1515 "devDependencies": {
src/felt.tsadded+85−0View file
@@ -0,0 +1,85 @@
1+// The felt of the stand as a viscoelastic material: a spring in parallel
2+// with a chain of Maxwell arms (a spring and a dashpot in series), whose
3+// relaxation times are spaced two per decade from 0.5 Hz to above the audio
4+// band and whose weights grow geometrically. Such a chain has a loss factor
5+// that is nearly flat over the band, which is what a single loss factor
6+// eta means, and unlike a dashpot it is causal and can be stepped in time:
7+// each arm carries one state per node, its extension p, which relaxes
8+// toward zero at its own rate while the node moves.
9+//
10+// Every quantity here is normalized to the stiffness at the reference
11+// frequency (the bounce of the plate on the stand), where the storage
12+// stiffness of the chain is 1, so that the felt modulus given by the user is
13+// the modulus at the frequency where it matters most. Below the reference
14+// the felt is softer and above it stiffer, as the Kramers-Kronig relations
15+// require of any material with a loss; at rest its stiffness is kInf.
16+
17+export interface FeltChain {
18+ /** Relaxation time of each arm [s]. */
19+ tau: Float64Array;
20+ /** Stiffness of each arm, per unit of the stiffness at the reference frequency. */
21+ k: Float64Array;
22+ /** Relaxed (static) stiffness, per unit of the stiffness at the reference frequency. */
23+ kInf: number;
24+ /** The loss factor the chain was built for. */
25+ eta: number;
26+}
27+
28+/** Lowest arm frequency [Hz]; the highest is FELT_ARMS_PER_DECADE arms per decade above it. */
29+export const FELT_LOWEST_HZ = 0.5;
30+export const FELT_ARMS_PER_DECADE = 2;
31+/** Number of arms: 0.5 Hz to 158 kHz. */
32+export const FELT_ARMS = 12;
33+
34+/**
35+ * The chain for a loss factor eta, normalized at the frequency fRef [Hz]. A
36+ * chain of arms with frequencies in the ratio r and stiffnesses in the ratio
37+ * 1 + c approximates the power law K(i omega) ~ (i omega)^(2 delta / pi) of a
38+ * material with the constant loss angle delta = atan(eta), whose modulus
39+ * rises by the factor r^(2 delta / pi) per arm; c is that factor minus one.
40+ */
41+export function feltChain(eta: number, fRef: number): FeltChain {
42+ const r = Math.pow(10, 1 / FELT_ARMS_PER_DECADE);
43+ const c = eta > 0 ? Math.pow(r, (2 * Math.atan(eta)) / Math.PI) - 1 : 0;
44+ const nArms = c > 0 ? FELT_ARMS : 0;
45+ const tau = new Float64Array(nArms), k = new Float64Array(nArms);
46+ let S = 1;
47+ for (let i = 0; i < nArms; i++) {
48+ tau[i] = 1 / (2 * Math.PI * FELT_LOWEST_HZ * Math.pow(r, i));
49+ k[i] = c * S;
50+ S *= 1 + c;
51+ }
52+ const raw: FeltChain = { tau, k, kInf: 1, eta };
53+ const scale = 1 / feltModulus(raw, 2 * Math.PI * Math.max(fRef, 1e-6)).storage;
54+ for (let i = 0; i < nArms; i++) k[i] *= scale;
55+ return { tau, k, kInf: scale, eta };
56+}
57+
58+/** Storage and loss stiffness of the chain at the angular frequency omega, per unit of the reference stiffness. */
59+export function feltModulus(ch: FeltChain, omega: number): { storage: number; loss: number } {
60+ let storage = ch.kInf, loss = 0;
61+ for (let i = 0; i < ch.tau.length; i++) {
62+ const x = omega * ch.tau[i], d = 1 + x * x;
63+ storage += (ch.k[i] * x * x) / d;
64+ loss += (ch.k[i] * x) / d;
65+ }
66+ return { storage, loss };
67+}
68+
69+/**
70+ * The constants of the exact time step of the arms for a step dt: with the
71+ * node moving linearly by dw over the step, the extension of arm i goes from
72+ * p to p decay_i + dw growth_i, where decay_i = exp(-dt / tau_i) and
73+ * growth_i = (tau_i / dt)(1 - decay_i). The force of the arms on the node is
74+ * -k_node sum_i k_i p_i.
75+ */
76+export function feltStepConstants(ch: FeltChain, dt: number): { decay: Float64Array; growth: Float64Array } {
77+ const n = ch.tau.length;
78+ const decay = new Float64Array(n), growth = new Float64Array(n);
79+ for (let i = 0; i < n; i++) {
80+ const x = dt / ch.tau[i];
81+ decay[i] = Math.exp(-x);
82+ growth[i] = -Math.expm1(-x) / x;
83+ }
84+ return { decay, growth };
85+}
src/main.tsmodified+3−3View file
@@ -543,7 +543,7 @@ class App {
543543 } else lines.push("Flat plate (no profile): the linear model is bending only");
544544 if (pl.stand && sr) {
545545 const pad = pl.stand.pad;
546- this.form.setNotes("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 in all, bounce of the rigid plate at ${sr.bounceHz.toFixed(0)} Hz, ${sr.nodes.length} pad nodes, preload ${(sr.preloadTop * 1e3).toFixed(2)} mm on top and ${(sr.preloadBottom * 1e3).toFixed(2)} mm below. The pads act on the plate nodes within the annulus whether or not the outline has a hole there.`);
546+ this.form.setNotes("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 in all, bounce of the rigid plate at ${sr.bounceHz.toFixed(0)} Hz, ${sr.nodes.length} pad nodes, preload ${(sr.preloadTop * 1e3).toFixed(2)} mm on top and ${(sr.preloadBottom * 1e3).toFixed(2)} mm below. The felt is viscoelastic: its modulus is the value at the bounce, and at rest it keeps ${(100 * sr.chain.kInf).toFixed(0)}% of it. The pads act on the plate nodes within the annulus whether or not the outline has a hole there.`);
547547 lines.push(`On the stand (see the Stand parameters): bounce ${sr.bounceHz.toFixed(0)} Hz`);
548548 } else {
549549 this.form.setNotes("Stand", pl.stand ? "The felt annulus covers no plate node: move the rod or enlarge the pads." : "Free plate: the stand is off.");
@@ -744,7 +744,7 @@ class App {
744744 const p = this.form.get();
745745 this.sim?.destroy();
746746 const fit = fitSigma(pl);
747- const loss: FullLoss = fullLoss(fit, standRates(pl), inPlaneRate(pl));
747+ const loss: FullLoss = fullLoss(fit, inPlaneRate(pl));
748748 this.sim = new Simulation(this.device, pl, p, loss);
749749 this.rec = { audio: new Float32Array(0), samples: 0, fs: p.fs, frames: [new Float32Array(pl.grid.n)], fps: p.fps, source: "full" };
750750 this.resetLoudness();
@@ -1036,7 +1036,7 @@ class App {
10361036 const d = compileDesign(pr.source);
10371037 const pl = evaluatePlate(d, 0.012, standFromParams(this.form.get()));
10381038 const r = await runSelfTest(this.device, pl, this.form.get());
1039- lines.push(`${id}: max relative error ${r.maxRel.toExponential(2)} (${r.n} nodes)`);
1039+ lines.push(`${id}: max relative error ${r.maxRel.toExponential(2)}, felt state error ${r.padRel.toExponential(2)} (${r.n} nodes)`);
10401040 }
10411041 return lines.join("\n");
10421042 }
src/modal/synth.tsmodified+7−4View file
@@ -16,6 +16,7 @@
1616 import { type LinearSystem } from "../model/assemble.ts";
1717 import { type ModeSet } from "./modes.ts";
1818 import { type Plate, type SimParams, STICK_RELEASE, listenerTerms, standRates, strikePatch, nodeMass } from "../params.ts";
19+import { feltModulus } from "../felt.ts";
1920
2021 /** What the synthesis needs of a plate and its modes; plain arrays so that it can be posted to a worker. */
2122 export interface ModalModel {
@@ -58,9 +59,10 @@ export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalM
5859 // Hysteretic damping: a loss factor eta on a store of energy U gives a mode
5960 // the rate sigma = eta U / omega. The pads hold the energy phi' K_pad phi / 2
6061 // 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.
62+ // mass-normalized mode; the metal has its loss factor, and the felt the
63+ // loss stiffness of its chain at the mode's frequency, which is eta times
64+ // the storage stiffness there. The thermoelastic and the radiation rates
65+ // need the shapes on the grid and are filled in by damping.ts.
6466 const padSigma = new Float64Array(M), sigmaMat = new Float64Array(M);
6567 const sr = standRates(pl);
6668 for (let m = 0; m < M; m++) {
@@ -79,7 +81,8 @@ export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalM
7981 }
8082 const omega = 2 * Math.PI * modes.freqs[m];
8183 if (!(omega > 0)) continue;
82- padSigma[m] = sr ? (sr.eta * e) / (2 * omega) : 0;
84+ // The felt's chain has the loss stiffness L(omega) per unit of its stiffness at the bounce (felt.ts).
85+ padSigma[m] = sr ? (feltModulus(sr.chain, omega).loss * e) / (2 * omega) : 0;
8386 sigmaMat[m] = Math.max(0, pl.mat.eta * (omega / 2 - e / (2 * omega)));
8487 }
8588 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 };
src/model/reference.tsmodified+89−26View file
@@ -23,6 +23,7 @@ import {
2323 standRates,
2424 } from "../params.ts";
2525 import { type SigmaFit } from "../loss.ts";
26+import { feltStepConstants } from "../felt.ts";
2627
2728 export interface Fields {
2829 u: Float64Array;
@@ -58,22 +59,90 @@ export function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float
5859 }
5960
6061 /**
61- * Force of the pads at a node with transverse displacement w. Without
62- * lift-off the pads are a linear spring of rate kw. With lift-off each pad
63- * pushes only while compressed: the bottom pad by (preloadBottom - w), the
64- * top pad by (preloadTop + w); the constant part of the two, which balances
65- * the weight of the plate, is subtracted so that w = 0 is the rest state.
62+ * Relaxed (static) force of the pads at a node with transverse displacement
63+ * w: the felt's relaxed stiffness kInf times the stiffness at the bounce
64+ * frequency, the rest of the felt's force being carried by the arms of its
65+ * chain (addPadMemory). Without lift-off the pads are a linear spring of
66+ * rate kInf kw. With lift-off each pad pushes only while compressed: the
67+ * bottom pad by (preloadBottom - w), the top pad by (preloadTop + w); the
68+ * constant part of the two, which balances the weight of the plate, is
69+ * subtracted so that w = 0 is the rest state.
6670 */
6771 export function padForce(sr: StandRates, k: number, w: number, liftoff: boolean): number {
68- if (!liftoff) return -sr.kw[k] * w;
69- const kb = sr.kwBottom[k];
70- const hasTop = sr.kw[k] > kb + 1e-300;
72+ const s = sr.chain.kInf;
73+ if (!liftoff) return -s * sr.kw[k] * w;
74+ const kb = s * sr.kwBottom[k];
75+ const hasTop = sr.kw[k] > sr.kwBottom[k] + 1e-300;
7176 const bottom = kb * Math.max(sr.preloadBottom - w, 0);
7277 const top = hasTop ? kb * Math.max(sr.preloadTop + w, 0) : 0;
7378 const rest = kb * sr.preloadBottom - (hasTop ? kb * sr.preloadTop : 0);
7479 return bottom - top - rest;
7580 }
7681
82+/** Number of pads in contact with the node at transverse displacement w (both pads without lift-off). */
83+export function padsInContact(sr: StandRates, k: number, w: number, liftoff: boolean): number {
84+ const hasTop = sr.kw[k] > sr.kwBottom[k] + 1e-300;
85+ if (!liftoff) return hasTop ? 2 : 1;
86+ return (sr.preloadBottom - w > 0 ? 1 : 0) + (hasTop && sr.preloadTop + w > 0 ? 1 : 0);
87+}
88+
89+/**
90+ * The extensions of the arms of the felt's chain at every node, for the
91+ * three fields: pad[(field * nArms + arm) * n + k] with field 0, 1, 2 for
92+ * u, v, w. Zero is the relaxed felt, which is the initial state; the layout
93+ * is that of the GPU buffer.
94+ */
95+export type PadState = Float64Array;
96+export function padState(pl: Plate, sr: StandRates | null): PadState {
97+ return new Float64Array(3 * (sr ? sr.chain.tau.length : 0) * pl.grid.n);
98+}
99+
100+/**
101+ * Adds to the accelerations `out` the force of the arms of the felt's
102+ * chain, -k_node sum_i k_i p_i on each field, with the transverse force
103+ * scaled by the number of pads in contact (the shear of the felt is applied
104+ * whether or not a pad has lifted, as its static part is).
105+ */
106+export function addPadMemory(pl: Plate, sr: StandRates | null, f: Fields, pad: PadState, liftoff: boolean, out: Fields): void {
107+ if (!sr) return;
108+ const g = pl.grid, n = g.n, ch = sr.chain, nArms = ch.tau.length;
109+ for (const k of sr.nodes) {
110+ const mass = nodeMass(pl, k % g.nx, Math.floor(k / g.nx));
111+ if (mass === 0) continue;
112+ let su = 0, sv = 0, sw = 0;
113+ for (let a = 0; a < nArms; a++) {
114+ su += ch.k[a] * pad[a * n + k];
115+ sv += ch.k[a] * pad[(nArms + a) * n + k];
116+ sw += ch.k[a] * pad[(2 * nArms + a) * n + k];
117+ }
118+ const nc = padsInContact(sr, k, f.w[k], liftoff);
119+ out.u[k] -= (sr.ku[k] * su) / mass;
120+ out.v[k] -= (sr.ku[k] * sv) / mass;
121+ out.w[k] -= (nc * sr.kwBottom[k] * sw) / mass;
122+ }
123+}
124+
125+/**
126+ * Advances the arm extensions over one step of dt in which the fields moved
127+ * by `dw` (the increments of u, v, w) from the displacements `f`: exact for
128+ * a linear motion within the step. At a node that had left every pad at the
129+ * start of the step the transverse arms relax at once.
130+ */
131+export function stepPadMemory(pl: Plate, sr: StandRates | null, dt: number, dw: Fields, f: Fields, pad: PadState, liftoff: boolean): void {
132+ if (!sr) return;
133+ const n = pl.grid.n, ch = sr.chain, nArms = ch.tau.length;
134+ const { decay, growth } = feltStepConstants(ch, dt);
135+ for (const k of sr.nodes) {
136+ const lifted = padsInContact(sr, k, f.w[k], liftoff) === 0;
137+ for (let a = 0; a < nArms; a++) {
138+ pad[a * n + k] = pad[a * n + k] * decay[a] + dw.u[k] * growth[a];
139+ pad[(nArms + a) * n + k] = pad[(nArms + a) * n + k] * decay[a] + dw.v[k] * growth[a];
140+ const iw = (2 * nArms + a) * n + k;
141+ pad[iw] = lifted ? 0 : pad[iw] * decay[a] + dw.w[k] * growth[a];
142+ }
143+ }
144+}
145+
77146 /**
78147 * Undamped accelerations (out.u, out.v, out.w) of every node for the
79148 * displacement fields f, measured from the stress-free shape pl.w0. Bending
@@ -162,7 +231,7 @@ export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false
162231 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
163232 F[0] = 0; F[1] = 0; F[2] = -hx * hy * bend;
164233 if (sr) {
165- F[0] -= sr.ku[k] * u[k]; F[1] -= sr.ku[k] * v[k];
234+ F[0] -= sr.chain.kInf * sr.ku[k] * u[k]; F[1] -= sr.chain.kInf * sr.ku[k] * v[k];
166235 F[2] += padForce(sr, k, w[k], liftoff);
167236 }
168237 cellForce(i, j, 0, F);
@@ -176,28 +245,24 @@ export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false
176245 }
177246 }
178247
179-/** The loss terms of the full simulation: the fitted sigma model on the transverse motion, an in-plane rate, and the pads' dashpots. */
248+/** The loss terms of the full simulation on the plate itself: the fitted sigma model on the transverse motion and an in-plane rate. The felt's loss is in its chain (addPadMemory). */
180249 export interface FullLoss {
181250 sigma0: number; // [1/s], applied as 2 sigma0 on all components
182251 sigma1: number; // [m^2/s], applied as 2 sigma1 times the Laplacian of the transverse velocity
183252 sigmaIn: number; // [1/s], extra rate on the in-plane motion
184- /** Dashpot rate of the pads per unit of pad stiffness [s]: c_k = padDamping * kw_k. */
185- padDamping: number;
186253 }
187254
188-/** The loss terms of the full simulation from the fitted sigma model and the stand. */
189-export function fullLoss(fit: SigmaFit, sr: StandRates | null, inPlaneRate: number): FullLoss {
190- const padDamping = sr && sr.bounceHz > 0 ? sr.eta / (2 * Math.PI * sr.bounceHz) : 0;
191- return { sigma0: fit.sigma0, sigma1: fit.sigma1, sigmaIn: inPlaneRate, padDamping };
255+/** The loss terms of the full simulation from the fitted sigma model. */
256+export function fullLoss(fit: SigmaFit, inPlaneRate: number): FullLoss {
257+ return { sigma0: fit.sigma0, sigma1: fit.sigma1, sigmaIn: inPlaneRate };
192258 }
193259
194260 /**
195261 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
196262 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
197- * transverse component, -2 sigmaIn v on the in-plane components, and the
198- * dashpots of the pads over their annulus.
263+ * transverse component, and -2 sigmaIn v on the in-plane components.
199264 */
200-export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields, sr: StandRates | null = standRates(pl)): void {
265+export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields): void {
201266 const g = pl.grid;
202267 const { sigma0, sigma1, sigmaIn } = loss;
203268 const { nx, ny, hx, hy } = g;
@@ -213,16 +278,14 @@ export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields,
213278 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
214279 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
215280 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
216- const cw = sr ? (loss.padDamping * sr.kw[k]) / mass : 0;
217- const cu = sr ? (loss.padDamping * sr.ku[k]) / mass : 0;
218- out.u[k] -= (2 * (sigma0 + sigmaIn) + cu) * vel.u[k];
219- out.v[k] -= (2 * (sigma0 + sigmaIn) + cu) * vel.v[k];
220- out.w[k] += -(2 * sigma0 + cw) * c + 2 * sigma1 * lap;
281+ out.u[k] -= 2 * (sigma0 + sigmaIn) * vel.u[k];
282+ out.v[k] -= 2 * (sigma0 + sigmaIn) * vel.v[k];
283+ out.w[k] += -2 * sigma0 * c + 2 * sigma1 * lap;
221284 }
222285 }
223286 }
224287
225-/** Discrete potential energy (bending + membrane + pads, linear pads) whose negative gradient is the force in `accelerations`. */
288+/** Discrete potential energy (bending + membrane + the relaxed pads, linear) whose negative gradient is the force in `accelerations`. */
226289 export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = standRates(pl)): number {
227290 const g = pl.grid, m = pl.mat, w0 = pl.w0;
228291 const { nx, ny, hx, hy } = g;
@@ -275,7 +338,7 @@ export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = st
275338 }
276339
277340 let Us = 0;
278- if (sr) for (const k of sr.nodes) Us += 0.5 * (sr.ku[k] * (u[k] * u[k] + v[k] * v[k]) + sr.kw[k] * w[k] * w[k]);
341+ if (sr) for (const k of sr.nodes) Us += 0.5 * sr.chain.kInf * (sr.ku[k] * (u[k] * u[k] + v[k] * v[k]) + sr.kw[k] * w[k] * w[k]);
279342 return Ub + Um + Us;
280343 }
281344
src/params.tsmodified+15−6View file
@@ -3,6 +3,8 @@
33 // 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
6+import { type FeltChain, feltChain } from "./felt.ts";
7+
68 export interface SimParams {
79 // Strike
810 strikeFracX: number; // strike point as a fraction of the bounding-box width from the left edge
@@ -356,8 +358,11 @@ export function thicknessRange(pl: Plate): { min: number; max: number } {
356358 * The pads as per-node spring rates. `kw[k]` is the transverse stiffness of
357359 * the pads at node k (both pads summed), `ku[k]` the in-plane (shear)
358360 * stiffness, both in N/m; `kwBottom` is the bottom pad alone, which the full
359- * simulation needs for the lift-off rule. `eta` is the loss factor, and
360- * `preloadTop`, `preloadBottom` are the pad compressions at rest.
361+ * simulation needs for the lift-off rule. All three are the stiffness at the
362+ * bounce frequency; `chain` is the viscoelastic model of the felt (felt.ts),
363+ * normalized to that stiffness, whose `kInf` is the relaxed fraction that
364+ * holds a static load. `eta` is the loss factor, and `preloadTop`,
365+ * `preloadBottom` are the pad compressions at rest.
361366 */
362367 export interface StandRates {
363368 kw: Float64Array;
@@ -365,6 +370,7 @@ export interface StandRates {
365370 kwBottom: Float64Array;
366371 nodes: number[];
367372 eta: number;
373+ chain: FeltChain;
368374 preloadTop: number;
369375 preloadBottom: number;
370376 /** Total transverse stiffness of the pads [N/m]. */
@@ -401,7 +407,7 @@ export function standRates(pl: Plate): StandRates | null {
401407 nodeArea += g.hx * g.hy * wgt;
402408 }
403409 }
404- const empty = { kw, ku, kwBottom, nodes, eta: pad.lossFactor, preloadTop: 0, preloadBottom: 0, totalKw: 0, bounceHz: 0 };
410+ const empty = { kw, ku, kwBottom, nodes, eta: pad.lossFactor, chain: feltChain(0, 1), preloadTop: 0, preloadBottom: 0, totalKw: 0, bounceHz: 0 };
405411 if (nodes.length === 0 || nodeArea === 0) return empty;
406412 const scale = annulus / nodeArea;
407413 let totalKw = 0;
@@ -413,10 +419,13 @@ export function standRates(pl: Plate): StandRates | null {
413419 totalKw += kw[k];
414420 }
415421 const { mass } = massProperties(pl);
422+ const bounceHz = Math.sqrt(totalKw / mass) / (2 * Math.PI);
423+ const chain = feltChain(pad.lossFactor, bounceHz);
416424 const preloadTop = st.topPad ? st.preload : 0;
417- // The bottom pad carries the top pad's preload plus the weight of the plate.
418- const preloadBottom = preloadTop + (mass * GRAVITY) / (kArea * annulus);
419- return { kw, ku, kwBottom, nodes, eta: pad.lossFactor, preloadTop, preloadBottom, totalKw, bounceHz: Math.sqrt(totalKw / mass) / (2 * Math.PI) };
425+ // The bottom pad carries the top pad's preload plus the weight of the plate,
426+ // held by the relaxed stiffness of the felt.
427+ const preloadBottom = preloadTop + (mass * GRAVITY) / (kArea * annulus * chain.kInf);
428+ return { kw, ku, kwBottom, nodes, eta: pad.lossFactor, chain, preloadTop, preloadBottom, totalKw, bounceHz };
420429 }
421430
422431 // ---------- Patches ----------
src/sim/selftest.tsmodified+24−6View file
@@ -1,9 +1,10 @@
11 // In-browser check of the WGSL step kernel against the float64 reference:
22 // one step from a smooth random state of a few thicknesses amplitude, with
3-// the losses and the stand (lift-off included), compared as accelerations.
3+// the losses and the stand (lift-off and the felt's chain included),
4+// compared as accelerations, and the update of the felt's state.
45
56 import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode, standRates } from "../params.ts";
6-import { accelerations, addDamping, fields, fullLoss } from "../model/reference.ts";
7+import { accelerations, addDamping, addPadMemory, fields, fullLoss, padState, stepPadMemory } from "../model/reference.ts";
78 import { fitSigma, inPlaneRate } from "../loss.ts";
89 import { Simulation } from "./simulation.ts";
910
@@ -11,6 +12,8 @@ export interface SelfTestResult {
1112 n: number;
1213 /** Largest acceleration error over all nodes and components relative to the largest reference acceleration. */
1314 maxRel: number;
15+ /** Largest error of the felt's arm extensions after the step relative to their largest value (0 on a free plate). */
16+ padRel: number;
1417 dt: number;
1518 }
1619
@@ -38,7 +41,7 @@ export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParam
3841 const hm = meanThickness(pl);
3942 const fit = fitSigma(pl);
4043 const sr = standRates(pl);
41- const loss = fullLoss(fit, sr, inPlaneRate(pl));
44+ const loss = fullLoss(fit, inPlaneRate(pl));
4245 const sim = new Simulation(device, pl, params, loss);
4346 const amp = 2.5 * hm;
4447 const f = { u: smoothField(pl, 0.02 * amp, 1), v: smoothField(pl, 0.02 * amp, 2), w: smoothField(pl, amp, 3) };
@@ -48,17 +51,24 @@ export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParam
4851 pos32[k] = f.u[k]; pos32[n + k] = f.v[k]; pos32[2 * n + k] = f.w[k];
4952 vel32[k] = vel.u[k]; vel32[n + k] = vel.v[k]; vel32[2 * n + k] = vel.w[k];
5053 }
54+ // The felt's arms start with smooth random extensions of a fraction of the amplitude.
55+ const pad = padState(pl, sr);
56+ const nArms = sr ? sr.chain.tau.length : 0;
57+ for (let a = 0; a < 3 * nArms; a++) pad.set(smoothField(pl, 0.3 * amp, 20 + a), a * n);
58+ const pad32 = new Float32Array(pad);
5159 // Use the float32-rounded state as the reference input.
5260 for (let k = 0; k < n; k++) {
5361 f.u[k] = pos32[k]; f.v[k] = pos32[n + k]; f.w[k] = pos32[2 * n + k];
5462 vel.u[k] = vel32[k]; vel.v[k] = vel32[n + k]; vel.w[k] = vel32[2 * n + k];
5563 }
56- sim.writeState(pos32, vel32);
64+ pad.set(pad32);
65+ sim.writeState(pos32, vel32, pad32);
5766 await sim.stepOnce();
5867 const st = await sim.readState();
5968 const ref = fields(n);
6069 accelerations(pl, f, ref, params.padLiftoff, sr);
61- addDamping(pl, loss, vel, ref, sr);
70+ addDamping(pl, loss, vel, ref);
71+ addPadMemory(pl, sr, f, pad, params.padLiftoff, ref);
6272 let maxRef = 0, maxErr = 0;
6373 const dt = sim.dt;
6474 for (let k = 0; k < n; k++) {
@@ -66,6 +76,14 @@ export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParam
6676 maxRef = Math.max(maxRef, Math.abs(ref.u[k]), Math.abs(ref.v[k]), Math.abs(ref.w[k]));
6777 maxErr = Math.max(maxErr, Math.abs(au - ref.u[k]), Math.abs(av - ref.v[k]), Math.abs(aw - ref.w[k]));
6878 }
79+ // The felt's state after the step, from the GPU's own new velocities.
80+ let padRef = 0, padErr = 0;
81+ if (sr) {
82+ const dw = fields(n);
83+ for (let k = 0; k < n; k++) { dw.u[k] = dt * st.vel[k]; dw.v[k] = dt * st.vel[n + k]; dw.w[k] = dt * st.vel[2 * n + k]; }
84+ stepPadMemory(pl, sr, dt, dw, f, pad, params.padLiftoff);
85+ for (let i = 0; i < pad.length; i++) { padRef = Math.max(padRef, Math.abs(pad[i])); padErr = Math.max(padErr, Math.abs(pad[i] - st.pad[i])); }
86+ }
6987 sim.destroy();
70- return { n, maxRel: maxErr / maxRef, dt };
88+ return { n, maxRel: maxErr / maxRef, padRel: padRef > 0 ? padErr / padRef : 0, dt };
7189 }
src/sim/shaders/common.wgslmodified+13−2View file
@@ -35,13 +35,21 @@ struct Params {
3535 liftoff: u32, // 1 when the pads push only while compressed
3636 preloadTop: f32, // compression of the top pad at rest [m] (0 without a top pad)
3737 preloadBottom: f32, // compression of the bottom pad at rest [m]
38- padDamping: f32, // dashpot rate of the pads per unit of pad stiffness [s]
38+ kInf: f32, // relaxed stiffness of the felt, per unit of its stiffness at the bounce frequency
3939 hasTop: u32, // 1 when there is a top pad
40- pad0: u32,
40+ nArms: u32, // arms of the felt's chain (at most MAX_ARMS)
4141 pad1: u32,
4242 pad2: u32,
4343 pad3: u32,
44+ // The chain of the felt (felt.ts), four arms per vec4: decay factors
45+ // exp(-dt / tau) in [0..2], growth factors (tau / dt)(1 - decay) in
46+ // [3..5], and arm stiffnesses per unit of pad stiffness in [6..8].
47+ chain: array<vec4f, 9>,
4448 };
49+const MAX_ARMS = 12u;
50+fn armDecay(a: u32) -> f32 { return P.chain[a >> 2u][a & 3u]; }
51+fn armGrowth(a: u32) -> f32 { return P.chain[3u + (a >> 2u)][a & 3u]; }
52+fn armK(a: u32) -> f32 { return P.chain[6u + (a >> 2u)][a & 3u]; }
4553
4654 // Per-step parameters, addressed with a dynamic uniform offset.
4755 struct StepParams {
@@ -68,6 +76,9 @@ fn nodeY(j: i32) -> f32 { return (f32(j) - 0.5 * f32(P.ny - 1u)) * P.hy; }
6876 // stiffness at the node [N/m] (zero away from the pads; the top pad, when
6977 // present, has the same rate).
7078 @group(0) @binding(8) var<storage, read> nd: array<vec4f>;
79+// Extensions of the arms of the felt's chain: padQ[(field * nArms + arm) * n + k]
80+// with field 0, 1, 2 for u, v, w; zero is the relaxed felt.
81+@group(0) @binding(9) var<storage, read_write> padQ: array<f32>;
7182 fn nodeWeight(i: i32, j: i32) -> f32 {
7283 if (i < 0 || i >= i32(P.nx) || j < 0 || j >= i32(P.ny)) { return 0.0; }
7384 return nd[u32(j) * P.nx + u32(i)].y;
src/sim/shaders/step.wgslmodified+44−10View file
@@ -168,10 +168,13 @@ fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
168168 return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
169169 }
170170
171-// Force of the pads on a node with bottom-pad rate kb and transverse
172-// displacement w (same rule as padForce in reference.ts). The constant part
173-// that balances the weight of the plate is subtracted, so w = 0 is rest.
174-fn padForce(kb: f32, w: f32) -> f32 {
171+// Relaxed force of the pads on a node with bottom-pad rate kb (at the bounce
172+// frequency) and transverse displacement w (same rule as padForce in
173+// reference.ts); the rest of the felt's force is in the arms of its chain.
174+// The constant part that balances the weight of the plate is subtracted,
175+// so w = 0 is rest.
176+fn padForce(kb0: f32, w: f32) -> f32 {
177+ let kb = P.kInf * kb0;
175178 if (P.liftoff == 0u) {
176179 if (P.hasTop != 0u) { return -2.0 * kb * w; }
177180 return -kb * w;
@@ -180,6 +183,13 @@ fn padForce(kb: f32, w: f32) -> f32 {
180183 if (P.hasTop != 0u) { f -= kb * max(P.preloadTop + w, 0.0) - kb * P.preloadTop; }
181184 return f;
182185 }
186+// Number of pads in contact with a node at transverse displacement w.
187+fn padsInContact(w: f32) -> f32 {
188+ if (P.liftoff == 0u) { return select(1.0, 2.0, P.hasTop != 0u); }
189+ var nc = select(0.0, 1.0, P.preloadBottom - w > 0.0);
190+ if (P.hasTop != 0u && P.preloadTop + w > 0.0) { nc += 1.0; }
191+ return nc;
192+}
183193
184194 @compute @workgroup_size(16, 16)
185195 fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
@@ -249,15 +259,26 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve
249259 F += cellForce(i, j - 1, 3u);
250260 let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);
251261 let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
252- // The pads: a spring to ground on the transverse motion (compression-only
253- // with lift-off), the felt's shear on the in-plane motion, and a dashpot.
262+ // The pads: the relaxed spring to ground on the transverse motion
263+ // (compression-only with lift-off) and the felt's shear on the in-plane
264+ // motion, plus the arms of the felt's chain, whose extensions padQ carry
265+ // the memory of the felt (felt.ts, addPadMemory in reference.ts).
254266 let kb = nd[k].w;
255- if (P.standOn != 0u && kb > 0.0) {
267+ let onPad = P.standOn != 0u && kb > 0.0;
268+ var nc = 0.0;
269+ if (onPad) {
256270 var nPads = 1.0; if (P.hasTop != 0u) { nPads = 2.0; }
257271 let ku = nPads * kb * 0.3846153846; // shear: G / E = 1 / (2 (1 + 0.3))
258- F.x -= ku * pOld.x + P.padDamping * ku * vel.x;
259- F.y -= ku * pOld.y + P.padDamping * ku * vel.y;
260- F.z += padForce(kb, pOld.z) - P.padDamping * nPads * kb * vel.z;
272+ F.x -= P.kInf * ku * pOld.x;
273+ F.y -= P.kInf * ku * pOld.y;
274+ F.z += padForce(kb, pOld.z);
275+ nc = padsInContact(pOld.z);
276+ var arm = vec3f(0.0);
277+ for (var a = 0u; a < P.nArms; a++) {
278+ let base = a * n + k;
279+ arm += armK(a) * vec3f(padQ[base], padQ[P.nArms * n + base], padQ[2u * P.nArms * n + base]);
280+ }
281+ F -= vec3f(ku, ku, nc * kb) * arm;
261282 }
262283 let mass = P.rho * TH(i, j) * P.hx * P.hy * weight;
263284 var acc = F * (1.0 / mass);
@@ -323,6 +344,19 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve
323344 let vNew = vel + P.dt * acc;
324345 if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = vNew.z; }
325346 var pNew = pOld + P.dt * vNew;
347+ // The arms of the felt relax over the step while the node moves by dt vNew
348+ // (exact for a linear motion within the step); a node that has left every
349+ // pad relaxes its transverse arms at once.
350+ if (onPad) {
351+ let dw = P.dt * vNew;
352+ for (var a = 0u; a < P.nArms; a++) {
353+ let base = a * n + k;
354+ let iu = base; let iv = P.nArms * n + base; let iw = 2u * P.nArms * n + base;
355+ padQ[iu] = padQ[iu] * armDecay(a) + dw.x * armGrowth(a);
356+ padQ[iv] = padQ[iv] * armDecay(a) + dw.y * armGrowth(a);
357+ padQ[iw] = select(0.0, padQ[iw] * armDecay(a) + dw.z * armGrowth(a), nc > 0.0);
358+ }
359+ }
326360 // Compensated update of w: the increment dt * v can be far below the
327361 // float32 resolution of w, so the rounding error is carried separately.
328362 let y = P.dt * vNew.z + posIn[3u * n + k];
src/sim/simulation.tsmodified+30−11View file
@@ -23,6 +23,7 @@ import {
2323 listenerTerms,
2424 } from "../params.ts";
2525 import { type FullLoss, maxAngularFrequency } from "../model/reference.ts";
26+import { FELT_ARMS, feltStepConstants } from "../felt.ts";
2627 import commonWgsl from "./shaders/common.wgsl?raw";
2728 import stepWgsl from "./shaders/step.wgsl?raw";
2829 import audioWgsl from "./shaders/audio.wgsl?raw";
@@ -52,7 +53,7 @@ export interface SimStatus {
5253 stick: StickState;
5354 }
5455
55-const PARAMS_LEN = 40;
56+const PARAMS_LEN = 40 + 9 * 4; // Params: 40 scalars and the 9 vec4 of the felt's chain
5657 const STEP_LEN = 12;
5758 const STICK_LEN = 16;
5859 /** Nodes per workgroup along each axis in step.wgsl. */
@@ -100,6 +101,7 @@ export class Simulation {
100101 private hist: GPUBuffer;
101102 private stickBuf: GPUBuffer;
102103 private nd: GPUBuffer;
104+ private padQ: GPUBuffer;
103105 private lw: GPUBuffer;
104106 private audioOut: GPUBuffer;
105107 private frameBuf: GPUBuffer;
@@ -145,6 +147,9 @@ export class Simulation {
145147 }
146148 this.nd = f32(4 * n);
147149 device.queue.writeBuffer(this.nd, 0, packed);
150+ const nArms = this.stand ? this.stand.chain.tau.length : 0;
151+ if (nArms > FELT_ARMS) throw new Error("the felt's chain has more arms than the kernel holds");
152+ this.padQ = f32(3 * nArms * n);
148153 this.lw = f32(2 * n);
149154 this.audioOut = f32(this.samplesPerBatch);
150155 this.frameBuf = f32(n);
@@ -179,8 +184,13 @@ export class Simulation {
179184 pu[30] = st && st.nodes.length > 0 ? 1 : 0;
180185 pu[31] = params.padLiftoff ? 1 : 0;
181186 pf[32] = st ? st.preloadTop : 0; pf[33] = st ? st.preloadBottom : 0;
182- pf[34] = loss.padDamping;
187+ pf[34] = st ? st.chain.kInf : 1;
183188 pu[35] = plate.stand?.topPad ? 1 : 0;
189+ pu[36] = nArms;
190+ if (st) {
191+ const { decay, growth } = feltStepConstants(st.chain, this.dt);
192+ for (let a = 0; a < nArms; a++) { pf[40 + a] = decay[a]; pf[52 + a] = growth[a]; pf[64 + a] = st.chain.k[a]; }
193+ }
184194 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
185195 device.queue.writeBuffer(this.paramsBuf, 0, pf);
186196
@@ -195,7 +205,7 @@ export class Simulation {
195205 const pipeline = (module: GPUShaderModule, layout: GPUBindGroupLayout) =>
196206 device.createComputePipeline({ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }), compute: { module, entryPoint: "main" } });
197207 const stepLayout = device.createBindGroupLayout({
198- entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage"), entry(5, "storage"), entry(6, "storage"), entry(7, "storage"), entry(8, "read-only-storage")],
208+ entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage"), entry(5, "storage"), entry(6, "storage"), entry(7, "storage"), entry(8, "read-only-storage"), entry(9, "storage")],
199209 });
200210 this.stepPipe = pipeline(shader(stepWgsl), stepLayout);
201211 const audioLayout = device.createBindGroupLayout({
@@ -219,6 +229,7 @@ export class Simulation {
219229 { binding: 6, resource: { buffer: this.hist } },
220230 { binding: 7, resource: { buffer: this.stickBuf } },
221231 { binding: 8, resource: { buffer: this.nd } },
232+ { binding: 9, resource: { buffer: this.padQ } },
222233 ],
223234 }),
224235 );
@@ -410,31 +421,39 @@ export class Simulation {
410421 this.totalSteps++;
411422 await this.device.queue.onSubmittedWorkDone();
412423 }
413- /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw]. */
414- writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>): void {
424+ /** Length of the felt's state: 3 fields times the arms of the chain times the nodes (0 on a free plate). */
425+ get padStateLength(): number {
426+ return 3 * (this.stand ? this.stand.chain.tau.length : 0) * this.plate.grid.n;
427+ }
428+ /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw], and the extensions of the felt's arms (relaxed when omitted). */
429+ writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>, pad?: Float32Array<ArrayBuffer>): void {
415430 const n = this.plate.grid.n;
416431 if (pos.length !== 4 * n || vel.length !== 3 * n) throw new Error("state arrays have the wrong length");
417432 this.device.queue.writeBuffer(this.pos[this.cur], 0, pos);
418433 this.device.queue.writeBuffer(this.vel[this.cur], 0, vel);
434+ const np = this.padStateLength;
435+ if (pad !== undefined && pad.length !== np) throw new Error("the felt state has the wrong length");
436+ if (np > 0) this.device.queue.writeBuffer(this.padQ, 0, pad ?? new Float32Array(np));
419437 }
420- /** Reads back the full state as [u | v | w | wLo] and [du | dv | dw]. */
421- async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
422- const n = this.plate.grid.n;
423- const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
438+ /** Reads back the full state as [u | v | w | wLo], [du | dv | dw], and the felt's arm extensions. */
439+ async readState(): Promise<{ pos: Float32Array; vel: Float32Array; pad: Float32Array }> {
440+ const n = this.plate.grid.n, np = this.padStateLength;
441+ const stg = this.device.createBuffer({ size: (7 * n + np) * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
424442 const enc = this.device.createCommandEncoder();
425443 enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
426444 enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
445+ if (np > 0) enc.copyBufferToBuffer(this.padQ, 0, stg, 7 * n * 4, np * 4);
427446 this.device.queue.submit([enc.finish()]);
428447 await stg.mapAsync(GPUMapMode.READ);
429448 const all = new Float32Array(stg.getMappedRange()).slice();
430449 stg.unmap();
431450 stg.destroy();
432- return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
451+ return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n, 7 * n), pad: all.subarray(7 * n) };
433452 }
434453
435454 destroy(): void {
436455 this.inFlight = null;
437- for (const b of [...this.pos, ...this.vel, this.hist, this.stickBuf, this.nd, this.lw, this.audioOut, this.frameBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
456+ for (const b of [...this.pos, ...this.vel, this.hist, this.stickBuf, this.nd, this.padQ, this.lw, this.audioOut, this.frameBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
438457 }
439458 get fs(): number { return this.params.fs; }
440459 get fps(): number { return this.params.fps; }
src/tutorial/steps.tsmodified+1−1View file
@@ -635,7 +635,7 @@ ${live(`Window of ${N} samples at ${fs / 1000} kHz: $\\Delta f$ = ${(fs / N).toF
635635 return `
636636 <p>A cymbal rests on a felt pad around the rod of its stand, usually with a second pad above it under a wing nut. Felt is a soft spring: a pad of thickness $t$ and Young's modulus $E_p$ has a stiffness per unit area of $E_p/t$ in compression, and an annulus of felt between the radii $r_1$ and $r_2$ resists a vertical motion of the cymbal with the stiffness</p>
637637 $$ K = n_p\\,\\frac{E_p}{t}\\,\\pi\\left(r_2^2 - r_1^2\\right), $$
638-<p>with $n_p$ the number of pads. The cymbal as a whole then bounces on the felt at $f_b = \\frac{1}{2\\pi}\\sqrt{K/M}$, and it rocks, tilting about a diameter, against the torque of the annulus, whose stiffness is its second moment $K_\\theta = n_p (E_p/t)\\,\\pi(r_2^4 - r_1^4)/4$, at $f_r = \\frac{1}{2\\pi}\\sqrt{K_\\theta/I}$ with $I$ the moment of inertia of the cymbal about a diameter. These are modes that a free plate does not have, below every plate mode: the resonances of the stand, a few hertz for the rock and tens of hertz for the bounce. The felt is lossy, and its loss factor $\\eta$ adds to each mode a decay rate proportional to the fraction of the mode's energy stored in the felt, $\\sigma = \\eta\\,\\phi^T K_s\\phi/(2\\omega)$; the rock and the bounce die within a fraction of a second, while the plate modes, which hardly move the center, are barely affected.</p>
638+<p>with $n_p$ the number of pads. The cymbal as a whole then bounces on the felt at $f_b = \\frac{1}{2\\pi}\\sqrt{K/M}$, and it rocks, tilting about a diameter, against the torque of the annulus, whose stiffness is its second moment $K_\\theta = n_p (E_p/t)\\,\\pi(r_2^4 - r_1^4)/4$, at $f_r = \\frac{1}{2\\pi}\\sqrt{K_\\theta/I}$ with $I$ the moment of inertia of the cymbal about a diameter. These are modes that a free plate does not have, below every plate mode: the resonances of the stand, a few hertz for the rock and tens of hertz for the bounce. The felt is lossy, and its loss factor $\\eta$ adds to each mode a decay rate proportional to the fraction of the mode's energy stored in the felt, $\\sigma = \\eta\\,\\phi^T K_s\\phi/(2\\omega)$; the rock and the bounce die within a fraction of a second, while the plate modes, which hardly move the center, are barely affected. A loss factor that is the same at every frequency is not what a spring and a dashpot give (a dashpot's loss grows with frequency), so the felt is modeled as a viscoelastic solid, a chain of springs and dashpots with relaxation times spread over the band, which the full model steps in time and the linear model evaluates at each mode's frequency; its modulus is then the value at the bounce, and it is a little stiffer above and softer below.</p>
639639 <p>In the linear model the pads are linear springs. The full model lets them do what real felt does, push only while compressed, so that a hard hit can lift the cymbal off the bottom pad. The rod and its sleeve are not modeled; the shear of the felt is what holds the cymbal sideways. The stand is a setting of the setup rather than of the design: the Stand parameters beside the script put any plate on it or take it off, and set the felt's radii, thickness, modulus and loss factor, the preload of the wing nut and the position of the rod.</p>
640640 ${liveText ? live(liveText) : need("Turn the stand on (the Stand parameters, or the button below) to see it.")}
641641 <p>The rock and the bounce are found among the computed modes by their rigid share, the fraction of a mode's kinetic energy that is a rigid translation or a rigid tilt of the cymbal: $(\\sum_k m_k\\phi_k)^2/M$ for the translation, and the same with the tilt. The rock comes out near $f_r$: at a few hertz it lies far below the first plate mode, and the cymbal tilts nearly as a whole. The bounce comes out well below $f_b$, and the reason is instructive. The formula treats the cymbal as rigid, but a thin plate is not: the pads push on its center, and, loosely speaking, the compliance of the plate around the center adds in series with that of the felt, so the center yields and the rim follows late. The bounce mode is a rigid translation blended with the umbrella-shaped plate mode, one nodal circle and no nodal diameter, whose antinode is at the center where the pads act; and the higher plate modes that move the center take a share of the felt's damping, which shortens their decay. A thicker plate is stiffer at its center and its bounce comes closer to $f_b$: compare the 20 inch ride with the flat disc.</p>`;
tests/assemble.test.tsmodified+4−1View file
@@ -41,7 +41,10 @@ test("K x equals the linearized reference force on a bowed, tapered plate with a
4141 const g = withMask(makeGrid(0.2, 0.2, 0.0125), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
4242 const w0 = new Float64Array(g.n);
4343 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = 3e-3 * Math.sin(20 * nodeX(g, i)) * Math.cos(15 * nodeY(g, j));
44- const pl: Plate = { ...uniformPlate(g, bronze, 1e-3, w0), stand: { x: 0, y: 0, pad: { innerRadius: 0.012, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
44+ // A lossless felt: with a loss factor the reference's static pad force is the relaxed
45+ // stiffness, and the rest of the linear model's stiffness at the bounce sits in the
46+ // arms of the felt's chain (reference.test.ts checks those).
47+ const pl: Plate = { ...uniformPlate(g, bronze, 1e-3, w0), stand: { x: 0, y: 0, pad: { innerRadius: 0.012, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0 }, topPad: true, preload: 1e-3 } };
4548 for (let k = 0; k < g.n; k++) if (pl.h[k] > 0) pl.h[k] = 1e-3 * (1 + 0.4 * Math.abs(nodeX(g, k % g.nx)) / 0.1);
4649 assert.ok(!isFlatPlate(pl));
4750 const sys = assemble(pl);
tests/design.test.tsmodified+3−1View file
@@ -112,7 +112,9 @@ design.inside = @(x, y, r, theta) (r <= 0.15) & (r >= 0.0065);
112112 const kPad = (1e6 / 0.01) * annulus;
113113 assert.ok(Math.abs(sr.totalKw / (2 * kPad) - 1) < 1e-9, `total ${sr.totalKw} vs ${2 * kPad}`);
114114 const { mass } = massProperties(pl);
115- assert.ok(Math.abs(sr.preloadBottom - (0.5e-3 + (mass * 9.81) / kPad)) < 1e-12);
115+ // The weight rests on the relaxed felt, which is softer than the felt at the bounce.
116+ assert.ok(sr.chain.kInf < 1 && sr.chain.kInf > 0);
117+ assert.ok(Math.abs(sr.preloadBottom - (0.5e-3 + (mass * 9.81) / (kPad * sr.chain.kInf))) < 1e-12);
116118 assert.equal(sr.preloadTop, 0.5e-3);
117119 assert.ok(Math.abs(sr.bounceHz - Math.sqrt(sr.totalKw / mass) / (2 * Math.PI)) < 1e-9);
118120 const pl2 = evaluatePlate(d, 0.01, { ...stand, topPad: false });
tests/felt.test.tsadded+65−0View file
@@ -0,0 +1,65 @@
1+// The viscoelastic felt: its chain of Maxwell arms has the loss factor it
2+// was built for across the audio band, the storage stiffness it was
3+// normalized to, and an exact time step.
4+import { test } from "node:test";
5+import assert from "node:assert/strict";
6+import { feltChain, feltModulus, feltStepConstants, FELT_ARMS } from "../src/felt.ts";
7+
8+test("the loss factor of the chain is flat across the band and the storage is one at the reference", () => {
9+ for (const eta of [0.05, 0.3, 0.6]) {
10+ for (const fRef of [15, 96, 400]) {
11+ const ch = feltChain(eta, fRef);
12+ assert.equal(ch.tau.length, FELT_ARMS);
13+ const ref = feltModulus(ch, 2 * Math.PI * fRef);
14+ assert.ok(Math.abs(ref.storage - 1) < 1e-12);
15+ assert.ok(Math.abs(ref.loss / ref.storage / eta - 1) < 0.02, `eta ${eta} at the reference ${fRef} Hz: ${ref.loss / ref.storage}`);
16+ let worst = 0, prev = 0;
17+ for (let f = 5; f <= 10000; f *= 1.1) {
18+ const m = feltModulus(ch, 2 * Math.PI * f);
19+ worst = Math.max(worst, Math.abs(m.loss / m.storage / eta - 1));
20+ assert.ok(m.storage > prev, "the storage stiffness rises with frequency");
21+ prev = m.storage;
22+ }
23+ const tol = eta > 0.5 ? 0.1 : 0.06;
24+ assert.ok(worst < tol, `eta ${eta}, reference ${fRef} Hz: loss factor off by up to ${(100 * worst).toFixed(1)}% between 5 Hz and 10 kHz`);
25+ // Relaxed and glassy stiffness bracket the reference.
26+ const glassy = ch.kInf + ch.k.reduce((a, b) => a + b, 0);
27+ assert.ok(ch.kInf < 1 && glassy > 1);
28+ }
29+ }
30+ const none = feltChain(0, 100);
31+ assert.equal(none.tau.length, 0);
32+ assert.equal(none.kInf, 1);
33+ assert.deepEqual(feltModulus(none, 1000), { storage: 1, loss: 0 });
34+});
35+
36+test("the exact step of an arm reproduces its response to a ramp and to a sinusoid", () => {
37+ const ch = feltChain(0.3, 50);
38+ const a = 6; // an arm in the middle of the band
39+ const tau = ch.tau[a];
40+ // A ramp w = s t from a relaxed arm: p(t) = s tau (1 - exp(-t / tau)), whatever the step.
41+ for (const dt of [tau / 50, tau / 3, 2 * tau]) {
42+ const { decay, growth } = feltStepConstants(ch, dt);
43+ const s = 1e-3;
44+ let p = 0;
45+ const steps = 20;
46+ for (let i = 0; i < steps; i++) p = p * decay[a] + s * dt * growth[a];
47+ const exact = s * tau * (1 - Math.exp(-(steps * dt) / tau));
48+ assert.ok(Math.abs(p / exact - 1) < 1e-12, `dt = ${dt / tau} tau: ${p} vs ${exact}`);
49+ }
50+ // A sinusoid w = A sin(omega t), stepped finely: the work per cycle against
51+ // the arm's force -k p is pi k A^2 omega tau / (1 + omega^2 tau^2).
52+ const omega = 1 / tau, A = 1e-4, steps = 2000;
53+ const dt = (2 * Math.PI) / omega / steps;
54+ const { decay, growth } = feltStepConstants(ch, dt);
55+ let p = 0, work = 0;
56+ for (let i = 0; i < 6 * steps; i++) {
57+ const w0 = A * Math.sin(omega * i * dt), w1 = A * Math.sin(omega * (i + 1) * dt);
58+ const f0 = ch.k[a] * p;
59+ p = p * decay[a] + (w1 - w0) * growth[a];
60+ const f1 = ch.k[a] * p;
61+ if (i >= 5 * steps) work += 0.5 * (f0 + f1) * (w1 - w0);
62+ }
63+ const expected = (Math.PI * ch.k[a] * A * A * (omega * tau)) / (1 + omega * omega * tau * tau);
64+ assert.ok(Math.abs(work / expected - 1) < 1e-3, `work per cycle ${work} vs ${expected}`);
65+});
tests/gpu.test.tsmodified+45−1View file
@@ -57,9 +57,10 @@ test("the WGSL step agrees with the float64 reference on several presets", async
5757 const text = await page.evaluate(() => (window as unknown as { cymbalsim: { selftest: () => Promise<string> } }).cymbalsim.selftest());
5858 console.log(text);
5959 for (const line of text.split("\n")) {
60- const m = /max relative error ([0-9.e+-]+)/.exec(line);
60+ const m = /max relative error ([0-9.e+-]+), felt state error ([0-9.e+-]+)/.exec(line);
6161 assert.ok(m, line);
6262 assert.ok(Number(m![1]) < 2e-4, line);
63+ assert.ok(Number(m![2]) < 1e-5, line);
6364 }
6465 assert.deepEqual(errors, []);
6566 });
@@ -323,6 +324,49 @@ test("a light hit sounds the same in the linear and the full model", async () =>
323324 });
324325 });
325326
327+test("the felt damps a light hit on the aluminum square alike in the linear and the full model", async () => {
328+ // The square on the default stand puts a share of many modes' energy into the felt.
329+ // What the stand adds to the decay of each band, measured as the band's drop between
330+ // 0.25 and 0.55 s with the stand on minus the same drop with the stand off, must agree
331+ // between the models: the difference cancels the fitted plate loss of the full model
332+ // against the per-mode loss of the linear one, and leaves the felt. (The first 0.25 s
333+ // are left out: the modes that hold a large share of their stiffness in the felt die
334+ // within it, and the linear model's first-order rates over-estimate their damping.)
335+ const fs = 48000, N = 8192;
336+ const bands: [number, number][] = [[50, 200], [200, 600], [600, 2000], [2000, 6000]];
337+ const bandEnergy = (x: number[], t0: number, lo: number, hi: number) => {
338+ const n0 = Math.round(t0 * fs);
339+ const re = new Float64Array(N), im = new Float64Array(N);
340+ for (let i = 0; i < N; i++) re[i] = x[n0 + i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / N));
341+ fftInPlace(re, im);
342+ let e = 0;
343+ for (let k = Math.round((lo / fs) * N); k < Math.round((hi / fs) * N); k++) e += re[k] * re[k] + im[k] * im[k];
344+ return e;
345+ };
346+ const drops = (x: number[]) => bands.map(([lo, hi]) => 10 * Math.log10(bandEnergy(x, 0.25, lo, hi) / bandEnergy(x, 0.55, lo, hi)));
347+ const result: Record<string, { linear: number[]; full: number[] }> = {};
348+ let skipped = false;
349+ for (const stand of ["1", "0"]) {
350+ await withPage(`preset=aluminum-square&duration=0.8&strikeVelocity=0.05&spacing=0.012&fps=60&standOn=${stand}`, async (page, errors) => {
351+ if (!(await hasGpu(page))) { console.log("no WebGPU in this Chromium; skipping"); skipped = true; return; }
352+ const grab = () => page.evaluate(() => { const a = (window as unknown as { cymbalsim: { app: { rec: { audio: Float32Array; samples: number } } } }).cymbalsim.app; return Array.from(a.rec.audio.subarray(0, a.rec.samples)); });
353+ await page.evaluate(() => (window as unknown as { cymbalsim: { runLinear: () => Promise<void> } }).cymbalsim.runLinear());
354+ const linear = await grab();
355+ await page.evaluate(() => (window as unknown as { cymbalsim: { runFull: () => Promise<void> } }).cymbalsim.runFull());
356+ const full = await grab();
357+ assert.equal(linear.length, full.length);
358+ result[stand] = { linear: drops(linear), full: drops(full) };
359+ assert.deepEqual(errors, []);
360+ });
361+ if (skipped) return;
362+ }
363+ bands.forEach(([lo, hi], b) => {
364+ const dl = result["1"].linear[b] - result["0"].linear[b], df = result["1"].full[b] - result["0"].full[b];
365+ console.log(`${lo} to ${hi} Hz: decay over 0.3 s added by the stand, linear ${dl.toFixed(1)} dB (${result["0"].linear[b].toFixed(1)} free), full ${df.toFixed(1)} dB (${result["0"].full[b].toFixed(1)} free)`);
366+ assert.ok(Math.abs(dl - df) < 4, `${lo} to ${hi} Hz: the stand adds ${dl} dB in the linear model and ${df} dB in the full model`);
367+ });
368+});
369+
326370 function fftInPlace(re: Float64Array, im: Float64Array): void {
327371 const n = re.length;
328372 for (let i = 1, j = 0; i < n; i++) {
tests/reference.test.tsmodified+53−9View file
@@ -4,9 +4,10 @@
44 import { test } from "node:test";
55 import assert from "node:assert/strict";
66 import { makeGrid, withMask, type Grid, type Plate, nodeMasses, standRates, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
7-import { uniformPlate, accelerations, addDamping, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
7+import { uniformPlate, accelerations, addDamping, addPadMemory, padState, stepPadMemory, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
88 import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
99 import { fullLoss } from "../src/model/reference.ts";
10+import { feltModulus } from "../src/felt.ts";
1011
1112 const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
1213
@@ -111,22 +112,62 @@ test("forces are the exact gradient of the discrete energy (flat and bowed, tape
111112 }
112113 });
113114
114-test("the loss terms and the pad dashpots only remove energy", () => {
115+test("the loss terms only remove energy", () => {
115116 const g = makeGrid(0.1, 0.08, 0.01);
116- const pl = withStand(uniformPlate(g, steel, 1e-3));
117- const sr = standRates(pl);
118- const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, sr, 20);
117+ const pl = uniformPlate(g, steel, 1e-3);
118+ const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, 20);
119119 const masses = nodeMasses(pl);
120120 for (const seed of [1, 2, 3]) {
121121 const vel = randomFields(pl, 1, seed);
122122 const acc = fields(g.n);
123- addDamping(pl, loss, vel, acc, sr);
123+ addDamping(pl, loss, vel, acc);
124124 let power = 0;
125125 for (let k = 0; k < g.n; k++) power += masses[k] * (acc.u[k] * vel.u[k] + acc.v[k] * vel.v[k] + acc.w[k] * vel.w[k]);
126126 assert.ok(power < 0, `power ${power}`);
127127 }
128128 });
129129
130+test("the felt's chain, stepped in time, dissipates the work of its loss stiffness per cycle", () => {
131+ const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
132+ const pl = withStand(uniformPlate(g, steel, 1e-3));
133+ const sr = standRates(pl)!;
134+ const masses = nodeMasses(pl);
135+ // A rigid translation w = A sin(omega t) of the plate, stepped finely; the work done
136+ // against the arms over a cycle, once the transient of the arms has passed, is
137+ // pi L(omega) K A^2 with L the loss stiffness per unit of K (eta at the bounce).
138+ for (const f of [sr.bounceHz, 2000]) {
139+ const omega = 2 * Math.PI * f, A = 1e-4;
140+ const steps = 400, dt = 1 / (f * steps);
141+ const pad = padState(pl, sr);
142+ const x = fields(g.n), dw = fields(g.n), acc = fields(g.n);
143+ const setW = (t: number, out: Fields) => { for (let k = 0; k < g.n; k++) out.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? A * Math.sin(omega * t) : 0; };
144+ let work = 0;
145+ const cycles = 4;
146+ for (let s = 0; s < cycles * steps; s++) {
147+ const t0 = s * dt, t1 = (s + 1) * dt;
148+ setW(t0, x); setW(t1, dw);
149+ for (let k = 0; k < g.n; k++) dw.w[k] -= x.w[k];
150+ acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
151+ addPadMemory(pl, sr, x, pad, false, acc);
152+ // Work of the plate against the arms over the step, with the force at mid-step (trapezoid).
153+ let F0 = 0;
154+ for (let k = 0; k < g.n; k++) F0 += masses[k] * acc.w[k];
155+ stepPadMemory(pl, sr, dt, dw, x, pad, false);
156+ setW(t1, x);
157+ acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
158+ addPadMemory(pl, sr, x, pad, false, acc);
159+ let F1 = 0;
160+ for (let k = 0; k < g.n; k++) F1 += masses[k] * acc.w[k];
161+ if (s >= (cycles - 1) * steps) work -= 0.5 * (F0 + F1) * dw.w[sr.nodes[0]];
162+ }
163+ const L = feltModulus(sr.chain, omega).loss;
164+ const expected = Math.PI * L * sr.totalKw * A * A;
165+ assert.ok(work > 0);
166+ assert.ok(Math.abs(work / expected - 1) < 0.02, `at ${f.toFixed(0)} Hz: work per cycle ${work} vs ${expected}`);
167+ if (f === sr.bounceHz) assert.ok(Math.abs(L / sr.eta - 1) < 0.02, `loss stiffness at the bounce ${L} vs eta ${sr.eta}`);
168+ }
169+});
170+
130171 test("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
131172 const g = makeGrid(0.1, 0.1, 0.01);
132173 const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
@@ -152,14 +193,17 @@ test("the stand resists a rigid translation with its total stiffness and a tilt
152193 accelerations(pl, f, acc, false, sr);
153194 let F = 0;
154195 for (let k = 0; k < g.n; k++) F += masses[k] * acc.w[k];
155- assert.ok(Math.abs(-F / dz / sr.totalKw - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${sr.totalKw}`);
196+ // The relaxed felt holds kInf of the stiffness at the bounce; the arms of the chain hold the rest.
197+ const kInf = sr.chain.kInf;
198+ assert.ok(kInf > 0.2 && kInf < 0.8, `relaxed fraction ${kInf}`);
199+ assert.ok(Math.abs(-F / dz / (kInf * sr.totalKw) - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${kInf * sr.totalKw}`);
156200 // Tilt about y: w = theta x; restoring moment = theta * sum_k kw_k x_k^2 (bending of a plane is
157201 // zero; the von Karman membrane strain of a rotation is theta^2 / 2, so a tiny tilt keeps it negligible).
158202 const th = 1e-7;
159203 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? th * nodeX(g, k % g.nx) : 0;
160204 accelerations(pl, f, acc, false, sr);
161205 let M = 0, I2 = 0;
162- for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
206+ for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += kInf * sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
163207 assert.ok(Math.abs(-M / th / I2 - 1) < 1e-6, `tilt stiffness ${-M / th} vs ${I2}`);
164208 // With lift-off, a small motion within the preload is the same linear spring; a large upward motion unloads the bottom pad.
165209 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 1e-5 : 0;
@@ -172,7 +216,7 @@ test("the stand resists a rigid translation with its total stiffness and a tilt
172216 let Flo = 0;
173217 for (let k = 0; k < g.n; k++) Flo += masses[k] * lo.w[k];
174218 // Bottom pad unloaded: only the top pad (compressed by preload + 10 mm) pushes down, minus the rest balance.
175- const kb = sr.totalKw / 2;
219+ const kb = (kInf * sr.totalKw) / 2;
176220 const expected = -(kb * (sr.preloadTop + 0.01) - kb * sr.preloadTop) - kb * sr.preloadBottom;
177221 assert.ok(Math.abs(Flo / expected - 1) < 1e-6, `lift-off force ${Flo} vs ${expected}`);
178222 });