4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1function [K, M] = cymbal_shell(G, n, E, nu, rho, Lref)
2% CYMBAL_SHELL Free-free thin-shell system for one circumferential order n.
3%
4% The cymbal is a shell of revolution, so its modes separate exactly in the
5% azimuth: every mode has n nodal diameters and is
6% u(s, th) = U(s) cos(n th) (meridional, along the surface toward the rim)
7% v(s, th) = V(s) sin(n th) (circumferential)
8% w(s, th) = W(s) cos(n th) (along the outward normal)
9% with a partner rotated by pi / 2n. A strike on the meridian th = 0 excites
10% only the cos member of each pair, which is why the model solves one member
11% per (n, k) and says so in the README. The 1-D energy functional below is the
12% th integral of the Love-Kirchhoff shell energy; every quadratic term carries
13% cos^2 or sin^2, so the integral contributes the same constant factor (pi, or
14% 2 pi at n = 0) to K and M and drops out of the eigenvalues. At n = 0 the V
15% field decouples and gives the axisymmetric torsional modes, which have no
16% normal motion and therefore radiate nothing.
17%
18% Strains (Kraus, with the curvatures signed so that a sphere of radius R has
19% 1 / R_s = 1 / R_theta = +1 / R for the outward normal):
20% eps_s = U' + W cs
21% eps_th = (n V + U r') / r + W ct
22% gamma = -n U / r + V' - V r' / r (engineering shear)
23% beta_s = U cs - W' (rotation of the normal)
24% beta_th = n W / r + V ct
25% kap_s = beta_s' = U' cs + U cs' - W''
26% kap_th = (n beta_th + beta_s r') / r
27% kap_sth = -n beta_s / r + beta_th' - beta_th r' / r (the full twist,
28% factor 2 included)
29% and the energy densities are the standard C = E t / (1 - nu^2) membrane and
30% D = E t^3 / (12 (1 - nu^2)) bending forms, integrated with weight r ds.
31%
32% Elements have two nodes and four degrees of freedom each: U and V are linear
33% (C0), W is cubic Hermite (C1) with dW/ds as its second freedom, so beta_s and
34% its derivative are both well defined. Four-point Gauss quadrature per
35% element evaluates the tabulated geometry at the quadrature points rather
36% than freezing it per element, so a curved meridian is integrated rather than
37% approximated by conical frusta.
38%
39% Rotary inertia rho t^3 / 12 (beta_s^2 + beta_th^2) is kept. It shifts the
40% retained frequencies by well under a percent, but it is what gives the
41% dW/ds freedoms mass and so makes M positive definite, which is what lets the
42% eigensolve reduce through chol(M). The dW/ds freedoms are then scaled by
43% Lref (a congruence, so the eigenvalues are untouched) to bring them onto the
44% same numerical scale as the displacements.
45%
46% Returns K and M, each 4 nn x 4 nn, symmetric, M positive definite. The
47% returned system is in the SCALED freedoms; multiply an eigenvector by
48% diag(1, 1, 1, Lref) per node to recover physical dW/ds.
50nn = G.nn;
51ndof = 4 * nn;
52K = zeros(ndof);
53M = zeros(ndof);
55% four-point Gauss on [0, 1]
56gx = [0.5 - 0.5 * 0.8611363115940526, 0.5 - 0.5 * 0.3399810435848563, ...
57 0.5 + 0.5 * 0.3399810435848563, 0.5 + 0.5 * 0.8611363115940526];
58gw = [0.3478548451374538, 0.6521451548625461, ...
59 0.6521451548625461, 0.3478548451374538] / 2;
61for e = 1:nn - 1
62 s0 = G.snode(e);
63 s1 = G.snode(e + 1);
64 L = s1 - s0;
65 Ke = zeros(8);
66 Me = zeros(8);
67 for g = 1:4
68 xi = gx(g);
69 sg = s0 + xi * L;
70 rr = interp1(G.s, G.r, sg);
71 tt = interp1(G.s, G.thk, sg);
72 dr = interp1(G.s, G.drds, sg);
73 cs = interp1(G.s, G.cs, sg);
74 ct = interp1(G.s, G.ct, sg);
75 csp = interp1(G.s, G.csp, sg);
76 ctp = interp1(G.s, G.ctp, sg);
78 % shape functions, DOF order [U1 V1 W1 W1' U2 V2 W2 W2']
79 NU = [1 - xi, 0, 0, 0, xi, 0, 0, 0];
80 NUp = [-1 / L, 0, 0, 0, 1 / L, 0, 0, 0];
81 NV = [0, 1 - xi, 0, 0, 0, xi, 0, 0];
82 NVp = [0, -1 / L, 0, 0, 0, 1 / L, 0, 0];
83 h1 = 1 - 3 * xi^2 + 2 * xi^3;
84 h2 = xi - 2 * xi^2 + xi^3;
85 h3 = 3 * xi^2 - 2 * xi^3;
86 h4 = -xi^2 + xi^3;
87 d1 = -6 * xi + 6 * xi^2;
88 d2 = 1 - 4 * xi + 3 * xi^2;
89 d3 = 6 * xi - 6 * xi^2;
90 d4 = -2 * xi + 3 * xi^2;
91 e1 = -6 + 12 * xi;
92 e2 = -4 + 6 * xi;
93 e3 = 6 - 12 * xi;
94 e4 = -2 + 6 * xi;
95 NW = [0, 0, h1, L * h2, 0, 0, h3, L * h4];
96 NWp = [0, 0, d1 / L, d2, 0, 0, d3 / L, d4];
97 NWpp = [0, 0, e1 / L^2, e2 / L, 0, 0, e3 / L^2, e4 / L];
99 Bes = NUp + cs * NW;
100 Bet = (n * NV + dr * NU) / rr + ct * NW;
101 Bg = -n * NU / rr + NVp - dr * NV / rr;
102 Bbs = cs * NU - NWp;
103 Bbt = n * NW / rr + ct * NV;
104 Bks = csp * NU + cs * NUp - NWpp;
105 Bkt = (n * Bbt + dr * Bbs) / rr;
106 Bbtp = n * (NWp / rr - dr * NW / rr^2) + ctp * NV + ct * NVp;
107 Bkst = -n * Bbs / rr + Bbtp - dr * Bbt / rr;
109 Cm = E * tt / (1 - nu^2);
110 Db = E * tt^3 / (12 * (1 - nu^2));
111 wq = gw(g) * L * rr;
113 Ke = Ke + wq * (Cm * (Bes' * Bes + Bet' * Bet ...
114 + nu * (Bes' * Bet + Bet' * Bes) ...
115 + (1 - nu) / 2 * (Bg' * Bg)) ...
116 + Db * (Bks' * Bks + Bkt' * Bkt ...
117 + nu * (Bks' * Bkt + Bkt' * Bks) ...
118 + (1 - nu) / 2 * (Bkst' * Bkst)));
119 Me = Me + wq * (rho * tt * (NU' * NU + NV' * NV + NW' * NW) ...
120 + rho * tt^3 / 12 * (Bbs' * Bbs + Bbt' * Bbt));
121 end
122 d = [4 * e - 3, 4 * e - 2, 4 * e - 1, 4 * e, 4 * e + 1, 4 * e + 2, 4 * e + 3, 4 * e + 4];
123 K(d, d) = K(d, d) + Ke;
124 M(d, d) = M(d, d) + Me;
125end
127% scale the dW/ds freedoms so every freedom is numerically comparable
128sc = ones(ndof, 1);
129sc(4:4:ndof) = Lref;
130Sd = sc * sc';
131K = (K .* Sd + (K .* Sd)') / 2;
132M = (M .* Sd + (M .* Sd)') / 2;
133end