Model and numerical method#

This note describes the model behind the simulator: a thin stainless steel sheet, rectangular or elliptical, hanging from a thread, that is either pushed at a point of its rim, held, and released, or struck with a drumstick. The sheet rings like a gong. Its vibration amplitude, several thicknesses, is beyond the range where a linear plate model is accurate, so we solve the geometrically nonlinear plate equations of Foppl and von Karman, keep the in-plane inertia, discretize on a regular grid, and step explicitly on the GPU. Section 1 gives the continuous model, Sections 2 and 3 the suspension and the two excitations, Sections 4 and 5 the discretization in space and time, Section 6 the radiation model, and Section 7 the limitations.

1. Plate model#

The sheet occupies the rectangle $[-L_x/2, L_x/2] \times [-L_y/2, L_y/2]$, or the ellipse inscribed in it, and has thickness $h$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$. 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 initial shape $w_0(x,y)$ that represents the bow of the real sheet (below). 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, the equations of motion are

with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. Bending acts on the deviation $w$ from the initial shape; the membrane forces act through the total slope. For $w_0 = 0$ these are the von Karman equations. The quadratic terms in the strains couple bending to stretching and produce the gong-like behavior: the pitch glide of a strongly excited plate and the transfer of energy from the low modes toward high frequencies. The in-plane inertia terms $\rho h u_{tt}$ and $\rho h v_{tt}$ 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 a time step set by the in-plane wave speed (Section 5).

Initial shape. Real sheet metal is never flat; rolling and handling leave a bow of a few millimeters. A bowed sheet is a shallow shell, and its behavior differs qualitatively from a flat plate: the flat plate has a symmetric, purely hardening nonlinearity, whereas a shell has an asymmetric one, may snap between two curved states, and oscillates more slowly at large amplitude as it passes through the soft flat configuration. This is the mechanism behind oil-canning, and it appears to be what produces the slow wobble of a real sheet. Three parabolic shapes are offered: curved across the width, curved across the height, and a dish with the same curvature both ways. The sag is prescribed at the center of the bounding box; for the two cylindrical shapes it is measured relative to the straight edges, and for the dish relative to the corners, with the edge midpoints at half the sag. The default is a 3 mm bow across the width. A positive push moves the rim point in the direction of the bow and a negative push against it; the two cases sound different, as they do on a real sheet.

Losses. The two-parameter loss model is the standard one in plate sound synthesis. A bending mode of angular frequency $\omega$ and wavenumber $k$ decays at the rate $\sigma_0 + \sigma_1 k^2$, and 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 high frequencies decay faster in proportion to frequency. The defaults, $\sigma_0 = 1.15\ \mathrm{s^{-1}}$ and $\sigma_1 = 5\times10^{-4}\ \mathrm{m^2/s}$, give 60 dB decay times of about 4.9 s at 100 Hz, 1.8 s at 1 kHz, and 0.5 s at 5 kHz. The in-plane motion carries an additional frequency-independent loss $\sigma_{\mathrm{in}}$ (default $20\ \mathrm{s^{-1}}$). The stretching modes of the sheet lie at several kilohertz, couple linearly to the transverse motion once the sheet is curved, and would otherwise ring for the whole low-frequency decay time; a material loss factor of about $10^{-3}$ at 5 kHz gives a rate of this order. None of these values are measured; they are chosen to sound plausible, and all three are exposed in the interface.

2. Suspension and rigid-body motion#

The whole edge is free: the bending moment and the Kirchhoff shear vanish, and so does the in-plane traction. Rather than imposing these conditions explicitly, the discretization of Section 4 is derived from a discrete energy, and free edges are then its natural boundary conditions.

The thread through the hole near one corner is not modeled as a constraint. At audio frequencies it is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the average position of the plate. We therefore treat the plate as free and remove the six rigid-body modes (three translations, three rotations) from the displacement and velocity fields by a mass-weighted projection every few time steps. Gravity is dropped: the in-plane stress it induces is about $\rho g L \approx 35$ kPa, a strain of $2\times10^{-7}$.

The same projection defines what holding the sheet means during a push. The finger prescribes the displacement of a small patch, and the reaction to the finger force is distributed over the plate as an inertial load, which is the inertia-relief equilibrium of a free body, so the push produces pure bending of the prescribed amount. This is a modeling choice: a sheet hanging freely from a thread and pushed 10 mm at a corner would mostly swing as a pendulum and bend by only a millimeter or two. The same choice applies to the drumstick, whose force would otherwise set the whole sheet swinging.

