f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 1# Theory and numerical method
3This note describes the model behind the simulator: a thin stainless steel sheet, suspended by a thread, pushed at a corner, held, and released. The sheet rings like a gong, which means that its vibration amplitude (about ten thicknesses here) is far beyond the range where a linear plate model is accurate. We therefore solve the geometrically nonlinear plate equations of Foppl and von Karman with the in-plane inertia retained, discretize them on a regular grid, and step them explicitly on the GPU.
5## 1. Plate model
7The sheet occupies the rectangle $[-L_x/2, L_x/2] \times [-L_y/2, L_y/2]$ and has thickness $h$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$. Its motion is described by three displacement fields, the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$. With the membrane strains
9$$
10\varepsilon_{xx} = u_x + \tfrac12 w_x^2, \qquad
11\varepsilon_{yy} = v_y + \tfrac12 w_y^2, \qquad
12\gamma_{xy} = u_y + v_x + w_x w_y,
13$$
15the 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}$, where $C = Eh/(1-\nu^2)$ and $G = Eh/(2(1+\nu))$. The equations of motion are
17$$
18\rho h\, u_{tt} = \partial_x N_{xx} + \partial_y N_{xy}, \qquad
19\rho h\, v_{tt} = \partial_x N_{xy} + \partial_y N_{yy},
20$$
22$$
23\rho h\, w_{tt} = -D\nabla^4 w + \partial_x(N_{xx} w_x + N_{xy} w_y) + \partial_y(N_{xy} w_x + N_{yy} w_y) - 2\sigma_0\rho h\, w_t + 2\sigma_1 \rho h\, \nabla^2 w_t ,
24$$
26with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. The quadratic terms in the strains couple bending to stretching. They are what produce the gong-like behaviour: the pitch glide of a strongly struck plate, and the cascade of energy from the low modes that the push excites toward high frequencies. The in-plane inertia terms $\rho h u_{tt}$ and $\rho h v_{tt}$ are usually dropped in favour of an Airy stress function, on the grounds that in-plane waves are much faster than bending waves. We keep them, as requested, at the cost of a much smaller time step (Section 5).
28The loss terms follow the two-parameter model that is standard in plate sound synthesis. A mode of angular frequency $\omega$ decays at the rate $\sigma_0 + \sigma_1 k^2$ with $k^2 = \omega/\kappa$ and $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster. The defaults ($\sigma_0 = 1.15\ \mathrm{s^{-1}}$, $\sigma_1 = 1.08\times10^{-4}\ \mathrm{m^2/s}$) give a 60 dB decay time of about 6 s at low frequencies and about 1 s near 10 kHz. These are not measured values; both parameters are exposed in the interface.
30## 2. Boundary conditions and suspension
32All four edges are 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 the natural boundary conditions of that energy. This avoids the intricate ghost-point formulas that free edges otherwise require, and it guarantees a symmetric stiffness operator.
34The thread through the hole near one corner is not modelled as a constraint. At audio frequencies the thread is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the plate's average position. 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, applied 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}$, which is negligible.
36The same projection defines what "holding" the sheet means during the push. The finger prescribes the displacement of a small patch at the far corner. The reaction to the finger force is distributed over the plate as an inertial load (this is the inertia-relief equilibrium of a free body), so the push produces pure bending of the prescribed amount. This is a modelling 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 about 1.5 mm.
38## 3. Excitation
b9a5a7eSelf-equilibrated release force and loudness normalizationJeremy Magland 40The 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 because 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 and the plate is not accelerated as a whole. The release time can be set to zero to recover the instantaneous case.
42## 4. Spatial discretization
44The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $h_x$ and $h_y$ (about 4 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 and its linear part has a symmetric stiffness matrix.
46**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 they exist (not on the edges perpendicular to the direction of differencing), and $w_{xy}$ by the four-point difference at cell centres, with trapezoidal quadrature weights at edges and corners. 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, which yields the free-edge conditions. The three rigid modes $w = 1, x, y$ lie exactly in the null space.
48**Membrane.** The cells are treated as bilinear (Q4) finite elements with 2 by 2 Gauss integration, using the von Karman strains above. This is the standard element for this purpose; one-point integration would be cheaper but admits hourglass modes, which the nonlinear coupling would excite. The force on a node from each of its four cells is computed once per cell within a GPU workgroup and gathered by the nodes.
50**Mass.** The mass matrix is lumped with trapezoidal weights: a full cell mass in the interior, one half on an edge, one quarter at a corner, so that the total mass is exact.
52**Validation.** The tests in `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, and that the first five frequencies of a completely free square plate agree with the reference values of Leissa and Gorman ($\lambda^2 = 13.47, 19.60, 24.27, 34.80, 34.80$) to within 1 percent at 20 cells per side.
54## 5. Time integration
56The system is advanced with the symplectic Euler (Stormer-Verlet) 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. Its 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 4 mm spacing, $\Delta t \approx 0.65\ \mu$s, roughly six times smaller than bending alone would require, and 4 s of sound takes about 6 million steps. The step is then rounded down so that an integer number of steps (32 at 4 mm and 48 kHz) makes one audio sample.
58Single precision needs one precaution. With $w \approx 10$ mm 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).
60The state is stored twice and the kernel reads one copy and writes the other. Each step is one compute dispatch of 16 by 16 workgroups; a per-step uniform, addressed by a dynamic offset, carries the prescribed corner displacement, the damping, and the audio bookkeeping, so a whole movie frame of steps is encoded in one pass.
62## 6. Radiation to the listener
64The default listener model is the Rayleigh integral for a baffled plate,
66$$
67p(\mathbf r, t) = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \ddot w_k(t - r_k/c),
68$$
70where the sum runs over the nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s. The transverse acceleration 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 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. The "dipole" option is a heuristic for this: each element is weighted by the cosine of its angle to the listener, and the signal is high-pass filtered with a first-order filter at $c/(2\sqrt{L_x L_y})$, about 460 Hz for the default sheet. An accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator.
b9a5a7eSelf-equilibrated release force and loudness normalizationJeremy Magland 72The pressure is reported in pascals, and the peak sound pressure level is displayed. Playback and the downloaded WAV file are normalized either to their peak or, by default, to the loudness of the sustained sound (the RMS after the first 50 ms is scaled to -16 dBFS and a soft clipper bounds the onset), since a signal dominated by the lowest modes has a large peak-to-RMS ratio.
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 73
74## 7. Limitations
76The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; the corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The loss model is empirical and its parameters are not calibrated against a real sheet. The free-edge treatment near the edges is only first-order accurate. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. 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.