71e94accymsim: nonlinear FEM cymbal simulation packageJeremy Magland 1# The cymsim model: equations, discretization and method
e30c940save .cymsim fileJeremy Magland 3This document states what cymsim solves and how. Section 2 gives the continuous model, section 3 the geometry, section 4 the spatial discretization, section 5 the time integration, section 6 the excitation, section 7 the radiation model, section 8 the diagnostics and section 9 the records kept for other programs. Section 10 lists the limitations we are aware of. Each section names the module that implements it. Units are SI throughout.
71e94accymsim: nonlinear FEM cymbal simulation packageJeremy Magland 4
5## 1. Overview
7A crash cymbal is a thin, shallow, curved bronze shell. Struck hard, it does not simply ring at its linear modes: the transverse motion is many times the thickness, membrane stretching couples the modes, and energy injected at low frequency by the stroke climbs to high frequency over tens of milliseconds. This upward cascade is what gives a crash its build-up and its noisy, bright sustain, and it is entirely a nonlinear effect. A model that is meant to reproduce it must therefore keep the geometric nonlinearity of the shell, resolve frequencies up to the audio band, and run long enough for the cascade to develop.
9cymsim does this with the smallest set of ingredients we found sufficient: the Marguerre-von Karman shallow-shell equations (the von Karman plate with an initial curvature), a flat-facet triangular finite-element discretization, an explicit leapfrog integrator with a stability-based sub-step rate, a rigid exciter coupled through a power-law contact spring, and a Rayleigh-integral listener. Everything from the element energies onward is written with JAX so that the internal force is the exact gradient of the discrete strain energy and the whole time step compiles to one kernel.
11## 2. Governing equations
13Implemented in [elements.py](../src/cymsim/elements.py) (strains, energies) and [shell.py](../src/cymsim/shell.py) (damping, mass).
15### 2.1 Kinematics
17Let $\Omega$ be the disk of radius $R$ in the $(x, y)$ plane. The undeformed mid-surface is the graph $z = z_0(x, y)$ of a *shallow* surface, meaning $|\nabla z_0| \ll 1$; we warn when $\max |\nabla z_0|$ exceeds 0.4. The displacement of a mid-surface point has in-plane components $u(x, y, t)$, $v(x, y, t)$ along $x$ and $y$ and a transverse component $w(x, y, t)$ along $z$. The Marguerre-von Karman membrane strains are
19$$
20\varepsilon_{xx} = u_{,x} + \tfrac12 w_{,x}^2 + z_{0,x}\, w_{,x}, \qquad
21\varepsilon_{yy} = v_{,y} + \tfrac12 w_{,y}^2 + z_{0,y}\, w_{,y}, \qquad
22\gamma_{xy} = u_{,y} + v_{,x} + w_{,x} w_{,y} + z_{0,x}\, w_{,y} + z_{0,y}\, w_{,x},
23$$
25and the curvatures are those of the plate, $\kappa = (w_{,xx},\, w_{,yy},\, 2 w_{,xy})$. Setting $z_0 = 0$ recovers the von Karman plate. Dropping the quadratic terms $\tfrac12 w_{,x}^2$, $\tfrac12 w_{,y}^2$, $w_{,x} w_{,y}$ gives the *linear shallow shell* with the same geometry; the configuration flag `numerics.nonlinear = false` does exactly this and is the reference against which the cascade is judged. The terms in $z_0$ are the linear coupling between transverse motion and membrane stretch that a curved shell has and a flat plate lacks; they are what make the bell and the bow stiffen the low modes.
27### 2.2 Constitutive relations and strain energy
29The material is isotropic with Young's modulus $E$, Poisson ratio $\nu$ and density $\rho$; the thickness $h(x, y)$ varies. Membrane forces and bending moments are $N = C \varepsilon$ and $M = D \kappa$ with
31$$
32C = \frac{E h}{1 - \nu^2} \begin{pmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \tfrac{1 - \nu}{2} \end{pmatrix}, \qquad
33D = \frac{E h^3}{12 (1 - \nu^2)} \begin{pmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \tfrac{1 - \nu}{2} \end{pmatrix}.
34$$
36The strain energy is
38$$
39U[u, v, w] = \frac12 \int_\Omega \left( \varepsilon^{\mathsf T} C\, \varepsilon + \kappa^{\mathsf T} D\, \kappa \right) dA ,
40$$
42the first term being the membrane energy and the second the bending energy. Their ratio during a run is one of the diagnostics of how nonlinear the motion is.
44### 2.3 Equations of motion and damping
46With $\rho h$ the mass per unit area and $\delta U / \delta(\cdot)$ the variational derivative,
48$$
49\rho h\, \ddot u = -\frac{\delta U}{\delta u} - 2 \sigma_0 \rho h\, \dot u, \qquad
50\rho h\, \ddot v = -\frac{\delta U}{\delta v} - 2 \sigma_0 \rho h\, \dot v, \qquad
51\rho h\, \ddot w = -\frac{\delta U}{\delta w} - 2 \sigma_0 \rho h\, \dot w + 2 \sigma_1 \rho h\, \nabla^2 \dot w + f_c ,
52$$
54where $f_c$ is the contact pressure of the exciter (section 6). Rotary inertia is neglected in the continuum model (it is reintroduced in modified form by the discretization, section 4.6). The losses follow the two-parameter form used throughout the sound-synthesis literature (Bilbao 2009): $\sigma_0$ gives a frequency-independent decay rate, so that the 60 dB decay time of the low modes is $T_{60} \approx 6.91 / \sigma_0$, and the $\sigma_1$ term gives a wave of wavenumber $k$ the additional decay rate $2 \sigma_1 k^2$. For bending waves $k^2 = \omega \sqrt{\rho h / D}$, so this rate grows linearly with frequency and sets how far up the cascade can reach before losses win. We make no attempt to model the physical loss mechanisms (thermoelastic, viscoelastic, radiation) separately; $\sigma_0$ and $\sigma_1$ are tuning parameters that lump them.
56### 2.4 Support
58The shell is free at the rim. A real cymbal rests on a felt washer at the bell, which we do not model; instead the six rigid-body motions of the shell are removed by projection at every step (section 5.3). This is the simplest device that keeps the cymbal from drifting away under the impulse of the stroke.
60## 3. Geometry
62Implemented in [geometry.py](../src/cymsim/geometry.py); parameters in `Profile` and `Imperfections`.
64The profile is written in polar coordinates $(r, \phi)$ on the projected disk. The *shoulder* is the junction between the bell and the bow; its radius $r_{sh}(\phi) = R_{sh} (1 + \epsilon_{oval} \cos(2\phi + \phi_{oval}))$ carries an $m = 2$ ovality. The bow coordinate $s = \mathrm{clip}\big((R - r) / (R - r_{sh}(\phi)),\, 0,\, 1\big)$ runs from 0 at the rim to 1 at the shoulder, and $u = 1 - s$. With $r_b$ the distance from the (offset) bell centre and $b(r_b) = \tfrac12 (1 + \cos \pi t)$, $t = \mathrm{clip}\big((r_b - R_{top}) / (r_{sh} - R_{top}),\, 0,\, 1\big)$, a $C^1$ ramp that is 1 on the bell plateau and 0 beyond the shoulder, the mid-surface height is
66$$
67z_0 = H_{bow} \left( 1 - u^{p} \right) + H_{bell}\, b(r_b) + u \sum_m a_m \cos(m \phi + \phi_m) + \sum_i A_i \exp\!\left( -\frac{|x - x_i|^2}{2 s_i^2} \right),
68$$
70bow, bell, rim warp and hammer dents in that order. The bow exponent $p = 1$ gives a cone and $p > 1$ a bow that is curved near the rim and flat at the shoulder. The warp amplitudes are $a_m = 2 a_{warp} / m$ for the listed modes with random phases; it vanishes at the shoulder and grows linearly to the rim, so the mean rim height stays $z_0 = 0$. The dents are Gaussian bumps of random sign (a cymbal is hammered from both sides), with centres uniform over the annulus between the shoulder and the rim and with jittered depth and width. The thickness is
72$$
73h = h_{rim} + (h_{sh} - h_{rim})\, s + (h_{bell} - h_{sh})\, b(r_b),
74$$
76a linear taper along the bow and a smooth step into the bell. All random draws come from `numpy.random.default_rng(seed)` in a fixed order, so a run is reproducible from its configuration. The imperfections matter: a perfectly axisymmetric shell has degenerate mode pairs and symmetric coupling rules, and its cascade is markedly less rich than that of a real cymbal.
78## 4. Spatial discretization
80### 4.1 Mesh
82Implemented in [mesh.py](../src/cymsim/mesh.py). The finite elements live on the flat projected disk; the curvature enters only through $z_0$ in the strains. Nodes are placed on concentric rings of spacing $\Delta r = R / K$ with $K = \mathrm{round}(R / \ell)$ for a target element size $\ell$; ring $k$ carries $n_k = \mathrm{round}(2\pi r_k / \Delta r) \approx 6k$ nodes and alternate rings are rotated by half a step. Consecutive rings are stitched by an angular merge: walking both rings in increasing angle, the next triangle uses whichever ring's next node comes first, which yields $n_{k-1} + n_k$ triangles per annulus and a valid triangulation. A Laplacian smoothing pass (`mesh.smoothing`, default 1) then relaxes the interior nodes with the rim nodes held on the circle. The mean element quality (inradius to circumradius ratio, 1 for an equilateral triangle) is about 0.96 and the minimum about 0.7. More passes improve the mean shape but shrink the smallest elements and thereby the stable time step, by about 10 percent for ten passes; the explicit step is set by the smallest elements in the thick bell (section 5.2), so a mesh graded coarser in the bell would be a worthwhile refinement. The rim is the polygon through the outer ring, so the mesh area is slightly below $\pi R^2$.
84### 4.2 Degrees of freedom
86Each node $i$ carries five unknowns $q_i = (u_i, v_i, w_i, \theta_{x,i}, \theta_{y,i})$ with the DKT convention $\theta_x = w_{,y}$ and $\theta_y = -w_{,x}$. The state of the shell is the array $q$ of shape $(N, 5)$.
88### 4.3 Membrane part: constant-strain triangle with the discrete initial surface
90On each triangle $u$, $v$ and $w$ are interpolated linearly from their nodal values, so their gradients are constant: $\partial_x \phi = \sum_i G_{x,i} \phi_i$ with $G_x$, $G_y$ the usual CST operators. The initial surface is represented in the same space, that is, $z_{0,x}$ and $z_{0,y}$ on an element are the CST gradients of the nodal heights $z_0(x_i, y_i)$, not the analytic slope at the centroid. This choice is deliberate. A rigid tilt of the shell about the $x$ axis by an angle $a$ displaces a surface point by $v = -a z_0$ and $w = a y$ (section 4.7); its linearized strain is $\varepsilon_{yy} = v_{,y} + z_{0,y} w_{,y} = -a\, (G z_0) + a\, z_{0,y}$, which vanishes exactly only if $z_{0,y}$ is the CST gradient of the nodal heights. With the analytic slope, the mismatch on a dented, coarse element left a rigid tilt with a spurious restoring force of order $10^{-4}$ of a typical stiffness; with the discrete slope the residual is at round-off (about $10^{-16}$ in float64). The membrane energy of element $e$ is $U_{m,e} = \tfrac12 A_e\, \varepsilon^{\mathsf T} C_e\, \varepsilon$ with $C_e$ evaluated at the centroid thickness.
92### 4.4 Bending part: the Discrete Kirchhoff Triangle
94Bending uses the DKT element of Batoz, Bathe and Ho (1980). Rotations $\beta_x, \beta_y$ are interpolated quadratically over the triangle from the three nodal $(w, \theta_x, \theta_y)$ triples, with the mid-side values fixed by imposing the Kirchhoff condition (rotation equals slope) discretely along each edge. The resulting curvature operator $B(\xi, \eta)$, with $\kappa = B\, U_e$ for the nine bending DOFs $U_e$ of the element, is linear in the area coordinates, so the element stiffness
96$$
97K_e = \int_{A_e} B^{\mathsf T} D_e\, B\; dA
98$$
100is integrated exactly by the three-point rule. The bending energy is $U_{b,e} = \tfrac12 U_e^{\mathsf T} K_e U_e$. DKT does not define a transverse displacement field inside the element, only the nodal values, which is why the membrane strains use the CST interpolation of $w$; the pairing of a DKT bending triangle with a CST membrane triangle is the classical flat shell triangle, and the Marguerre terms supply the curvature coupling between them. The element stiffness is computed once in float64 by automatic differentiation of the shape functions and checked at build time: a rigid tilt gives zero force (residual about $10^{-17}$), the state $w = x^2/2$ has energy $\tfrac12 D A_e$ to $10^{-10}$, and the $15 \times 15$ tangent of the full element energy has exactly six zero eigenvalues (three translations and three rotations). `cymsim check` prints these.
102### 4.5 Internal force by automatic differentiation
104The discrete strain energy is $U_h(q) = \sum_e \big( U_{m,e}(q_e) + U_{b,e}(q_e) \big)$ with $q_e$ the $(3, 5)$ gather of the element's nodal values. The internal nodal force is
106$$
107f_{int}(q) = -\frac{\partial U_h}{\partial q},
108$$
110obtained by reverse-mode differentiation (`jax.grad`) of the summed energies; the gather, the per-element evaluation and the scatter are all inside the differentiated function. The force is therefore the exact gradient of the discrete energy, so the semi-discrete system is conservative up to the damping terms, and the energy diagnostic is meaningful. The tangent stiffness needed for the stability bound is the Hessian of the same function at $q = 0$.
112### 4.6 Lumped mass and rotational inertia
114The mass is lumped: element $e$ contributes $m_e = \rho h_e A_e / 3$ to the $u$, $v$ and $w$ entries of each of its nodes. The rotational entries receive
116$$
117I_e = m_e \left( \frac{h_e^2}{12} + f A_e \right),
118$$
120where $h_e^2 / 12$ is the physical rotary inertia per unit mass and $f A_e$ (default $f = 1/64$) is an added inertia. The addition is a numerical device: with the physical value alone, the rotational DOFs of DKT have eigenfrequencies far above any bending wave the mesh can represent and would dictate an unusably small time step. The price is that bending waves whose wavelength is comparable to the element size are slowed slightly, by a relative amount of order $(k h_{eff})^2 / 24$ with $h_{eff}^2 = h_e^2 + 12 f A_e$, which for 6 mm elements is a few percent at 20 kHz and well below one percent at 5 kHz. These are also the waves the mesh resolves least well, so we consider the trade acceptable; $f$ can be reduced at the cost of a smaller time step. A second, optional scaling `inplane_mass_scale` multiplies the in-plane mass; values above 1 slow the fast in-plane waves ($c_L \approx 3.8$ km/s in bronze) and relax the time step further, but they also alter the in-plane dynamics that mediate the nonlinear coupling, so the default is 1.
122### 4.7 Rigid-body modes of the curved shell
124The free shell has six rigid-body modes. Three are translations. The in-plane rotation is $(u, v) = (-y, x)$. A tilt about the $x$ axis by an angle $a$ carries the surface point $(x, y, z_0)$ to first order to $(x,\, y - a z_0,\, z_0 + a y)$, so the mode is $v = -z_0$, $w = y$, $\theta_x = 1$; likewise the tilt about $y$ is $u = z_0$, $w = -x$, $\theta_y = 1$. The in-plane components proportional to $z_0$ are what distinguish the curved shell from the plate and are essential: without them a tilt stretches the shell. The six vectors are made $M$-orthonormal by Gram-Schmidt. At build time we verify that the tangent stiffness annihilates all six to round-off. Note that this holds for the linearized operator only: the Marguerre-von Karman strains are invariant under a finite rotation only to first order, so a finite tilt of amplitude $a$ leaves a residual force of order $a^2$. This is a known property of the shallow-shell approximation, not of the discretization, and is harmless here because the projection of section 5.3 keeps the rigid amplitudes at zero.
126### 4.8 The damping operator
128The $\sigma_1$ term is discretized in weak form with the CST gradients. The nodal force on $w$ is
130$$
131f_{\sigma_1, i} = 2 \sigma_1 \sum_{e \ni i} (\rho h A)_e\, G_{e,i}^{\mathsf T} \big( G_e\, \dot w_e \big),
132$$
134which is the assembly of $2 \sigma_1 \rho h \int \nabla \dot w \cdot \nabla \phi_i\, dA$. Its power $\dot w^{\mathsf T} f_{\sigma_1} = 2 \sigma_1 \sum_e (\rho h A)_e |\nabla \dot w|^2 \geq 0$ is non-negative, so the term only ever removes energy, and a uniform velocity is undamped by it. The $\sigma_0$ term is diagonal and acts on all five DOFs.
136## 5. Time integration
138Implemented in [simulate.py](../src/cymsim/simulate.py).
140### 5.1 Leapfrog with implicit $\sigma_0$
142With $M$ the diagonal lumped mass, one sub-step of size $\Delta t$ is
144$$
145v^{n+1} = P \left[ \frac{v^n + \Delta t\, M^{-1} \big( f_{int}(q^n) + f_c(q^n, z_s^n) - f_{\sigma_1}(v^n) \big)}{1 + 2 \sigma_0 \Delta t} \right], \qquad
146q^{n+1} = q^n + \Delta t\, v^{n+1},
147$$
149where $P$ is the rigid-mode projection of section 5.3. This is the symplectic Euler scheme in $(q, v)$, equivalent to the central-difference (Stormer-Verlet) scheme in $q$, second-order accurate and non-dissipative for the conservative part. The $\sigma_0$ term is treated implicitly, which costs nothing because it is diagonal; the $\sigma_1$ term is evaluated explicitly at $v^n$, which is accurate enough because $2 \sigma_1 k^2 \Delta t \ll 1$ for every wavenumber the mesh can carry at the default values. The exciter is advanced by the same symplectic Euler step (section 6.3).
151### 5.2 Stability and the choice of the time step
e30c940save .cymsim fileJeremy Magland 153For the linear undamped system the scheme is stable when $\Delta t \leq 2 / \omega_{max}$, with $\omega_{max}^2$ the largest eigenvalue of $M^{-1} K$. Instead of solving that global eigenproblem we use the element bound (Irons; see Hughes 2000): $\lambda_{max}(M^{-1} K) \leq \max_e \lambda_{max}(M_e^{-1} K_e)$, where $K_e$ is the $15 \times 15$ tangent stiffness of the element at rest and $M_e$ its share of the lumped mass. The bound is computed in float64 at build time and the limiting element (the thinnest, near the rim) is reported. The nonlinear membrane stiffening raises the instantaneous frequencies during large motion, so we take only a fraction `cfl_safety` (default 0.6) of the bound. The sub-step count per audio sample is $m = \lceil 1 / (f_s\, c_{safety}\, \Delta t_{bound}) \rceil$ and $\Delta t = 1 / (m f_s)$, so that the radiation buffers decimate exactly to the audio rate. The run is extended by $W/2 = 16$ audio samples beyond the requested duration so that the last samples of the decimated records (section 9) are final. If the energy diagnostic grows during a run the step is too large; the remedy is a smaller `cfl_safety`. The run aborts if the energy becomes non-finite.
71e94accymsim: nonlinear FEM cymbal simulation packageJeremy Magland 154
155### 5.3 Rigid-mode projection
157After each velocity update, $v \leftarrow v - \sum_{k=1}^{6} r_k \big( r_k^{\mathsf T} M v \big)$ with the $M$-orthonormal rigid modes $r_k$ of section 4.7. This removes the linear and angular momentum that the stroke imparts, which is what a stand would do, and it also removes the slow rigid drift that round-off would otherwise accumulate. The projection changes only the rigid part of the velocity and hence removes only rigid kinetic energy; the vibrational energy is untouched. We emphasize that this is a proxy for the support, not a model of it: a felt washer at the bell would also add a local stiffness and some damping.
159### 5.4 Precision and compilation
161Set-up (element stiffness, eigenvalue bound, rigid modes) is done in float64; the time loop runs in float32 by default (`numerics.dtype`), which is what a GPU does quickly. At the default parameters the displacement increment $\Delta t\, v$ is about $10^{-3}$ of the displacement itself, so float32 loses a few of its 24 bits per step in $q$ without visible effect on the energy balance; float64 is available for verification. The step is compiled once with `jax.jit` as a `fori_loop` over `steps_per_block` sub-steps, and the Python loop runs over blocks, reporting energies between them.
163## 6. Excitation
e19a540Add a prescribed-force exciter and point probes of the normal acceleration; accept Python 3.9 in the YAML loaderJeremy Magland 165Implemented in [exciter.py](../src/cymsim/exciter.py); parameters in `Hit`, `StickExciter`, `MalletExciter`, `ForceExciter`.
71e94accymsim: nonlinear FEM cymbal simulation packageJeremy Magland 166
167### 6.1 Contact law
169Each exciter is a rigid body of mass $m$ that moves only along $z$, with face height $z_s(t)$, coupled to the shell through $P$ springs at patch points $p$. The surface height under point $p$ is $w_p + o_p$, where $w_p = \sum_i \lambda_{p,i} w_i$ is the barycentric interpolation of the nodal transverse displacements in the containing element and $o_p$ is the height of the initial surface above the exciter's face plane, the best-fit plane through the patch points' $z_0$ values. The penetration and the force at point $p$ are
171$$
172\delta_p = \max(w_p + o_p - z_s,\, 0), \qquad F_p = \frac{k}{P}\, \delta_p^{\alpha}, \qquad F = \sum_p F_p .
173$$
175The plate receives $-F_p$ on the $w$ DOFs of the containing element, distributed with the weights $\lambda_{p,i}$, and the exciter obeys $m \ddot z_s = F$. The exponent $\alpha = 1.5$ is the Hertz law for a hard sphere on a hard surface; $\alpha \approx 2$ to $3$ fits felt. The contact spring is conservative (potential $k \delta^{\alpha+1} / (P(\alpha + 1))$), so the exciter rebounds with whatever energy the shell has not absorbed; there is no contact hysteresis, which slightly overestimates the rebound of a felt mallet. The residual offsets $o_p$ mean that the face touches the curved surface at a few points first, which is what a rigid face on a curved shell does.
177### 6.2 Contact patches
179A `StickExciter` contacts along a line of length $\ell$ through the strike point, at an angle `patch_angle` from the radial direction, with $P$ equally spaced springs and a given total stiffness $k$. A `MalletExciter` contacts on a disk: one spring at the centre and six on a hexagon of the head radius, and its stiffness is set from a target compression $\delta_c$ at full stop on a rigid wall, $\tfrac12 m v^2 = k \delta_c^{\alpha + 1} / (\alpha + 1)$. The set-up log reports, for each hit, the contact duration the exciter would have against a rigid wall,
181$$
182\tau = \frac{2 c(\alpha)\, \delta_{max}}{v}, \qquad c(\alpha) = \int_0^1 \frac{dx}{\sqrt{1 - x^{\alpha + 1}}},
183$$
185which for $\alpha = 1.5$ is Hertz's $2.94\, \delta_{max} / v$, together with the frequency $1 / (2\tau)$ below which the stroke itself injects energy. For the default mallet this is about 120 Hz, so everything the listener hears above a few hundred hertz has come up the cascade. On the shell the actual contact is longer and softer than on a rigid wall because the surface gives way.
187### 6.3 Schedule of hits
189A run has a list of hits, each with its own time, position, velocity and exciter. Every exciter exists for the whole run as a row of fixed-shape arrays (hits $\times$ largest patch, padded with zero weights and a very negative offset so that padded points never touch), which lets one compiled step serve any number of hits. Before its release step $n_k = \mathrm{round}(t_k / \Delta t)$ an exciter is inactive: it produces no force and is held far above the shell. At step $n_k$ it is placed at $z_s = \max_p (w_p + o_p) + \text{gap}$, just above the *current* surface, with velocity $-v_k$. This matters for a train of hits, where the surface under a later hit is already moving by several thickness units. After the rebound the exciter simply moves away at constant speed; there is no gravity and no hand. Two consequences follow: an exciter never strikes twice on its own, and if the shell were to overtake a slowly rebounding exciter a second contact would occur (a real player pulls the stick back). The summary of each hit reports the peak force, the total time with $F > 0$ and the rebound velocity.
e19a540Add a prescribed-force exciter and point probes of the normal acceleration; accept Python 3.9 in the YAML loaderJeremy Magland 191### 6.4 Prescribed force
193A `ForceExciter` replaces the contact by a force given in advance, a raised cosine of peak $F_0$ and length $T$ applied at one point from the hit time $t_k$:
195$$
196F(t) = \tfrac12 F_0 \left(1 - \cos \frac{2\pi (t - t_k)}{T}\right), \qquad t_k \le t \le t_k + T,
197$$
199distributed to the containing element's $w$ DOFs with the barycentric weights and pushing the surface down. The force does not depend on the motion, so there is no rebound, and the hit's velocity is ignored. This is not a model of a stick or a mallet. It is the excitation used by modal plate models (Ducceschi and Touzé 2015, among others), and it exists so that cymsim can be run on exactly their problems.
e30c940save .cymsim fileJeremy Magland 203Implemented in [radiation.py](../src/cymsim/radiation.py) and [decimation.py](../src/cymsim/decimation.py); parameters in `Listener`.
205The pressure at the listener point $x_L$ from a baffled planar radiator is given by the Rayleigh integral, which we discretize with one point source per element:
207$$
e30c940save .cymsim fileJeremy Magland 208p(x_L, t) = \frac{\rho_{air}}{2\pi} \sum_e \frac{A_e\, \bar a_e\!\left( t - r_e / c \right)}{r_e},
e30c940save .cymsim fileJeremy Magland 211where $\bar a_e$ is the mean over the element's three nodes of the *normal acceleration* $a_i = n_i \cdot (\ddot u_i, \ddot v_i, \ddot w_i)$, $n_i$ being the unit normal of the rest surface at node $i$ (the area-weighted mean of the adjacent elements' constant slopes, pointing to the $+z$ side), and $r_e$ the distance from the element centroid, placed on the curved surface at height $z_0$, to the listener. The acceleration is taken from the integrator directly: for the leapfrog scheme $(v^{n+1} - v^n) / \Delta t$ is the exact central second difference of $q$ at step $n$, so no differentiation of stored data is involved. This is exact for a rigid piston in an infinite baffle and an approximation here: the surface is shallow rather than flat, and the air does not load the shell (radiation damping is part of what $\sigma_0$ and $\sigma_1$ lump). At the bell, where the slope reaches 0.39, the normal component differs from $\ddot w$ by several percent, which is why we use it rather than $\ddot w$.
71e94accymsim: nonlinear FEM cymbal simulation packageJeremy Magland 212
213A cymbal has no baffle. Its back face radiates with the opposite sign, and that sound reaches the listener by diffracting around the rim; at low frequency the two contributions nearly cancel, which is why an unbaffled plate is a weak low-frequency source. We estimate the back contribution with the same sum but with $r_e$ replaced by the shortest path over the rim, $r_e^{back} = \min_j \big( |x_e - x_j^{rim}| + |x_j^{rim} - x_L| \big)$ over 720 rim points, weighted by `back_diffraction_gain`, and subtract it: $p_{unbaffled} = p_{front} - g\, p_{back}$. This is a Kirchhoff-style edge-diffraction estimate, not a boundary-element solution, and it ignores shadowing of the front by the shell itself. Both signals are returned; the baffled one is a cleaner measure of the shell motion, the unbaffled one sounds more like a cymbal in a room.
e30c940save .cymsim fileJeremy Magland 215Retarded times are handled at sub-step resolution: the delay of element $e$ in sub-steps, $d_e = r_e / (c \Delta t)$, is split into an integer part and a fraction, and each element's contribution is added to two adjacent slots of a ring buffer with linear interpolation. After each block the host reads the slots of the finished sub-steps, clears them and decimates them to the audio rate with the kernel of section 9, the same kernel that filters the acceleration record. Since delay and filtering are both linear and time-invariant, the stored listener signal follows from the stored accelerations up to the linear-interpolation error at the sub-step scale, which is negligible; the test suite rebuilds the baffled signal from the record to $1.5 \times 10^{-4}$ relative rms.
71e94accymsim: nonlinear FEM cymbal simulation packageJeremy Magland 216
217## 8. Diagnostics
219Implemented in [simulate.py](../src/cymsim/simulate.py) (energies) and [diagnostics.py](../src/cymsim/diagnostics.py) (spectra). During the run the kinetic, membrane, bending and exciter energies are logged at `n_progress` points together with $\max_i |w_i| / h_i$, the *nonlinearity meter*: below about 1 the run is effectively linear, at 3 to 5 the cascade is clearly developed. The sum of the four energies can only decrease between the release of a hit and the next one, which is a direct check of the integrator. After the run, an STFT of the unbaffled pressure (256-point windows, 32-sample hop, about 0.7 ms) gives per-octave band energy envelopes from 60 Hz to 16 kHz; the *onset* of a band is the first time it reaches 20 dB below its own peak, and the spectral centroid is tracked against time. In a cascade the onsets are staggered with the high bands last and the centroid climbs over tens of milliseconds; in the linear reference run the bands rise together and the centroid is flat.
223Implemented in [decimation.py](../src/cymsim/decimation.py), [simulate.py](../src/cymsim/simulate.py) and [cymsim_file.py](../src/cymsim/cymsim_file.py); parameters in `Output`; file layout in [format.md](format.md).
225The listener model of section 7 is deliberately simple, and a natural use of the simulation is to feed the shell's motion to a different acoustic model: another listener position, a room, or a solver of the acoustic wave equation with the shell as a vibrating boundary. What such a model needs is the normal acceleration of the surface at the audio rate, and the question is how to obtain it without noise. Differentiating stored displacements twice is the wrong route: a second difference multiplies noise by $\omega^2$, so the top of the band suffers most, it folds whatever the shell does above the audio Nyquist frequency back into the band, and in float32 the three-term difference of nearly equal numbers loses most of its digits. We therefore take the acceleration from the integrator, where it is exact for the scheme (section 7), and the only remaining step is the reduction from the sub-step rate, typically 2.4 MHz, to the audio rate.
227That reduction has to be a proper low-pass filter, because the cascade puts real energy above 20 kHz and acceleration weights it by $\omega^2$. We use a Kaiser-windowed sinc $h(d)$ spanning $W = 32$ audio samples ($Wm$ sub-steps) with $\beta = 8$ and its $-6$ dB point at $0.46 f_s$: flat to about $0.38 f_s$, $-23$ dB at $f_s / 2$, below $-80$ dB from $0.54 f_s$. It is applied inside the compiled time loop through $W$ accumulator slots per node: at sub-step $j_0 m + kk$ the node adds $a_i\, h(kk + (W/2 - 1 - w)\, m)$ to slot $w$, which belongs to audio sample $j_0 + w - W/2 + 1$; when the sub-step count crosses a sample boundary the oldest slot is complete, is written to the block's output and is reused for the newest sample. The cost is $W$ multiply-adds per node per sub-step, small next to the element force evaluation. The kernel is linear-phase and the read-out removes its delay, so sample $j$ is centred on $t = j / f_s$; the price is that the newest $W/2$ samples are not final, which is why the run is extended by $W/2$ samples (section 5.2) and why a partial save stops $W/2$ samples short of the current time.
229The record is kept only on request (`output.record_acceleration`, default off), because it is large: one value per node per audio sample, stored as `int16` against one `float32` scale per sample, is $2 N f_s$ bytes per simulated second, about 0.4 GB per second at the default mesh. The per-sample scale keeps the quantisation noise at $1.5 \times 10^{-5}$ of each sample's own peak through the 90 dB or more of a stroke's decay; the noise is incoherent across nodes and changes the built-in Rayleigh sum by less than $10^{-4}$ relative rms. Transverse displacements, in contrast, are kept only at a low frame rate (`output.displacement_rate`, default 60 Hz), in float32 and exactly at the frame times. At 60 Hz they are a strobe of the motion, useful for display and for the slow quasi-static deflection, but not a signal from which anything should be differentiated.
e19a540Add a prescribed-force exciter and point probes of the normal acceleration; accept Python 3.9 in the YAML loaderJeremy Magland 231For comparisons at a single point, `output.probes` lists points $(x, y)$ on the projected domain where the normal acceleration is kept at the audio rate in float64: the nodal values are interpolated with the barycentric weights of the containing element, held for one block at the sub-step rate, and decimated on the host with the same kernel. A probe placed on a node reproduces that node's column of the acceleration record. Probe signals are returned in the result and its `.npz` file but are not yet written to the `.cymsim` file.
e30c940save .cymsim fileJeremy Magland 233Both records, the two listener signals, the configuration, the mesh with its nodal normals, the energy log, the hit summaries and the integrator state are written to an HDF5 file with the extension `.cymsim`, described in [format.md](format.md). The time series are resizable, so `cymsim run -o run.cymsim --save-every 60` writes a consistent partial file every minute and extends it, and a run that is interrupted leaves a loadable file. The web companion writes and reads the same format.
235## 10. Limitations and future work
237The model is a shallow shell, so it is only as good as the assumption $|\nabla z_0| \ll 1$; the default profile reaches a slope of 0.39 at the bell, which is at the edge of what shallow-shell theory is meant for, and a steep bell would call for a proper curved-shell formulation. The elements are flat facets with the curvature entering through $z_0$ alone, and the membrane strain uses the linear interpolation of $w$ while bending uses DKT rotations; this is a standard pairing but not a fully consistent one, and its convergence with element size has not been studied here. The lumped rotational inertia is inflated for the sake of the time step (section 4.6). The support is a rigid-mode projection rather than a model of the stand. The damping has two parameters and no physics behind them. The contact is one-dimensional along $z$, with no friction, no tangential force, no hysteresis and no hand. The radiation is a Rayleigh sum with a heuristic back-side correction.
239The natural next steps are, roughly in order of expected benefit: a convergence study of the cascade statistics with element size and precision; a bell support with local stiffness and damping in place of the projection; a contact model with hysteretic losses for felt; and a proper unbaffled radiation model. A train of hits already works structurally, but the behaviour of the contact model against a strongly moving surface deserves its own validation.
243- Batoz, J.-L., Bathe, K.-J., Ho, L.-W. (1980). A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15, 1771-1812.
244- Marguerre, K. (1938). Zur Theorie der gekrummten Platte grosser Formanderung. Proceedings of the 5th International Congress for Applied Mechanics, 93-101.
245- Bilbao, S. (2009). Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics. Wiley. (Plate and shell schemes; the $\sigma_0$, $\sigma_1$ loss model.)
246- Chaigne, A., Touze, C., Thomas, O. (2005). Nonlinear vibrations and chaos in gongs and cymbals. Acoustical Science and Technology 26, 403-409.
247- Touze, C., Bilbao, S., Cadot, O. (2012). Transition scenario to turbulence in thin vibrating plates. Journal of Sound and Vibration 331, 412-433.
248- Ducceschi, M., Touze, C. (2015). Modal approach for nonlinear vibrations of damped impacted plates: application to sound synthesis of gongs and cymbals. Journal of Sound and Vibration 344, 313-331.
249- Hughes, T. J. R. (2000). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover. (Element eigenvalue bound for explicit time stepping.)
250- Hairer, E., Lubich, C., Wanner, G. (2006). Geometric Numerical Integration. Springer. (Symplectic Euler and Stormer-Verlet.)
251- Pierce, A. D. (1989). Acoustics: An Introduction to Its Physical Principles and Applications. Acoustical Society of America. (Rayleigh integral.)
255| symbol | meaning | configuration / code |
256|---|---|---|
257| $E, \nu, \rho$ | Young's modulus, Poisson ratio, density | `Material.E`, `.nu`, `.rho` |
258| $\sigma_0, \sigma_1$ | loss parameters | `Material.sigma0`, `.sigma1` |
259| $z_0(x, y)$, $h(x, y)$ | mid-surface height, thickness | `CymbalProfile.z0`, `.thickness` |
260| $C_e, D_e, K_e$ | membrane stiffness, bending rigidity, DKT stiffness of element $e$ | `ShellModel.Cm_e`, `.D_e`, `.Ke` |
261| $G_x, G_y$ | CST gradient operators | `Mesh.Gx`, `.Gy` |
262| $M$, $I_e$, $f$ | lumped mass, rotational inertia, its scaling factor | `ShellModel.M`, `.I_e`, `Numerics.rot_inertia_factor` |
263| $r_k$ | rigid-body modes | `ShellModel.Rm` |
264| $\omega_{max}$, $\Delta t$, $m$ | stability bound, sub-step, sub-steps per sample | `ShellModel.omega_max`, `Simulator.dt`, `.substeps` |
265| $k, \alpha, P$ | contact stiffness, exponent, springs per patch | `ContactPatch.k_total`, `.alpha`, `.n_points` |
266| $z_s, v_s$ | exciter face height and velocity | state entries `zs`, `vs` |
267| $p_{front}, p_{back}$ | Rayleigh front and diffracted back pressure | `SimulationResult.pressure_rayleigh`, `.pressure_unbaffled` |
e30c940save .cymsim fileJeremy Magland 268| $n_i$, $a_i$ | nodal normal of the rest surface, normal acceleration | `ShellModel.normal_node`, `SimulationResult.acceleration` |
269| $h(d)$, $W$ | decimation kernel, its span in audio samples | `decimation.decimation_kernel`, `DECIMATION_SPAN` |