1function basis = triangle_basis(P)
2% TRIANGLE_BASIS Mode set of the open triangle: FEM assembly and eigensolve.
3%
4% The planar geometry splits the frame into two exactly decoupled families
5% (see triangle_fem.m): in-plane modes (bending + stretching) and
6% out-of-plane modes (bending + torsion, converted into each other at the
7% two closed corners). Each family gives a generalized eigenproblem
8% K phi = om^2 M phi, reduced to symmetric standard form through the
9% Cholesky factor of the consistent mass matrix and diagonalized by svd
10% rather than eig: the reduced matrix is symmetric positive semidefinite,
11% so its singular vectors are its eigenvectors, and svd keeps them
12% orthonormal even for the near-degenerate pairs this geometry produces
13% (eig's nonsymmetric path returns complex, non-orthogonal vectors there).
14% The eigenvectors then come back mass-normalized: phi' M phi = I. The
15% three rigid-body modes of each family sit at zero frequency and are
16% dropped along with anything below min_mode_hz; modes above max_mode_hz
17% are truncated (the broadband part they carry is reinstated by the direct
18% term in triangle_radiate.m).
19%
20% For radiation and excitation the basis stores, per retained mode, the
21% displacement at every node (NUx, NUy, NUz; out-of-plane modes have only
22% NUz, in-plane modes only NUx, NUy) and the transverse displacement at
23% every element center (EUpx, EUpy, EUpz), with the axial component removed:
24% axial motion of a thin rod displaces almost no air. Element centers,
25% tangents, and lengths are kept alongside so the radiation stage can treat
26% each element as a point dipole.
27%
28% Damping: zeta = 0.5 (loss_factor + loss_slope om) plus radiation damping
29% from the incoherent sum of the element dipole powers,
30% zeta_rad = rho0 om^3 sum_e (Vdip_e^2 |u_perp,e|^2) / (24 pi c^3),
31% Vdip_e = 2 pi a^2 l_e (the factor 2 is the added-mass coefficient of a
32% transversally moving cylinder). The incoherent sum ignores the
33% cancellation between anti-phase sections, which matters only below the
34% coincidence frequency (about 1.9 kHz for an 8 mm steel rod); above it the
35% flexural waves are supersonic and the estimate is fair.
37geo = triangle_geometry(P);
38nn = geo.nn;
39nel = nn - 1;
40E = P.youngs_modulus;
41rho = P.density;
42nu = P.poisson_ratio;
43d = P.rod_diameter_m;
44c = P.speed_of_sound;
45rho0 = P.air_density;
46nodes = geo.nodes;
48% ---- element geometry ----
49ex = zeros(nel, 1);
50ey = zeros(nel, 1);
51tx = zeros(nel, 1);
52ty = zeros(nel, 1);
53el = zeros(nel, 1);
54for e = 1:nel
55 dx = nodes(1, e + 1) - nodes(1, e);
56 dy = nodes(2, e + 1) - nodes(2, e);
57 el(e) = sqrt(dx^2 + dy^2);
58 tx(e) = dx / el(e);
59 ty(e) = dy / el(e);
60 ex(e) = (nodes(1, e) + nodes(1, e + 1)) / 2;
61 ey(e) = (nodes(2, e) + nodes(2, e + 1)) / 2;
62end
64% ---- eigensolve per family, collect retained modes ----
65fs = [];
66fams = [];
67NUx = [];
68NUy = [];
69NUz = [];
70EUpx = [];
71EUpy = [];
72EUpz = [];
73omLo = 2 * pi * P.min_mode_hz;
74for family = 0:1
75 [K, M] = triangle_fem(nodes, family, E, rho, nu, d);
76 Rc = chol(M);
77 Am = (Rc') \ (K / Rc);
78 Am = (Am + Am') / 2;
79 [Uv, Sv, ~] = svd(Am);
80 om2 = flipud(diag(Sv)); % ascending
81 Vv = fliplr(Uv);
82 om = sqrt(max(0, om2)) * P.stiffness_trim;
83 keep = find(om > omLo & om < 2 * pi * P.max_mode_hz);
84 if length(keep) == 0
85 continue;
86 end
87 Phi = Rc \ Vv(:, keep); % mass-normalized shapes, 3 nn x nk
88 nk = length(keep);
89 if family == 0
90 ux = Phi(1:3:3 * nn, :);
91 uy = Phi(2:3:3 * nn, :);
92 uz = zeros(nn, nk);
93 eux = (ux(1:nn - 1, :) + ux(2:nn, :)) / 2;
94 euy = (uy(1:nn - 1, :) + uy(2:nn, :)) / 2;
95 ax = eux .* (tx * ones(1, nk)) + euy .* (ty * ones(1, nk));
96 upx = eux - ax .* (tx * ones(1, nk));
97 upy = euy - ax .* (ty * ones(1, nk));
98 upz = zeros(nel, nk);
99 else
100 ux = zeros(nn, nk);
101 uy = zeros(nn, nk);
102 uz = Phi(1:3:3 * nn, :);
103 upx = zeros(nel, nk);
104 upy = zeros(nel, nk);
105 upz = (uz(1:nn - 1, :) + uz(2:nn, :)) / 2;
106 end
107 fs = [fs; om(keep) / (2 * pi)];
108 fams = [fams; family * ones(nk, 1)];
109 NUx = [NUx, ux];
110 NUy = [NUy, uy];
111 NUz = [NUz, uz];
112 EUpx = [EUpx, upx];
113 EUpy = [EUpy, upy];
114 EUpz = [EUpz, upz];
115end
116n = length(fs);
117if n == 0
118 error('No rod modes between min_mode_hz and max_mode_hz -- check geometry/material.');
119end
121[fs, order] = sort(fs);
122fs = reshape(fs, n, 1);
123fams = reshape(fams(order), n, 1);
124NUx = NUx(:, order);
125NUy = NUy(:, order);
126NUz = NUz(:, order);
127EUpx = EUpx(:, order);
128EUpy = EUpy(:, order);
129EUpz = EUpz(:, order);
131% ---- damping: material + radiation ----
132omega = 2 * pi * fs;
133a = d / 2;
134Vdip = 2 * pi * a^2 * el; % nel x 1 dipole volume per element
135up2 = EUpx.^2 + EUpy.^2 + EUpz.^2; % nel x n
136S2 = (Vdip.^2).' * up2; % 1 x n
137zetaRad = rho0 * omega.^3 .* reshape(S2, n, 1) / (24 * pi * c^3);
138zeta = 0.5 * (P.loss_factor + P.loss_slope * omega) + P.radiation_loss * zetaRad;
140basis.n = n;
141basis.f = fs;
142basis.omega = omega;
143basis.zeta = zeta;
144basis.zetaRad = zetaRad;
145basis.family = fams;
146basis.nodeX = nodes(1, :);
147basis.nodeY = nodes(2, :);
148basis.sfrac = geo.sfrac;
149basis.totalLen = geo.totalLen;
150basis.NUx = NUx;
151basis.NUy = NUy;
152basis.NUz = NUz;
153basis.elemX = ex;
154basis.elemY = ey;
155basis.elemTanX = tx;
156basis.elemTanY = ty;
157basis.elemLen = el;
158basis.EUpx = EUpx;
159basis.EUpy = EUpy;
160basis.EUpz = EUpz;
161end