3. Excitation#

Push and release. The simulation has three phases. During the push (default 0.1 s) the transverse displacement of the nodes within the finger radius follows a smooth ramp to the push amount, and heavy artificial damping ($\sigma_0 = 300\ \mathrm{s^{-1}}$) removes the waves that the ramp launches. During the hold the constraint stays in place until the non-rigid kinetic energy drops below a threshold or a maximum settle time elapses. At the release the clock is set to zero, the damping returns to its physical value, and the recording of audio and movie frames begins. The finger does not let go instantaneously: the reaction force that held each node of the patch is recorded at the end of the hold and then applied as an external force that falls smoothly to zero over the release time (20 ms by default). This matters for the sound. An instantaneous release is a step in the finger force, which excites all modes at once, and since the radiated pressure is proportional to the plate acceleration, the result is a nearly flat pressure spectrum, a cymbal-like click. A release over 20 ms confines the excitation to roughly the first hundred hertz, as a fingertip does. The rigid-body components of the recorded reaction are removed before it is applied, so that the finger force and the distributed holding reaction relax together. The release time can be set to zero to recover the instantaneous case.

Drumstick strike. The stick tip is a free mass $M$ (default 30 g, the effective mass of a stick in a wrist stroke rather than its full mass) arriving with speed $v_0$ (default 3 m/s) at a point of the sheet. The contact patch consists of the sheet nodes within the tip radius (5 mm) plus the nearest node. Its contact coordinate is the mass-weighted mean transverse displacement $w_c$ of these nodes, the penetration is $\delta = w_c - x$ with $x$ the tip position, and the force is the Hertz law $F = K\,\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ for a rounded tip and $K = 5\times10^{8}\ \mathrm{N/m^{1.5}}$ for wood on steel. The force is spread over the patch in proportion to the node masses. Because that distribution is exactly conjugate to the mean displacement, the contact has a potential energy and the coupled system is Hamiltonian apart from the losses. The stick obeys $M\ddot x = F$ and is stepped together with the plate; it may bounce and strike again if the sheet comes back to meet it. The recording starts at the impact, with the sheet at rest in its natural shape. Note that a thin sheet is a soft target: its driving-point impedance is $8\sqrt{D\rho h} \approx 59\ \mathrm{N\,s/m}$, so the contact lasts well under a millisecond and the peak force is of the order of a hundred newtons. The excitation is therefore broadband, unlike the push, and the sound is a short bright attack followed by the ring.

4. Spatial discretization#

The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $h_x$ and $h_y$ (8 mm by default). All three displacements live at the nodes. 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.

Bending. The continuous bending energy is $\tfrac{D}{2}\int [w_{xx}^2 + w_{yy}^2 + 2\nu w_{xx}w_{yy} + 2(1-\nu) w_{xy}^2]\,dA$. We evaluate $w_{xx}$ and $w_{yy}$ by three-point second differences at the nodes where both neighbors exist, and $w_{xy}$ by the four-point difference at cell centers, with trapezoidal quadrature weights (1 inside, 1/2 on an edge, 1/4 at a corner). Where a second difference does not exist, at the edges perpendicular to the direction of differencing, its term is simply absent from the energy. In the interior the resulting force is the standard 13-point biharmonic stencil; near the edges the stencil is modified exactly as the energy dictates. The three rigid modes $w = 1, x, y$ lie in the null space.

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. One-point integration would be cheaper but admits hourglass modes, which the nonlinear coupling would excite. On the GPU the force of each cell on its four nodes is computed once per cell within a workgroup and gathered by the nodes.

Mass. The mass matrix is lumped with the same trapezoidal weights, so that the total mass is exact.

Elliptical sheets. The same grid serves an elliptical sheet through a mask. A cell belongs to the sheet when its center lies inside the ellipse, and a node when it is a corner of such a cell; the edge is therefore a staircase of whole cells. Every rule above carries over unchanged: a second difference exists where the node and both of its neighbors belong to the sheet, a cell contributes membrane energy and the mixed bending term when its four corners do, and the weight of a node is the fraction of its four cells that belong to the sheet. Because the forces are still the gradient of the discrete energy, the staircase edge is free, the stiffness stays symmetric, and exactly three rigid modes remain in each block. Nodes outside the sheet are frozen and carry no mass. 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. The hole and the pushed patch are placed on the rim in the upper-left and lower-right directions, so the push remains asymmetric as it is at the corner of a rectangle.

