Model and numerical method#

This note describes the model behind cymsim: 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. It follows the note of the browser application cymbalsim, of which cymsim is a port; the differences of implementation are noted where they occur (Sections 3 and 4). The plate is described by a Design; 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. The full model steps the geometrically nonlinear plate equations of Foppl and von Karman explicitly in time; it 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.

1. Plate model#

The plate occupies a region $\Omega$ of the plane, given by the design as a bounding box and an outline function. 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. 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$. For a flat plate 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. A profile couples the two through its slopes: a curved plate stretches when it bends, which raises the low modes.

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

Profile. Only the slopes of the stress-free shape 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; Plate.profile_stats reports the largest slope.

Losses. 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. 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 of order $10^{-4}$. 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)$ and the relaxation time $\tau = h^2\rho c_p/(\pi^2 k)$, from the thermal expansion $\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$. 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). A plate with free edges radiates poorly below its coincidence frequency (about 17 kHz for 1 mm bronze): its modes have no net volume velocity and no dipole moment, and at a free edge the bending moment and the shear vanish. On the presets the radiation overtakes the internal friction only above about 5 kHz.

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$, 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 a staircase of whole cells, and a hole is a region whose cells are absent. 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. Deciding membership by cell centers keeps the area unbiased.

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 (from the mean thickness of the four corners):

Where $D$ is constant the interior stencil is the standard 13-point biharmonic. 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$.

Accuracy. 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, and the degenerate pairs, which the staircase splits, differ by about 1 percent. Grid dispersion 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$; for 1 mm bronze on the 8 mm grid these points are near 1.1 kHz and 3.6 kHz. Plate.dispersion_frequency reports them. 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 on the grid gives $M \ddot q + K q = f$, with $q$ the vector of unknowns (the transverse displacement of every plate node for a flat plate; $u$, $v$, $w$ interleaved for a plate with a profile), $M$ the diagonal lumped mass, and $K$ the symmetric stiffness, the Hessian of $U$ at the rest state; it is assembled as a sparse matrix from the stencil operators of Section 2, and the tests check that $Kq$ 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$; the rigid modes at $\omega = 0$ are dropped, since a hanging plate's rigid motion is taken up by the thread and does not radiate in the model below.

Here cymsim departs from the browser application, which slices the spectrum and solves each slice by shift-invert Lanczos in a pool of workers. With LAPACK at hand we form $A = M^{-1/2} K M^{-1/2}$ as a dense matrix and call the driver that computes the eigenpairs in a given interval, $(2\pi \cdot 1\,\mathrm{Hz})^2 < \lambda \le (2\pi f_{\max})^2$. A 16 inch cymbal at 8 mm has about 2150 nodes and 1900 bending modes below 16 kHz, computed in about 3 s on a 12-thread CPU or 0.3 s on an RTX 6000 Ada with torch.linalg.eigh; with a profile the same cymbal has 6500 unknowns and takes 25 to 40 s on the CPU, or 1.6 s on the GPU. The dense approach stores $n^2$ numbers and its cost grows as $n^3$, so it is the right tool up to roughly 10 000 unknowns; beyond that, spectrum slicing would be the remedy.

Modal equations. Expanding $q = \sum_m \phi_m q_m(t)$ turns the plate into independent oscillators $\ddot q_m + 2\sigma_m \dot q_m + \omega_m^2 q_m = \phi_m^T f(t)$, with $\sigma_m$ the sum of the three loss rates of the mode. The internal friction gives $\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, 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, 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 2-D FFTs of the zero-padded grid per mode, batched over the modes, rather than $N^2$ terms.

Strike. The stick tip is a free mass $M_s$ (30 g by default) 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, which stands for the compliance of a real tip. 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 unconditionally stable for the plate and only the stick is integrated explicitly. Once the contact has been open for 5 ms the stick is withdrawn. 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}$.

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

with $c = 343$ m/s; the common delay is discarded and only the spread of the delays across the plate is kept. For a mode in free decay the delay is a phase, 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 in blocks as a complex matrix product of the powers $e^{s_m n/f_s}$ with the coefficients. Before the free decay, while the modes are being driven by the stick, the delayed sum is evaluated directly from a table of the modal accelerations at twice the audio rate, with linear interpolation; 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 there. The tests check that the two join without a seam.

