/ proof-of-concept / cymbalsim
Sign in
proof-of-concept / cymbalsim
540490c
Branches
main
cymbalsim / THEORY.md
137 lines · 33.6 KBCodeBlameHistory
3This note describes the model behind cymbalsim: a thin metal plate of arbitrary outline and spatially varying thickness, flat or with a stress-free curved profile, either hanging freely or sitting on a cymbal stand, struck once with a drumstick. The sound is computed in one of two ways. The *linear model* takes the equations of the plate linearized about its rest shape, computes their modes once per design, and synthesizes the strike from the modes; this takes a few seconds and needs no GPU. The *full model* steps the geometrically nonlinear plate equations of Foppl and von Karman explicitly on the GPU; this takes tens of seconds per second of sound and captures what the linear model cannot, the pitch glide and the transfer of energy toward high frequencies that turn the initial tone of a hard hit into a wash. Both use the same spatial discretization, so the linear model is exactly the small-amplitude limit of the full one, and a light hit sounds the same either way.
5Section 1 gives the plate model, Section 2 its discretization on the grid, Section 3 the stand, Section 4 the linear model and the modal synthesis, Section 5 the full model, Section 6 the checks that were made, and Section 7 the limitations. The discretization and the GPU kernel follow those of steelsim41, the earlier simulator of a scripted cymbal, simplified where the modal approach allows it; what carries over is repeated here so that this note stands alone.
7## 1. Plate model
9The plate occupies a region $\Omega$ of the plane, given by the design script as a bounding box and an outline function (a rectangle, a disc or an ellipse with or without a center hole, or any other test on the point). It has a thickness $h(x,y)$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$; the material constants are uniform, the thickness need not be. Its motion is described by the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$, all measured from a stress-free shape $w_0(x,y)$, the profile of the design, which is zero for a flat plate. The membrane strains in the Marguerre form are
11$$
12\varepsilon_{xx} = u_x + \tfrac12\left[(w_0 + w)_x^2 - w_{0,x}^2\right], \qquad
13\varepsilon_{yy} = v_y + \tfrac12\left[(w_0 + w)_y^2 - w_{0,y}^2\right], \qquad
14\gamma_{xy} = u_y + v_x + (w_0 + w)_x (w_0 + w)_y - w_{0,x} w_{0,y},
15$$
17so that $u = v = w = 0$ carries no stress. The membrane forces per unit length are $N_{xx} = C(\varepsilon_{xx} + \nu\varepsilon_{yy})$, $N_{yy} = C(\varepsilon_{yy} + \nu\varepsilon_{xx})$, and $N_{xy} = G\gamma_{xy}$, with $C = Eh/(1-\nu^2)$ and $G = Eh/(2(1+\nu))$. Writing $W = w_0 + w$ for the total shape and $D = Eh^3/(12(1-\nu^2))$ for the local bending stiffness, the equations of motion are
19$$
20\rho h\, u_{tt} = \partial_x N_{xx} + \partial_y N_{xy}, \qquad
21\rho h\, v_{tt} = \partial_x N_{xy} + \partial_y N_{yy},
22$$
24$$
25\rho h\, w_{tt} = -\mathcal{B}[w] + \partial_x(N_{xx} W_x + N_{xy} W_y) + \partial_y(N_{xy} W_x + N_{yy} W_y) + \text{(losses)} ,
26$$
28where $\mathcal{B}[w]$ is the variable-stiffness bending operator, the Euler-Lagrange derivative of the bending energy
30$$
31U_b = \frac12 \int_\Omega D\left[w_{xx}^2 + w_{yy}^2 + 2\nu\, w_{xx}w_{yy} + 2(1-\nu)\, w_{xy}^2\right] dA ,
32$$
34which reduces to $D\nabla^4 w$ where $D$ is constant. Bending acts on the deviation $w$ from the stress-free shape; the membrane forces act through the total slope. For $w_0 = 0$ and constant $h$ these are the von Karman equations.
36**The linear part.** Dropping the terms of second order in the displacements, the strains become $\varepsilon_{xx} = u_x + w_{0,x} w_x$, $\varepsilon_{yy} = v_y + w_{0,y} w_y$, $\gamma_{xy} = u_y + v_x + w_{0,x} w_y + w_{0,y} w_x$, and the transverse equation keeps $\partial_x(N_{xx} w_{0,x} + N_{xy} w_{0,y}) + \partial_y(\ldots)$ with the linear $N$. For a flat plate ($w_0 = 0$) this decouples completely: the transverse motion is the Kirchhoff plate $\rho h w_{tt} = -\mathcal{B}[w]$, and the in-plane motion is neither excited by a transverse strike nor heard, so the linear model of a flat plate is a bending problem alone. A profile couples the two through its slopes: a curved plate stretches when it bends, which raises the low modes, and the coupling is what makes a shell differ from a plate.
38**The nonlinear part.** The quadratic terms in the strains couple bending to stretching at finite amplitude. A flat plate has a symmetric, hardening nonlinearity; a curved plate has an asymmetric one, may snap between curved states, and softens as it passes through the flat configuration. These terms produce the pitch glide of a strongly struck cymbal and the cascade of energy from the low modes toward high frequencies. They are kept by the full model and dropped by the linear one.
40**Thickness field.** A real cymbal is lathed to a taper, thick at the bell and thin at the edge. The thickness enters in three places: the mass per unit area $\rho h$, the membrane stiffnesses $C$ and $G$, proportional to $h$, and the bending stiffness $D$, proportional to $h^3$. A plate of twice the thickness has twice the bending frequencies and the same in-plane frequencies. The thickness must be positive over the plate and should vary slowly on the scale of the grid.
42**Profile.** The stress-free shape is given by the script as a function of position. Only its slopes enter the strains, so its mean is immaterial, but its tilt is not: the von Karman strains are not invariant under rotation, and a tilted stress-free shape would add a spurious membrane coupling linear in the tilt. The mass-weighted least-squares plane is therefore subtracted from every profile. The treatment is that of a shallow shell, valid while the slopes are small compared with one; the interface reports the largest slope and warns above 0.1. The bell of a real crash rises 25 mm over a radius of 50 mm, a slope of about 0.5, so the curved preset uses a flatter bell, and this is the main point of tension between the model and a real cymbal (Section 7).
44**Losses.** The design gives the decay time $T_{60}(f)$, the time in which a partial at frequency $f$ decays by 60 dB, as a number, a table interpolated in log-log coordinates, or a function of frequency. The decay rate of a partial is then $\sigma(f) = \ln(1000)/T_{60}(f)$, so that its amplitude decays as $e^{-\sigma t}$. The linear model applies this rate to each mode at its own frequency (Section 4). The full model needs a loss operator rather than a rate per mode, and uses the standard two-parameter model of plate synthesis, a rate $\sigma_0 + \sigma_1 k^2$ for a bending wave of wavenumber $k$, with $\sigma_0$ and $\sigma_1$ fitted to the design's curve (Section 5). None of these values are measured; the presets are chosen to sound plausible, and the fit is reported next to the requested curve so that the difference is visible.
46## 2. Spatial discretization
48The bounding box is covered by a regular grid of $n_x \times n_y$ nodes with spacings $h_x$ and $h_y$ (8 mm by default), and the outline selects the part of the grid that is the plate. All three displacements live at the nodes, and so does the thickness. The forces are the negative gradient of a discrete potential energy $U = U_b + U_m + U_s$ (bending, membrane, stand), so the semi-discrete system is Hamiltonian, its linear part has a symmetric stiffness matrix, and the free edge needs no special formulas.
50**Outline.** A cell belongs to the plate when its center lies inside the outline, and a node when it is a corner of such a cell; the edge is therefore a staircase of whole cells, and a hole is simply a region whose cells are absent. Every rule below refers to this membership: a second difference exists where the node and both of its neighbors belong to the plate, a cell contributes membrane energy and the mixed bending term when its four corners do, and the weight $\omega_k$ of a node is the fraction of its four cells that belong to the plate (1 inside, 1/2 on an edge, 1/4 at a corner of a full rectangle). Deciding membership by cell centers rather than by nodes keeps the area unbiased; a node-based mask loses about half a cell all around the edge and raises every frequency by several percent.
52**Bending.** We evaluate $w_{xx}$ and $w_{yy}$ by three-point second differences at the nodes and $w_{xy}$ by the four-point difference at cell centers, and integrate the bending energy with the node weights, the bending stiffness $D_k$ at the nodes and $D_c$ at the cells, where the cell thickness is the mean of its four corners:
54$$
55U_b = \frac{h_x h_y}{2}\left[\sum_k \omega_k D_k \left(a_k^2 + b_k^2 + 2\nu a_k b_k\right) + 2(1-\nu)\sum_c D_c\, c_c^2\right],
56$$
58with $a_k$, $b_k$ the second differences at node $k$ (absent where undefined) and $c_c$ the mixed difference at cell $c$. Where $D$ is constant the interior stencil is the standard 13-point biharmonic; with a varying $D$ the same gradient produces the discrete form of the variable-stiffness operator, and the edges are handled as the energy dictates. The three rigid modes $w = 1, x, y$ lie in the null space for any thickness field and any outline.
60**Membrane.** The cells are bilinear (Q4) finite elements with 2 by 2 Gauss integration, using the strains of Section 1 with $w_0$ and $w$ interpolated bilinearly and the stiffnesses of the cell's mean thickness. One-point integration would be cheaper but admits hourglass modes, which the nonlinear coupling would excite.
62**Mass.** The mass matrix is lumped: the mass of node $k$ is $m_k = \rho h_k h_x h_y \omega_k$, so the total mass is the exact integral of $\rho h$ over the staircase plate.
64**Accuracy.** Two errors limit the accuracy of the spatial scheme. The first is the treatment of the edge. For the rectangle the scheme converges at second order: the first seven frequencies of a completely free square plate are within 8 percent of Gorman's reference values at 10 cells per side and within 2 percent at 20. The staircase edge of a curved outline is only first-order accurate; the first nine frequencies of a free circular plate are within 3 percent of Leissa's values at 40 cells across the diameter, and the degenerate pairs, which the staircase splits, differ by about 1 percent. A 16 inch cymbal at 8 mm has about 50 cells across, so its low modes are placed to within a percent or two, and the splitting of the pairs is of the same order as the splitting that the hammering of a real cymbal produces. The second error is grid dispersion, which affects every partial regardless of the outline. A plane bending wave of wavenumber $k$ along a grid axis has the discrete frequency ratio $(2 - 2\cos kh)/(kh)^2$, so a partial comes out low by 3 percent at $kh \approx 0.60$ and by 10 percent at $kh \approx 1.12$. Since the bending frequency at a given $k$ is proportional to the thickness, these points move up with the thickness: for 1 mm bronze on the 8 mm grid the 3 percent point is near 1.1 kHz and the 10 percent point near 3.6 kHz. The interface reports the two frequencies at the mean thickness of the design. The high partials are therefore audibly flattened at 8 mm, in both models alike, and a finer grid is a matter of patience rather than of the method.
66## 3. The stand
68A cymbal sits on a stand: a rod through the center hole, with a felt pad below the cymbal and usually another above it under a wing nut. We model each pad as an annulus of felt of inner radius $r_1$, outer radius $r_2$ and thickness $t$, with Young's modulus $E_p$ and loss factor $\eta$. Its stiffness per unit area in compression is $E_p/t$, and its stiffness per unit area in shear is $G_p/t$ with $G_p = E_p/(2(1+\nu_p))$ and $\nu_p = 0.3$; the shear is what holds the cymbal in place sideways, and it stands in for the plastic sleeve of the rod, which is not modeled (Section 7). The stand energy is
70$$
71U_s = \frac12 \sum_{k \in \text{pads}} \left[ k_k\, w_k^2 + g_k (u_k^2 + v_k^2) \right],
72$$
351af76Interactive tutorial on the physics and the mathematics of the simulatorJeremy Magland 74where 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 script.
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.
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.
80## 4. The linear model
82**Modes.** Collecting the linear part of Section 1 on the grid of Section 2 gives the matrix equation $M \ddot q + K q = f$, with $q$ the vector of unknowns (the transverse displacement of every plate node for a flat plate, and $u$, $v$, $w$ interleaved for a plate with a profile), $M$ the diagonal lumped mass, and $K$ the symmetric stiffness, the Hessian of $U_b + U_m + U_s$ at the rest state. The tests check that $K q$ equals the linearized force of the full model to rounding. The modes are the solutions of $K\phi = \omega^2 M \phi$, normalized so that $\phi^T M \phi = 1$; a free plate has three (flat) or six (with a profile) rigid modes at $\omega = 0$, which are computed and dropped, since a hanging plate's rigid motion is taken up by the thread and does not radiate in the model of Section 4 below. A plate on a stand has none.
84The stiffness is banded, with a half-bandwidth of about twice the number of nodes across the shorter side of the grid (times three with a profile), and its $LDL^T$ factorization at a shift, $K - \sigma M = LDL^T$, costs little. We use it in two ways. The number of negative pivots of $D$ is the number of modes with $\omega^2 < \sigma$ (Sylvester's law of inertia), which tells exactly how many modes lie in any interval. And the factorization applies $(K - \sigma M)^{-1}$, which turns the modes near $\sigma$ into the dominant eigenvalues of a symmetric operator that the Lanczos method finds in a few dozen steps. All modes up to a cutoff frequency (16 kHz by default) are computed by *spectrum slicing*: the interval is cut into slices of about 30 modes (100 with a profile) at boundaries equally spaced in frequency, refined where a slice is too full, the inertia at each boundary fixes the count in each slice, and each slice is solved by shift-invert Lanczos with full reorthogonalization from a random start at the slice's center, until the count of converged modes in the slice equals the count from the inertias; a slice that comes up short (a degenerate pair of a disc can be missed from an unlucky start) is restarted with a new random vector orthogonal to the modes found. The slices are independent and run in a pool of Web Workers. A 16 inch cymbal at 8 mm has about 2150 nodes and about 1900 bending modes below 16 kHz; they take about 2.5 s on a laptop with 12 threads. With a profile the same cymbal has 6500 unknowns and a band three times as wide, and the factorizations and solves become bound by memory traffic; the modes then take about 30 s on the same laptop.
86**Modal equations.** Expanding $q = \sum_m \phi_m q_m(t)$ turns the plate into independent oscillators,
88$$
89\ddot q_m + 2\sigma_m \dot q_m + \omega_m^2 q_m = \phi_m^T f(t),
90$$
92with $\sigma_m = \ln(1000)/T_{60}(f_m)$ from the design's decay times, plus the pad term of Section 3 on a stand. The listener signal and the movie are then sums over the modes.
02e2d5cDefault contact stiffness 1e7 N/m^1.5Jeremy Magland 94**Strike.** The stick tip is a free mass $M_s$ (30 g by default, the effective mass of a stick in a wrist stroke) arriving with speed $v_0$ at a point of the plate. The contact patch consists of the plate nodes within the tip radius (5 mm) plus the nearest one; its contact coordinate is the mass-weighted mean transverse displacement $w_c = \sum_k (m_k/M_p) w_k$ of these nodes, the penetration is $\delta = w_c - z$ with $z$ the tip position, and the force is the Hertz law $F = K\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ and $K = 10^7\ \mathrm{N/m^{1.5}}$ by default, a softer contact than the Hertz value for a hard wooden tip on metal (of order $10^8$ to $10^9$), which stands for the compliance of a real tip and gives a rounder attack; the value is a parameter of the strike. The force is spread over the patch in proportion to the node masses, which makes the pair conjugate: the modal force is $-F g_m$ and the contact coordinate is $\sum_m g_m q_m$ with the same gains $g_m = \sum_k (m_k/M_p)\, \phi_m(k)$. The stick obeys $M_s \ddot z = F$. The contact is stepped at 2 microseconds; over each step the force is held constant, and every oscillator is advanced by its exact solution under a constant force, so the step is accurate and unconditionally stable for the plate, and only the stick is integrated explicitly. Within the hit the stick may chatter, releasing and re-engaging as the plate runs ahead of it and comes back; once the contact has been open for 5 ms the stick is withdrawn, so that a plate rocking on its stand cannot strike it again tens of milliseconds later, as a drummer lifts the stick well before that. After that instant $t_1$ every mode decays freely, and its state is one complex amplitude, $q_m(t) = \mathrm{Re}\left[A_m e^{s_m (t - t_1)}\right]$ with $s_m = -\sigma_m + i\sqrt{\omega_m^2 - \sigma_m^2}$. A thin plate is a soft target: its driving-point impedance is $8\sqrt{D\rho h}$, about 60 N s/m for 1 mm bronze, so the contact lasts about a millisecond and the peak force is a few hundred newtons; the excitation is broadband, and the sound is a short bright attack followed by the ring. A hard hit in the linear model is simply a loud one; the model does not know that the plate would have stiffened.
96**Radiation.** The listener model is the Rayleigh integral for a baffled plate,
98$$
99p(\mathbf r, t) = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \ddot w_k(t - \tau_k), \qquad \tau_k = r_k / c,
100$$
102where the sum runs over the plate nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s; the common delay is discarded and only the spread of the delays across the plate is kept, up to about a millisecond for a 16 inch cymbal. For a mode in free decay the delay is a phase: $\ddot w_k(t - \tau_k) = \mathrm{Re}[\phi_m(k)\, s_m^2 A_m e^{s_m(t - t_1)} e^{-s_m \tau_k}]$, so the whole sum over the plate collapses to one complex gain per mode,
104$$
105G_m = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \phi_m(k)\, e^{-s_m \tau_k}, \qquad
106p(t) = \sum_m \mathrm{Re}\left[G_m s_m^2 A_m e^{s_m (t - t_1)}\right],
107$$
109which is exact for the free decay and includes the directivity of every mode. The signal is generated by the two-term recurrence of each damped sinusoid at the audio rate, about 1900 recurrences per sample. Before the free decay, while the modes are being driven by the stick, the delayed sum cannot be collapsed, and it is evaluated directly: the modal accelerations are tabulated at twice the audio rate during the contact phase, the nodal accelerations $\ddot w_k(t) = \sum_m \phi_m(k) \ddot q_m(t)$ are formed from them, and the Rayleigh sum is taken with linear interpolation in the table. The table path covers the output up to the moment at which every node's retarded time is past the last contact, about 2 ms after the impact, and the modal path takes over from there; the two join without a seam, and the tests check that they do. The linear interpolation in the table during the first two milliseconds is the one approximation of the linear path beyond the modal truncation itself. Moving the listener changes only the gains, so the sound follows a move in well under a second.
111**Movie.** The transverse displacement at the frame times is the modal sum $\sum_m \phi_m(k) q_m(t)$, with $q_m$ from the contact tables before $t_1$ and from the free decay after it. Its cost, 1900 modes times 2150 nodes times 480 frames, is the largest part of the synthesis and is split among the workers, after the audio, so that playback can start first.
113## 5. The full model
115The 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.
117**Losses.** The two-parameter model gives a bending mode of angular frequency $\omega$ and wavenumber $k$ the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes the high frequencies decay faster in proportion to frequency. The two parameters are fitted by non-negative least squares to the design's rates $\ln(1000)/T_{60}(f)$ on log-spaced frequencies from 50 Hz to 10 kHz, with $\kappa$ at the mean thickness; a decay curve that is not of this form (a table with a kink, for instance) is matched only approximately, and the interface reports the fitted decay times at 100 Hz, 1 kHz and 5 kHz next to the requested ones. The $\sigma_1$ term is implemented as $2\sigma_1$ times the graph Laplacian of the transverse velocity and does not depend on the thickness; where the plate is thicker, a given frequency corresponds to a smaller $k^2$, so the frequency-dependent loss is weaker there. The in-plane motion carries an additional frequency-independent rate, taken as the design's rate at 5 kHz: the stretching modes lie at several kilohertz, couple linearly to the transverse motion once the plate is curved, and would otherwise ring for the whole low-frequency decay time.
119**Free plate.** A sheet hanging from a thread is treated as free, and the reaction to the stick is distributed over the plate as the inertial load of a rigid body, which is the inertia-relief equilibrium of a free body: the force of the previous step is applied as the rigid acceleration field $a_k = F\,(1/M + y_k y_c/I_x + x_k x_c/I_y)$ that balances the patch load and its moments about the center, so the plate gains no net momentum or angular momentum. The one-step lag has zero net impulse. A plate on a stand takes the reaction through the pads. This is the counterpart of dropping the rigid modes in the linear model.
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 121**Single precision.** Three places in the kernel would otherwise turn rounding into audible noise. The first is the membrane strain of a curved design. Written directly, $\varepsilon_{xx} = u_x + \tfrac12[(w_0+w)_x^2 - w_{0,x}^2]$ subtracts two nearly equal numbers once $w$ is small compared with the profile, and the rounding error of the two squares is proportional to $w_{0,x}^2$, a constant: it does not fall as the vibration decays, and it acts as a broadband forcing that leaves a noise floor which never dies away. The strains are therefore written so that every term carries the deviation, $\varepsilon_{xx} = u_x + w_x(w_{0,x} + \tfrac12 w_x)$ and $\gamma_{xy} = u_y + v_x + w_{0,x}w_y + w_{0,y}w_x + w_x w_y$ with $w$ now the deviation alone, which is algebraically the same and has nothing large to cancel. On the 16 inch crash with a bell, whose largest profile slope is about 0.23, the direct form left the top octave at a constant level for the whole run, about 25 dB below the peak of a light hit; the rearranged form puts it more than 60 dB down and decaying. The second is the update of $w$. With $w$ of a few millimeters and $\Delta t \approx 2\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would again act as a broadband noise source. The update of $w$ is therefore compensated: the rounding error of each addition is carried in a separate array and added back at the next step. The third is the listener sum, which forms the acceleration as the difference of successive delayed velocities at the audio rate, the average acceleration over one sample interval, rather than sampling the instantaneous acceleration, whose per-step rounding error of order $1\ \mathrm{m/s^2}$ per node would go straight into the audio. What remains after all three is a noise floor proportional to the signal, roughly 75 dB below the peak in the top octave for a strongly curved design and 100 dB for a flat one.
123**The kernel.** The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch in which each workgroup of 256 threads advances a 12 by 12 tile of nodes: the threads load the tile with a two-node halo into shared memory, compute the membrane forces of the cells and the bending moments at the nodes once each, and then evaluate the update of the tile's nodes from shared memory. On grids of a few thousand nodes the cost of a step is the latency of one such workgroup rather than the arithmetic. A per-step uniform, addressed by a dynamic offset, carries the losses and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass; the retarded-time listener sum runs once per frame for all of the frame's samples from a ring buffer of the transverse velocities, and the following frame is encoded while the GPU works on the current one. A 16 inch cymbal at 8 mm runs at about 10 s per second of sound on the integrated GPU of a laptop.
125## 6. Validation
127The tests check, without a GPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the full model are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered, free and on a stand; that the loss terms and the pad dashpots only remove energy; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate on a stand; that the banded factorization solves and that its inertia counts the eigenvalues; that every mode of a small square and of a small disc, degenerate pairs included, matches a dense eigensolution to $10^{-7}$ with no mode missed; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that the stand resists a translation with its total stiffness and a tilt with the pads' second moment, and that the lift-off rule unloads the right pad; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, that the attack and the tail of the synthesized signal join, and that the parts computed by the workers sum to the whole.
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 129With 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.
131## 7. Limitations
133The linear model is linear: it has no pitch glide, no wash, no lift-off, and a hard hit is only a loud one. It is the right tool for the tone, the modes and the decay of a design, and for a first listen, and the full model is the right tool for the sound of a hard hit; the two agree for a light one. The modal truncation at the cutoff frequency removes the highest partials, which the grid places poorly anyway; the linear interpolation of the acceleration table during the first two milliseconds slightly blurs the attack.
135The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a strike that moves a 400 mm cymbal by a few millimeters is well within its range. The Marguerre treatment of the profile is a shallow-shell approximation, adequate for a bow of a few millimeters on a sheet and for a gentle cymbal bow, but strained by the bell of a real cymbal, whose slope is several times the range in which the approximation holds; the curved preset uses a flatter bell than a real cymbal, and the interface reports the largest slope of a design. A deep bell also behaves as a stiff cup that shifts the bell modes up in ways the shallow shell does not capture. The thickness field is treated as a slowly varying coefficient in a thin-plate model, which is appropriate for a lathed taper but not for a step or a groove on the scale of the grid. The hammering of a real cymbal, which leaves residual stresses and a slightly non-axisymmetric shape, is absent unless the script adds it. The loss model is empirical and its parameters are not calibrated against a real cymbal, and the full model matches the design's decay curve only as well as the two-parameter fit allows. Grid dispersion lowers the high partials by the amounts given in Section 2, and the staircase edge places the low frequencies within a few percent and splits modes that would be degenerate on a true disc.
137The 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.
↑↓ move↵openescclose