1function op = cymbal2_perm(P, G, lam2)
2% CYMBAL2_PERM Radial operators for one azimuthal wavenumber, dimensional.
3% lam2 is the (negative) eigenvalue of the circulant second-difference in
4% theta for this wavenumber, so every d^2/dtheta^2 becomes multiplication
5% by lam2. Returns the Laplacian, the biharmonic (as the exact discrete
6% composition Lap * Lap), and the radial-moment and Kelvin-Kirchhoff shear
7% rows used as boundary conditions at the washer (l = 3) and at the free
8% edge (l = Nr). The moment and shear expressions are those of Samejima
9% (2021), Eqs. (6)-(7) and (14)-(15).
11nrt = G.nrt;
12nu = P.poisson_ratio;
13Dfl = P.youngs_modulus * P.thickness_m^3 / (12 * (1 - nu^2));
15Rm = diag(G.Rd);
16R2 = diag(G.R2d);
17R3 = diag(G.Rd.^3);
19Lap = G.D2 + Rm * G.D1 + lam2 * R2;
20Bih = Lap * Lap;
22% Radial bending moment (per unit length, divided by -D):
23% u_rr + nu * (u_r / r + u_thth / r^2)
24Mom = G.D2 + nu * (Rm * G.D1 + lam2 * R2);
26% Kelvin-Kirchhoff shear: V = -D * [u_rrr + u_rr / r - u_r / r^2
27% + (nu - 3) / r^3 * u_thth + (2 - nu) / r^2 * u_r_thth]
28Sh = -Dfl * (G.D3 + Rm * G.D2 - R2 * G.D1 ...
29 + (nu - 3) * lam2 * R3 + (2 - nu) * lam2 * R2 * G.D1);
31op.Lap = Lap;
32op.Bih = Bih;
33op.MomW = Mom(3, :); % moment row at the washer circle
34op.MomE = Mom(nrt - 2, :); % moment row at the free edge (l = Nr)
35op.ShW = Sh(3, :);
36op.ShE = Sh(nrt - 2, :);
37op.Dflex = Dfl;
38end