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 we make every external load self-equilibrated: the reaction to the finger or 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 sheet then gains no net momentum or angular momentum, which stands in for the thread taking them up. 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, and a struck sheet would recoil and swing. Gravity is dropped: the in-plane stress it induces is about $\rho g L \approx 35$ kPa, a strain of $2\times10^{-7}$.

Inertia relief takes two forms here. During a push the finger prescribes the displacement of a small patch, and the reaction is not known in advance, so the six rigid-body modes (three translations, three rotations) are removed from the displacement and velocity fields by a mass-weighted projection every 16 time steps (about 20 microseconds) while the sheet is held; the push therefore produces pure bending of the prescribed amount. The drumstick force is known at every step, so its reaction is applied directly as the rigid acceleration field that balances the force and its moments about the center (Section 3). Once the sheet rings freely, no load remains that could move it as a whole, apart from rounding and the small non-invariance of the von Karman energy under finite rotation, and the projection is optional there and off by default. Note that a rigid rotation is not exactly a null mode of the von Karman energy: the strains contain $\tfrac12 w_x^2$, so a rotation $\theta$ about an in-plane axis produces a spurious strain $\theta^2/2$. This is the reason the loads have to be self-equilibrated, since a slow drift in rotation would otherwise stress the sheet.

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 sheet has settled or a maximum settle time elapses. The sheet counts as settled when its non-rigid kinetic energy drops below a threshold, or when that energy has stopped falling: in single precision the rounding of the force evaluation leaves a floor of a few times $10^{-8}$ J that can lie above any fixed threshold, so the hold also ends once the energy has not halved over the last 50 ms (or over three settle time constants, whichever is longer). With the default settle damping the release comes about 0.14 s into the hold. 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. At the same time the force is balanced over the whole sheet by the inertia-relief field $a_k = F\,(1/M_s + y_k y_c/I_x + x_k x_c/I_y)$, where $M_s$ is the mass of the sheet, $I_x$ and $I_y$ its moments of inertia about the in-plane axes, and $(x_c, y_c)$ the mass-weighted center of the patch, so that the total force and the two moments of the load vanish (Section 2). The field uses the force of the previous step, a lag of one microsecond whose net impulse over the contact is zero. 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. A series of hits is scheduled at fixed intervals with speeds that are constant, rising, or falling from hit to hit. At each scheduled instant the stick is replaced: the tip is put at the current contact coordinate $w_c$ of the moving sheet, wherever that is, with the scheduled speed toward the sheet and zero force, and the rebounding stick of the previous hit is discarded. Each hit therefore starts from touching contact exactly as the first does, and the sheet keeps whatever motion it has from the earlier hits. 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 90 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 (15 at 8 mm and 48 kHz) makes one audio sample.

In-plane mass scaling. Since the in-plane waves set the step and the in-plane inertia is a small effect (the usual Airy stress function formulation drops it altogether), the interface offers a factor $s \geq 1$ by which the in-plane inertia is scaled; the default is 1. The in-plane wave speed becomes $c_p/\sqrt{s}$ and the stable step grows by up to $\sqrt{s}$, until bending sets the limit at about $s = 30$ on the 8 mm grid. The trade-off is that the in-plane modes of the sheet, which start near 8 kHz for the default sheet, move down by the same factor and, once the sheet is curved, are coupled to the transverse motion and to the sound. Note that the scaling moves these modes in the opposite direction from the quasi-static limit, which places them at infinite frequency. A factor of 4 halves the run time and keeps the in-plane modes above 4 kHz; it is meant for quick previews at a finer spacing rather than for the final sound. The kinetic energy reported by the interface includes the factor, so the settle criterion is unaffected.

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 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, about 17 microseconds on an integrated GPU, rather than the arithmetic, so the shared phases, which avoid recomputing the stencil ingredients at every node, roughly halve the step time; the step time barely grows until the grid exceeds several thousand nodes. 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, and the following frame is encoded while the GPU works on the current one. The mass-weighted sums over the sheet (rigid-body components, kinetic energy, peak displacement) run in one workgroup every 16 steps while the sheet is held and once per frame otherwise, or every 32 steps when the optional projection of Section 2 is on.

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 once per movie frame for all the frame's samples, 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.