1793cecTime-domain nonlinear shell cymbal after Samejima (2021)Jeremy Magland 1% Validate the cymbal2 model. There is no external reference
2% implementation, so the checks are (a) closed-form limits the
3% discretization must reproduce (the free circular plate, rigid-body
4% motion on the washer, cantilever statics of the stick), (b) agreement
5% with the linear-FEM eigenfrequencies published in Samejima (2021) for
6% the modes that model resolves the same way, (c) structural identities
7% of the transform and of the assembled update, (d) exact reproducibility,
8% and (e) pinned end-to-end diagnostics from a reduced configuration,
9% recorded when the model was built and treated as a regression reference.
10% Run with: numbl run cymbal2_validate.m
12fails = 0;
13P = cymbal2_defaults();
15% ---- 1. flat limit: the completely free circular plate ----
16% With the curvature removed and the washer stiffness zeroed, the shell
17% must reproduce the free plate's frequency equation (the 2 x 2 J_n / I_n
18% determinant under zero radial moment and zero Kelvin-Kirchhoff shear).
19% m = 0 is excluded: its limit as the free hole closes is logarithmic, so
20% a finite hole shifts it far more than the m >= 1 modes.
21Q = P;
22Q.rise_m = 1e-9;
23Q.washer_stiffness = 0;
24Q.washer_radius_m = 0.010;
25Q.n_radial = 96;
26Dp = Q.youngs_modulus * Q.thickness_m^3 / (12 * (1 - Q.poisson_ratio^2));
27scal = sqrt(Dp / (Q.density * Q.thickness_m)) / Q.radius_m^2 / (2 * pi);
28errP = 0;
29for m = 1:4
30 f = cymbal2_modes(Q, m);
31 fe = f(f > 0.5);
32 lam = c2v_plate_roots(m, Q.poisson_ratio);
33 fc = lam.^2 * scal;
34 nc = min(3, min(length(fe), length(fc)));
35 errP = max(errP, max(abs(fe(1:nc) - fc(1:nc)') ./ fc(1:nc)'));
36end
37ok = errP < 0.04;
38fails = fails + (~ok);
39fprintf('flat plate : max freq err %.1e over m = 1..4 %s\n', errP, c2v_verdict(ok));
41% ---- 2. rigid limit: bounce and tilt on the washer ----
42% With the plate stiffened a hundredfold the lowest m = 0 and m = 1 modes
43% must approach the closed-form rigid-body values on the line spring:
44% bounce f^2 = KW 2 pi c / (rho h pi (a^2 - c^2)) / (2 pi)^2 and tilt
45% f^2 = 4 KW c^3 / (rho h (a^4 - c^4)) / (2 pi)^2.
46Q = P;
47Q.rise_m = 1e-9;
48Q.youngs_modulus = P.youngs_modulus * 100;
49Q.n_radial = 96;
50a = Q.radius_m; c = Q.washer_radius_m;
51rhoh = Q.density * Q.thickness_m;
52KW = Q.washer_stiffness;
53f0r = sqrt(KW * 2 * pi * c / (rhoh * pi * (a^2 - c^2))) / (2 * pi);
54f1r = sqrt(4 * KW * c^3 / (rhoh * (a^4 - c^4))) / (2 * pi);
55fa = cymbal2_modes(Q, 0); fa = fa(fa > 0.01);
56fb = cymbal2_modes(Q, 1); fb = fb(fb > 0.01);
57e0 = abs(fa(1) / f0r - 1);
58e1 = abs(fb(1) / f1r - 1);
59ok = e0 < 0.03 && e1 < 0.03;
60fails = fails + (~ok);
61fprintf('washer : bounce err %.1e, tilt err %.1e vs rigid closed form %s\n', e0, e1, c2v_verdict(ok));
63% ---- 3. the paper's linear FEM reference, where comparable ----
64% Samejima's Fig. 9 gives linear-FEM eigenfrequencies of the same shell.
65% For the clamped center, the (2,1), (3,1), (4,1) modes must agree to a
66% few percent ((1,1) is excluded: it is controlled by the support model,
67% where this shell's moment-free hole and their 2-D elements differ).
68Q = P;
69Q.support_clamped = 1;
70ref = [20.1, 43.3, 77.7];
71got = zeros(1, 3);
72for m = 2:4
73 f = cymbal2_modes(Q, m);
74 f = f(f > 1);
75 got(m - 1) = f(1);
76end
77errF = max(abs(got - ref) ./ ref);
78ok = errF < 0.06;
79fails = fails + (~ok);
80fprintf('paper FEM : clamped (2,1)(3,1)(4,1) = %.1f %.1f %.1f Hz (ref %.1f %.1f %.1f), err %.1e %s\n', ...
81 got(1), got(2), got(3), ref(1), ref(2), ref(3), errF, c2v_verdict(ok));
83% ---- 4. transform and update-structure identities ----
84% The azimuthal transform must be exactly invertible, and the physical
85% derivative pipeline must agree with the per-wavenumber operators: for a
86% field cos(m theta) g(r), the physical-space Laplacian must equal the
87% per-m radial Laplacian applied to g.
88G = cymbal2_grid(P);
89eT = max(max(abs(G.F * G.Gi - eye(G.nth))));
90gr = exp(-((G.r - 0.12) / 0.05).^2);
91m = 7;
92Uf = gr * cos(m * G.th');
93LapPhys = G.D2 * Uf + G.Rd .* (G.D1 * Uf) + G.R2d .* (Uf * G.T2);
94lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;
95op = cymbal2_perm(P, G, lam2);
96LapM = (op.Lap * gr) * cos(m * G.th');
97r2 = 2:G.nrt - 1;
98eL = max(max(abs(LapPhys(r2, :) - LapM(r2, :)))) / max(max(abs(LapM(r2, :))));
99ok = eT < 1e-10 && eL < 1e-10;
100fails = fails + (~ok);
101fprintf('operators : transform inverse %.1e, physical vs per-m Laplacian %.1e %s\n', ...
102 eT, eL, c2v_verdict(ok));
104% ---- 5. the assembled update is not amplifying ----
105% Every azimuthal column's two-step update must have spectral radius at
106% most 1 (the damping makes physical modes strictly contractive; boundary
107% constraints contribute zeros). The system builder truncates columns
108% that fail; at the defaults none may be.
109S = cymbal2_system(P);
110nrt = S.G.nrt;
111srMax = 0;
112for col = [1, 2, S.G.nth / 2, S.G.nth]
113 T = [S.AiB(:, :, col), S.AiC(:, :, col); eye(nrt), zeros(nrt)];
114 srMax = max(srMax, max(abs(eig(T))));
115end
116ok = S.nDropped == 0 && srMax <= 1 + 1e-9;
117fails = fails + (~ok);
118fprintf('stability : %d truncated columns, sampled spectral radius %.9f %s\n', ...
119 S.nDropped, srMax, c2v_verdict(ok));
121% ---- 6. stick statics ----
122% A static tip load on the clamped-hand stick must deflect by
123% F L^3 / (3 E I); the rotational-spring hand adds F L^2 / K_R times L.
124MBc = cymbal2_mallet(setfield(P, 'grip_clamped', 1));
125MBs = cymbal2_mallet(P);
126EI = P.stick_youngs * P.stick_inertia_m4;
127L = P.stick_length_m;
128Fl = 1.0;
129% static balance of the update: at a fixed point, (I - A^-1 B - A^-1 C) d
130% equals A^-1 ht^2 K d, so the deflection under a constant tip load F is
131% (I - ASiB - ASiC) \ (tipcol * ht^2 * F)
132htv = MBc.ht;
133dc = ((eye(size(MBc.ASiB, 1)) - MBc.ASiB - MBc.ASiC) \ (MBc.tipcol * htv^2)) * Fl;
134ds = ((eye(size(MBs.ASiB, 1)) - MBs.ASiB - MBs.ASiC) \ (MBs.tipcol * htv^2)) * Fl;
135wc = dc(MBc.tipi);
136ws = ds(MBs.tipi);
137wcRef = Fl * L^3 / (3 * EI);
138wsRef = wcRef + (Fl * L / P.grip_spring) * L; % grip rotation F L / K_R times the arm
139ec = abs(wc / wcRef - 1);
140es = abs(ws / wsRef - 1);
141ok = ec < 1e-6 && es < 1e-6;
142fails = fails + (~ok);
143fprintf('stick : clamped-hand tip %.3e (ref %.3e), spring-hand %.3e (ref %.3e) %s\n', ...
144 wc, wcRef, ws, wsRef, c2v_verdict(ok));
146% ---- 7. reproducibility and pinned diagnostics ----
147% Two identical renders of a reduced configuration must agree exactly,
148% and its diagnostics are pinned. The trajectory of a nonlinear shell is
149% sensitive, so the acoustic pins carry percent-level tolerances while
150% the contact pins are tight.
151Q = P;
152Q.sample_rate = 44100;
153Q.oversample = 1;
154Q.tail_s = 1.5;
155[y1, d1] = cymbal2_render(Q);
156[y2, d2] = cymbal2_render(Q);
157rms1 = sqrt(mean(y1.^2));
158rms2 = sqrt(mean(y2.^2));
159same = max(abs(y1 - y2)) == 0;
160ok = same ...
161 && d1.contactEvents == d2.contactEvents ...
162 && abs(d1.peakForceN - 311.14) < 4 ...
163 && abs(d1.contactVel(1) + 8.223) < 0.08 ...
164 && abs(rms1 - 0.0973) < 0.010;
165fails = fails + (~ok);
166fprintf('pinned : repro %d, peakF %.2f N (ref 311.14), cvel %.3f (ref -8.223), rms %.4f (ref 0.0973) %s\n', ...
167 same, d1.peakForceN, d1.contactVel(1), rms1, c2v_verdict(ok));
169% ---- 8. the crash decays ----
170n = length(y1);
171fs = Q.sample_rate;
172head = y1(1:round(0.2 * fs));
173tail = y1(n - round(0.2 * fs):n);
174dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
175ok = dropE < -3 && dropE > -50;
176fails = fails + (~ok);
177fprintf('decay : last/first 0.2 s %.1f dB %s\n', dropE, c2v_verdict(ok));
179if fails > 0
180 fprintf('%d CHECK(S) FAILED\n', fails);
181else
182 fprintf('all checks passed\n');
183end
185function lam = c2v_plate_roots(n, nu)
186% Roots of the free-edge frequency equation of a circular plate
187lam = [];
188prev = c2v_plate_det(n, nu, 0.30);
189for x = 0.30:0.01:20
190 cur = c2v_plate_det(n, nu, x + 0.01);
191 if prev * cur < 0
192 lo = x; hi = x + 0.01;
193 for it = 1:60
194 mid = (lo + hi) / 2;
195 if c2v_plate_det(n, nu, lo) * c2v_plate_det(n, nu, mid) < 0
196 hi = mid;
197 else
198 lo = mid;
199 end
200 end
201 lam = [lam, (lo + hi) / 2];
202 end
203 prev = cur;
204end
205end
207function d = c2v_plate_det(n, nu, x)
208J = besselj(n, x); I = besseli(n, x);
209Jp = besselj(n - 1, x) - (n / x) * J;
210Ip = besseli(n - 1, x) - (n / x) * I;
211Jpp = -Jp / x + (n^2 / x^2 - 1) * J;
212Ipp = -Ip / x + (n^2 / x^2 + 1) * I;
213a11 = Jpp + nu * (Jp / x - n^2 * J / x^2);
214a12 = Ipp + nu * (Ip / x - n^2 * I / x^2);
215a21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);
216a22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);
217d = a11 * a22 - a12 * a21;
218end
220function s = c2v_verdict(ok)
221if ok
222 s = 'ok';
223else
224 s = 'MISMATCH';
225end
226end