1793cecTime-domain nonlinear shell cymbal after Samejima (2021)Jeremy Magland 1function S = cymbal2_system(P)
2% CYMBAL2_SYSTEM Precompute everything the time loop needs. The linear
3% updates block-diagonalize over azimuthal wavenumber, so the implicit
4% shell solve and the Airy (stress function) solve are stored as one small
5% radial matrix per transform column: the time step becomes
6% Uhat+ = (A^-1 B) Uhat + (A^-1 C) Uhat- + A^-1 (ht^2/rho h) NLhat
7% with A, B, C the matrices of Samejima (2021) Eq. (20) restricted to one
8% wavenumber, the boundary-condition rows of Eqs. (22)-(30) built in, and
9% NL the explicit nonlinear + forcing terms evaluated in physical space.
10% Rows are equilibrated before inversion; the stored inverses absorb the
11% scaling so the loop applies them to raw right-hand sides.
13G = cymbal2_grid(P);
14nrt = G.nrt;
15Nr = G.Nr;
16nth = G.nth;
18h = P.thickness_m;
19rhoh = P.density * h;
20Eh = P.youngs_modulus * h;
21Rs = P.radius_m^2 / (2 * P.rise_m);
22ht = 1 / (P.sample_rate * P.oversample);
23g = P.g_param;
24phiw = P.phi_param;
25kap = 0; % set below from cymbal2_perm (needs Dflex)
26s1 = P.sigma1 * P.radius_m^2;
27s0 = P.sigma0;
29MW = P.washer_mass;
30RW = P.washer_damping;
31KW = P.washer_stiffness;
32% Spring term split as phiw/2 on the (k+1, k-1) average and (1 - phiw) on
33% step k. Note the paper prints the complementary split, which is unstable
34% for a massless washer (the roots of (1-phi)/2 (z^2+1) + phi z include
35% |z| > 1 whenever phi > 1/2); with the implicit weight on the average,
36% phi = 0.7 is unconditionally stable.
37AW = MW + (ht / 2) * RW + (phiw / 2) * ht^2 * KW;
38BW = 2 * MW - (1 - phiw) * ht^2 * KW;
39CW = -MW + (ht / 2) * RW - (phiw / 2) * ht^2 * KW;
41AiB = zeros(nrt, nrt, nth);
42AiC = zeros(nrt, nrt, nth);
43Ai = zeros(nrt, nrt, nth);
44Fi = zeros(nrt, nrt, nth);
45nDropped = 0;
47e1 = zeros(1, nrt); e1(1) = 1;
48e3 = zeros(1, nrt); e3(3) = 1;
49I = eye(nrt);
51for m = 0:nth / 2
52 lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;
53 op = cymbal2_perm(P, G, lam2);
54 kap = op.Dflex / rhoh;
56 Am = zeros(nrt, nrt);
57 Bm = zeros(nrt, nrt);
58 Cm = zeros(nrt, nrt);
59 if P.support_clamped ~= 0
60 l0 = 4;
61 else
62 l0 = 3;
63 end
64 for l = l0:Nr
65 Am(l, :) = (1 + s0 * ht) * I(l, :) + g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);
66 Bm(l, :) = 2 * I(l, :) - (1 - 2 * g) * kap * ht^2 * op.Bih(l, :);
67 Cm(l, :) = (s0 * ht - 1) * I(l, :) - g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);
68 end
69 if P.support_clamped ~= 0
70 Am(1, :) = e1;
71 Am(2, :) = e3;
72 Am(3, :) = G.D1(3, :);
73 else
74 % The shear side of the washer balance gets the same implicit split as
75 % the spring; the paper evaluates it at step k only, which feeds a
76 % spurious mode on the virtual lines when the washer is massless.
77 Am(1, :) = op.MomW;
78 Am(2, :) = AW * e3 - (phiw / 2) * ht^2 * op.ShW;
79 Bm(2, :) = BW * e3 + (1 - phiw) * ht^2 * op.ShW;
80 Cm(2, :) = CW * e3 + (phiw / 2) * ht^2 * op.ShW;
81 end
82 Am(Nr + 1, :) = op.MomE;
83 Am(Nr + 2, :) = op.ShE;
85 % Row equilibration, absorbed into the stored inverses
86 sc = 1.0 ./ max(abs(Am), [], 2);
87 Asc = sc .* Am;
88 Aisc = inv(Asc);
89 AiBm = Aisc * (sc .* Bm);
90 AiCm = Aisc * (sc .* Cm);
91 Aim = Aisc .* sc'; % inv(Am) applied to raw rhs
93 % On grids too coarse to resolve a wavenumber's inner boundary layer,
94 % the discretized boundary conditions can hand that wavenumber a
95 % negative stiffness; such columns are truncated rather than integrated.
96 % At the default grid none are.
97 sr = max(abs(eig([AiBm, AiCm; I, zeros(nrt)])));
98 if sr > 1 + 1e-7
99 AiBm = zeros(nrt, nrt);
100 AiCm = zeros(nrt, nrt);
101 Aim = zeros(nrt, nrt);
102 nDropped = nDropped + 1 + (m > 0 && m < nth / 2);
103 end
105 % Airy solve: biharmonic rows inside, phi = phi_r = 0 at both circles
106 Fm = zeros(nrt, nrt);
107 Fm(1, 1) = 1;
108 Fm(2, :) = G.D1(3, :);
109 Fm(3, 3) = 1;
110 for l = 4:Nr - 1
111 Fm(l, :) = op.Bih(l, :);
112 end
113 Fm(Nr, Nr) = 1;
114 Fm(Nr + 1, :) = G.D1(Nr, :);
115 Fm(Nr + 2, Nr + 2) = 1;
116 scf = 1.0 ./ max(abs(Fm), [], 2);
117 Fim = inv(scf .* Fm) .* scf';
119 if m == 0
120 cols = 1;
121 elseif m == nth / 2
122 cols = nth;
123 else
124 cols = [2 * m, 2 * m + 1];
125 end
126 for cc = 1:length(cols)
127 AiB(:, :, cols(cc)) = AiBm;
128 AiC(:, :, cols(cc)) = AiCm;
129 Ai(:, :, cols(cc)) = Aim;
130 Fi(:, :, cols(cc)) = Fim;
131 end
132end
134if nDropped > nth / 2
135 error(['cymbal2: the radial grid cannot resolve the inner boundary layer ' ...
136 '(need roughly washer_radius > 3 (a - c) / (n_radial - 3)); raise n_radial']);
137end
139% Strike and pickup nodes (radial fraction 0 at the washer, 1 at the rim)
140ls = 3 + round(min(1, max(0, P.strike_r)) * (Nr - 3));
141lp = 3 + round(min(1, max(0, P.pickup_r)) * (Nr - 3));
142js = 1 + mod(round(P.strike_theta_deg / 360 * nth), nth);
143jp = 1 + mod(round(P.pickup_theta_deg / 360 * nth), nth);
145S.G = G;
146S.AiB = AiB; S.AiC = AiC; S.Ai = Ai; S.Fi = Fi;
147S.rhoh = rhoh; S.Eh = Eh; S.Rs = Rs; S.ht = ht;
148S.ls = ls; S.lp = lp; S.js = js; S.jp = jp;
149S.eomRow0 = 3 + (P.support_clamped ~= 0); % first equation-of-motion row
150S.nDropped = nDropped; % truncated theta columns
151S.strikeRadiusM = G.r(ls);
152S.pickupRadiusM = G.r(lp);
153S.forceScale = 1 / (G.r(ls) * G.hr * G.hth); % point force -> pressure
154end