Movie. The transverse displacement at the frame times is the modal sum $\sum_m \phi_m(k) q_m(t)$, one matrix product of the frame-by-mode amplitudes with the shapes.

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 by default. 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. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system, computed at start-up as the largest eigenvalue of the three-field stiffness of 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: at 8 mm spacing $\Delta t \approx 1.9\ \mu$s, 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; we keep them because they make the whole system explicit and local, at the cost of the time step just described.

Losses. The full model needs a loss operator rather than a rate per mode and uses the two-parameter model of plate synthesis: a bending mode of angular frequency $\omega$ and wavenumber $k$ decays at the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$. 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 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), after Leppington, Broadbent and Heron, reduced by the factor $(f/f_c)^2$ that the free edges impose. 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. The in-plane motion carries only the internal friction, as a frequency-independent rate at the lowest in-plane resonance of the plate.

Free plate. The reaction to the stick is distributed over the plate as the inertial load of a rigid 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. 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 that does not fall as the vibration decays and acts as a broadband forcing. The strains are therefore written so that every term carries the deviation, $\varepsilon_{xx} = u_x + w_x(w_{0,x} + \tfrac12 w_x)$ with $w$ the deviation alone. 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$, so the update is compensated, with the rounding error of each addition 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 rather than sampling the instantaneous acceleration. What remains 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; dtype="float64" removes it at a cost in speed on consumer GPUs.

The kernel. Here cymsim departs from the browser application's hand-written WGSL kernel. The step is a pure function of PyTorch tensors over the $(n_y, n_x)$ grid: the second differences and the cell quantities are slices and pads of the whole arrays, the membrane forces of the four Gauss points are accumulated per cell corner and scattered to the nodes by four padded additions, the stick is a handful of scalar tensor operations, and the transverse velocity at the sample instants is written into a ring buffer from which the retarded-time listener sum is formed once per batch of one movie frame. As written the step is some two hundred small tensor operations, and on the CPU or a GPU alike it would be bound by their dispatch; torch.compile fuses them into a few kernels, which makes the step about six times faster on the CPU (0.17 ms for the 16 inch crash, about 90 s per second of sound). On a GPU the fused step is still bound by the Python call and the kernel launches (about 0.1 ms), so the steps of one audio sample are captured into a CUDA graph on static buffers and replayed, one launch per sample; the step then costs about 7 microseconds on an RTX 6000 Ada, 3.9 s per second of sound for the 16 inch crash at 8 mm, and hardly more at four times the nodes. The same code runs on CUDA, on Apple's MPS and on the CPU.

5. Validation#

The test suite checks, on the CPU: 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 kernel 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 kernel to $10^{-6}$ on a bowed, tapered, holed plate; that the three eigensolvers agree and return mass-normalized shapes that satisfy the eigen equation; 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 and the FFT convolution equals the direct double sum for every mode, that simply supported shapes reproduce Wallace's classical efficiencies to a few percent, 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, and that the attack and the tail of the synthesized signal join; that a single-precision step of the kernel from a smooth random state of a few thicknesses amplitude agrees with the double-precision step to single-precision rounding, that the energy is conserved without losses, and that the compiled step equals the uncompiled one; and that a light hit (0.05 m/s) on a small plate simulated both ways gives magnitude spectra below 6 kHz with a correlation above 0.98 and levels within 25 percent, the remaining difference coming from the fitted two-parameter loss against the per-mode decay times and the modal truncation.

6. Limitations#

The linear model is linear: it has no pitch glide, no wash, and a hard hit is only a loud one. The modal truncation at the cutoff frequency removes the highest partials, which the grid places poorly anyway. The von Karman model assumes moderate rotations, and the Marguerre treatment of the profile is a shallow-shell approximation, adequate for a bow of a few millimeters but strained by the bell of a real cymbal, whose slope is several times the range in which the approximation holds. The thickness field is treated as a slowly varying coefficient, appropriate for a lathed taper but not for a step on the scale of the grid. The hammered preset adds the shape of hammering but not its residual stresses. The loss factor of the metal and its thermal constants are typical published values rather than measurements of a particular alloy. The radiation loss is that of a plate in a rigid wall, an overstatement at low frequency for a plate in free air; the losses to the air's viscosity and to the support are absent; and the full model carries the losses only through the two-parameter fit to an average. Grid dispersion lowers the high partials by the amounts given in Section 2. The plate is free: neither a thread nor a stand is modeled, and air loading is not included. Single precision limits the dynamic range of the tail of the full model's sound to roughly 100 dB below the peak.