Accuracy. Two errors limit the accuracy of the spatial scheme, and it is worth separating them. 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 low by about 3.5 percent at 10 cells per side, 0.8 percent at 20, and 0.3 percent at 30, compared with the reference values of Gorman. The staircase edge of the ellipse is only first-order accurate; the frequencies of a free circular plate are within about 2 percent of Leissa's values at 24 cells across the diameter and within 1 percent at 60, and the degenerate pairs, which the staircase splits, differ by 1 to 2 percent. 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 1 percent at $kh \approx 0.35$, by 3 percent at $kh \approx 0.60$, and by 10 percent at $kh \approx 1.12$. For the default sheet this places the 3 percent point near 1.1 kHz on the 8 mm grid, 1.9 kHz on a 6 mm grid, and 4.2 kHz on a 4 mm grid, and the 10 percent point near 3.6, 6.5, and 15 kHz respectively. The high partials of a gong sound are therefore audibly flattened at 8 mm, and the interface reports these two frequencies for the chosen spacing. Since the cost scales with the inverse cube of the spacing (Section 5), the default is a compromise in favor of interactive use, and a finer grid is a matter of patience rather than of the method.

Validation. The tests check that both stiffness blocks are symmetric with exactly three rigid modes, that the forces are the exact gradient of the discrete energy at a strongly nonlinear state with and without a bow, that the loss terms only remove energy, that the free-plate frequencies converge as stated, and that the same holds on an elliptical mask. A further test runs one time step of the GPU kernel on a smooth random state of about 2.5 thicknesses amplitude, flat, bowed, and elliptical, and compares it with the float64 reference: the accelerations agree to about $10^{-6}$ relative, which is single-precision rounding.

5. Time integration#

The system is advanced 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 (Stormer-Verlet) 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. We estimate $\omega_{\max}$ at start-up by power iteration on a small grid with the same spacing and use 85 percent of the limit. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)} \approx 5200$ m/s, set this limit: at 8 mm spacing $\Delta t \approx 1.3\ \mu$s, roughly six times smaller than bending alone would require, and 4 s of sound take about 3 million steps; at 4 mm the step halves and the node count quadruples. The step is rounded down so that an integer number of steps (16 at 8 mm and 48 kHz) makes one audio sample.

Single precision needs one precaution. With $w$ of a few millimeters and $\Delta t \approx 1\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would 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 (Kahan summation). The self-test of Section 4 also checks this update.

The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch of 16 by 16 workgroups; a per-step uniform, addressed by a dynamic offset, carries the prescribed patch displacement, the damping, and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass.

6. Radiation to the listener#

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

where the sum runs over the nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s. The transverse velocity of every node is recorded in a ring buffer at the audio rate, and each output sample is a delayed, weighted sum over the plate, evaluated on the GPU with linear interpolation between samples; the acceleration is formed as the difference between successive delayed velocities, which is the average acceleration over one sample interval. This matters in single precision: the force evaluation carries a broadband rounding error of the order of $1\ \mathrm{m/s^2}$ per node, and sampling the instantaneous acceleration would put that noise straight into the audio, whereas the averaged form suppresses it by the square root of the number of steps per sample. The common delay is discarded; only the spread of delays across the plate (about 5 samples for a listener 1 m in front) is kept.

The plate is unbaffled in reality, so the sound from the 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. The pressure is reported in pascals and the peak sound pressure level is displayed. Playback and the WAV file use, by default, a fixed mapping from pressure to full scale (1 Pa at 0 dBFS, shifted by a user gain in decibels), so that runs with different excitations can be compared by ear; a soft clipper bounds the output and the resulting peak level is shown. Two normalized modes are also offered, to the peak or to the loudness of the sustained sound.

7. Limitations#

The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The Marguerre treatment of the initial shape is a shallow-shell approximation, adequate for a bow of a few millimeters, and the real shape of a given sheet is of course unknown; the three offered shapes are idealizations. The loss model is empirical and its parameters are not calibrated against a real sheet. Grid dispersion lowers the high partials by the amounts given in Section 4, and the staircase edge of an elliptical sheet places the low frequencies within a few percent and splits modes that would be degenerate on a true ellipse. The suspension is idealized as inertia relief, so the swinging of the sheet on its thread is absent. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included, and the radiation model is approximate, as discussed above. Finally, single precision limits the dynamic range of the tail of the sound to roughly 100 dB below the peak.