Model and numerical method#

This 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, hanging freely, struck once with a drumstick. The plate is described by a script; the strike and the listener are parameters of the setup. 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.

Section 1 gives the plate model, Section 2 its discretization on the grid, Section 3 the linear model and the modal synthesis, Section 3 the full model, Section 4 the checks that were made, and Section 5 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.

1. Plate model#

The 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

so 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

where $\mathcal{B}[w]$ is the variable-stiffness bending operator, the Euler-Lagrange derivative of the bending energy

which 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.

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.

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.

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.

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 6).

Losses. The decay of the sound is not prescribed. A partial of angular frequency $\omega$ decays as $e^{-\sigma t}$, and its decay time $T_{60} = \ln(1000)/\sigma$ is the time in which it falls by 60 dB; the rate $\sigma$ is the sum of three mechanisms, each with its own physics and its own constants in the script. Internal friction of the metal is a hysteretic loss: each cycle dissipates the fraction $2\pi\eta$ of the stored energy, so a mode decays at $\sigma = \eta\omega/2$, and the loss factor $\eta$ is a material constant on the same footing as $E$ and $\rho$, measured rather than derived, of order $10^{-4}$ for bronze, steel and aluminum. Thermoelastic damping comes from the heat that flows across the thickness between the compressed and the stretched face of a bent plate; Zener's analysis gives the loss factor $\eta_{te}(\omega, h) = \Delta\,\omega\tau/(1 + \omega^2\tau^2)$ with the relaxation strength $\Delta = [E\alpha^2 T/(\rho c_p)]\,(1+\nu)/(1-\nu)$, the biaxial form of the beam value, and the relaxation time $\tau = h^2\rho c_p/(\pi^2 k)$, the time heat takes to diffuse across the thickness, from the thermal expansion $\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$ (293 K unless the script says otherwise). The loss peaks at $\omega\tau = 1$, below 100 Hz for a millimeter of metal, and is small for the presets; it is what damps thin aluminum. Radiation into the air is the energy the plate loses as sound, which the linear model computes exactly for each mode from its shape (Section 3) and the full model takes from an analytic average (Section 4). It is smaller than one might expect. A plate with free edges radiates poorly below its coincidence frequency, the frequency at which the bending wavelength equals the wavelength of sound (about 17 kHz for 1 mm bronze): its modes are orthogonal to the rigid translation and rotations, so they have no net volume velocity and no dipole moment, and at a free edge the bending moment and the shear vanish, which removes the edge radiation that a simply supported plate has. The efficiency rises steeply toward coincidence, so on the presets the radiation overtakes the internal friction only in the top octaves, above about 5 kHz, and it is there that it makes the high partials die first; below that the decay of a free plate is set by the metal. The linear model adds the three rates for each mode, and the mode browser shows the decay time each would give alone. The full model needs a loss operator rather than a rate per mode and uses the standard two-parameter model of plate synthesis, a rate $\sigma_0 + \sigma_1 k^2$ for a bending wave of wavenumber $k$, fitted to the estimated curve (Section 4). The loss factors and thermal constants of the presets are typical published values, not measurements of a particular cymbal.

2. Spatial discretization#

The 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.

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.

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:

with $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.

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.

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.

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.

3. The linear model#

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 3 below. A plate on a stand has none.

The 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.

Modal equations. Expanding $q = \sum_m \phi_m q_m(t)$ turns the plate into independent oscillators,

with $\sigma_m$ the sum of the three loss rates of the mode. The internal friction gives $\eta U_m/\omega_m$ with $U_m$ the strain energy of the mode, $\omega_m^2/2$ for a mass-normalized mode, that is $\eta\omega_m/2$; the thermoelastic loss gives $\sum_k \eta_{te}(\omega_m, h_k)\,U_{b,k}/\omega_m$, the local loss factor weighted by the bending energy $U_{b,k}$ of the mode at each node (the same discrete bending energy as Section 2, with the cell terms shared among their corners), which reduces to $\eta_{te}\omega_m/2$ on a flat uniform plate. The radiation rate is the power the mode radiates over twice its energy. A baffled plate whose nodes move with the velocities $v_k$ radiates the power $v^T R v$ from one face, with the radiation resistance matrix $R_{jk} = (\rho_{\mathrm{air}}\omega^2/4\pi c)\,\Delta A_j \Delta A_k\, \mathrm{sinc}(k_a r_{jk})$ of Elliott and Johnson, $k_a = \omega/c$; both faces radiate, the velocity of a mass-normalized mode of amplitude $A$ is $\omega A\phi$ and its energy $\omega^2 A^2/2$, so

which for a small piston at low frequency reduces to the textbook resistance $\rho_{\mathrm{air}}\omega^2 S^2/(2\pi c)$ divided by the mass. On a regular grid the kernel depends only on the offset $j - k$, so the double sum is a convolution, and by Parseval it is a sum over the spectrum: two two-dimensional FFTs of the zero-padded grid per mode, about 5 million operations each for a 16 inch cymbal at 8 mm, split among the workers, rather than $N^2$ terms. The treatment of a node's area as a point source holds while $k_a h$ is small; near coincidence the mode shape itself is barely resolved by the grid, so nothing is lost there. The listener signal and the movie are then sums over the modes.

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.

Radiation. The listener model is the Rayleigh integral for a baffled plate,

where 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,

which 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.

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.

4. The full model#

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.

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 3 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.

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. This is the counterpart of dropping the rigid modes in the linear model.

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.

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.

5. Validation#

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; that the loss terms only remove energy; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate; 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 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.

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, 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 flat and on a tapered plate. 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 4. 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.

6. Limitations#

The linear model is linear: it has no pitch glide, no wash, 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.

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.

The plate is free: nothing supports it, and neither the thread a sheet hangs from nor the stand a cymbal sits on is modeled, so the losses and the low modes a support adds are absent. 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.