1function [f, W] = cymbal2_modes(P, m)
2% CYMBAL2_MODES Eigenfrequencies (Hz) of the linearized shell for one
3% azimuthal wavenumber m, with the membrane (curvature) stiffening and the
4% support boundary condition included; the washer is taken as its static
5% stiffness (its mass and damping are ignored here, as in the paper's
6% linear FEM reference). Used by the demo and the validation, not by the
7% renderer. W returns the mode shapes on the real radial nodes l = 3..Nr.
9G = cymbal2_grid(P);
10nrt = G.nrt;
11Nr = G.Nr;
12hth = G.hth;
13lam2 = (2 * cos(m * hth) - 2) / hth^2;
14op = cymbal2_perm(P, G, lam2);
16E = P.youngs_modulus;
17h = P.thickness_m;
18rhoh = P.density * h;
19Rs = P.radius_m^2 / (2 * P.rise_m);
21% Airy solve operator: Fm * phi = Z * (E h / Rs) * Lap * u, with
22% phi = phi_r = 0 at both circles and identity rows on the outermost
23% virtual lines.
24Fm = zeros(nrt, nrt);
25Fm(1, 1) = 1;
26Fm(2, :) = G.D1(3, :);
27Fm(3, 3) = 1;
28for l = 4:Nr - 1
29 Fm(l, :) = op.Bih(l, :);
30end
31Fm(Nr, Nr) = 1;
32Fm(Nr + 1, :) = G.D1(Nr, :);
33Fm(Nr + 2, Nr + 2) = 1;
34Zl = zeros(nrt, nrt);
35for l = 4:Nr - 1
36 Zl(l, :) = op.Lap(l, :);
37end
38Phimap = Fm \ ((E * h / Rs) * Zl);
40% Full stiffness on the equation-of-motion rows l = 3..Nr
41Kfull = zeros(nrt, nrt);
42Kterm = op.Dflex * op.Bih + (1 / Rs) * op.Lap * Phimap;
43for l = 3:Nr
44 Kfull(l, :) = Kterm(l, :);
45end
47% Boundary rows, static. With the washer the shell's equation of motion
48% holds on l = 3..Nr and the four virtual lines are eliminated; with the
49% rigid clamp u and u_r vanish at l = 3, so u(3) itself joins the
50% eliminated set and the equations of motion start at l = 4 (the innermost
51% virtual line is then never engaged).
52e3 = zeros(1, nrt); e3(3) = 1;
53if P.support_clamped ~= 0
54 Cmat = [e3; G.D1(3, :); op.MomE; op.ShE];
55 keep = 4:Nr;
56 elim = [2, 3, Nr + 1, Nr + 2];
57else
58 Cmat = [op.MomW; op.ShW - P.washer_stiffness * e3; op.MomE; op.ShE];
59 keep = 3:Nr;
60 elim = [1, 2, Nr + 1, Nr + 2];
61end
62X = -Cmat(:, elim) \ Cmat(:, keep);
63Kred = Kfull(keep, keep) + Kfull(keep, elim) * X;
65[V, Dm] = eig(Kred / rhoh);
66om2 = real(diag(Dm));
67[om2, idx] = sort(om2);
68om2 = max(om2, 0);
69f = sqrt(om2) / (2 * pi);
70W = real(V(:, idx));
71end