# Theory and numerical method

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

## 1. Plate model

The 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

$$
\varepsilon_{xx} = u_x + \tfrac12 w_x^2, \qquad
\varepsilon_{yy} = v_y + \tfrac12 w_y^2, \qquad
\gamma_{xy} = u_y + v_x + w_x w_y,
$$

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}$, where $C = Eh/(1-\nu^2)$ and $G = Eh/(2(1+\nu))$. The equations of motion are

$$
\rho h\, u_{tt} = \partial_x N_{xx} + \partial_y N_{xy}, \qquad
\rho h\, v_{tt} = \partial_x N_{xy} + \partial_y N_{yy},
$$

$$
\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 ,
$$

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

The loss terms follow the two-parameter model that is standard in plate sound synthesis. A bending 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 = 5\times10^{-4}\ \mathrm{m^2/s}$) give a 60 dB decay time of about 5 s at 100 Hz, 1.8 s at 1 kHz, and 0.5 s at 5 kHz. The value of $\sigma_1$ is chosen to match an estimate of the radiation damping of a thin steel sheet in air, which for this thickness corresponds to a loss factor of roughly $10^{-3}$ at 1 kHz; radiation, not the metal itself, is the main loss of a thin sheet at mid and high frequencies. The in-plane motion carries a separate frequency-independent loss $\sigma_{\mathrm{in}}$ (default $20\ \mathrm{s^{-1}}$). It matters because the stretching modes of the sheet lie at several kilohertz, couple linearly to transverse motion once the sheet is curved (Section 1, initial shape), 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 are measured values, and all three are exposed in the interface.

### Initial shape (imperfection)

Real sheet metal is never flat; rolling and handling leave a bow of a few millimetres. A bowed sheet is a shallow shell, and its behaviour 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 the wobble board, and it appears to be what produces the slow wobble of a real sheet. We therefore allow a stress-free initial shape $w_0(x,y)$ in the Marguerre form of the von Karman equations: the membrane strains use the total slope minus the stress-free slope,

$$
\varepsilon_{xx} = u_x + \tfrac12\left[(w_0 + w)_x^2 - w_{0,x}^2\right], \qquad
\varepsilon_{yy} = v_y + \tfrac12\left[(w_0 + w)_y^2 - w_{0,y}^2\right], \qquad
\gamma_{xy} = u_y + v_x + (w_0 + w)_x (w_0 + w)_y - w_{0,x} w_{0,y},
$$

the membrane forces act through the total slope $(w_0 + w)$, and bending acts on the deviation $w$ alone. The shape $w_0 = 0$, $w = 0$ carries no stress. Three shapes are offered, all parabolic with a prescribed sag of the centre relative to the edges: curved across the width, curved across the height, and a dish curved both ways. The default is a 3 mm bow across the width. A positive push moves the corner in the direction of the bow and a negative push against it; the two cases sound different, as they do on a real sheet.

## 2. Boundary conditions and suspension

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

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

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

## 3. Excitation

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

## 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$ (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.

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

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

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

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

## 5. Time integration

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

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

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

## 6. Radiation to the listener

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

$$
p(\mathbf r, t) = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \ddot w_k(t - r_k/c),
$$

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

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

## 7. Limitations

The 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 Marguerre treatment of the initial shape is likewise a shallow-shell approximation, adequate for a bow of a few millimetres, 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. 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.