1function [K, M] = triangle_fem(nodes, family, E, rho, nu, d)
2% TRIANGLE_FEM Free-free frame system for one mode family of a bent rod.
3%
4% The rod centerline is the open polyline `nodes` (2 x nn, metres, in the
5% x-y plane), discretized into nn - 1 two-node elements. Because the geometry
6% is planar, the six degrees of freedom per node decouple exactly into two
7% independent families of three, and each family is assembled separately:
8% family 0 (in-plane): ux, uy, theta_z (in-plane bending + stretching)
9% family 1 (out-of-plane): uz, theta_x, theta_y (out-of-plane bending + torsion)
10% The bends are where the families earn their keep: within one straight leg
11% out-of-plane bending and torsion are independent, but a change of direction
12% rotates the (torsion, bending-rotation) pair, so the corner converts one
13% into the other. That conversion is what the per-element rotation T applies,
14% and it is what makes a triangle sound different from a straight rod.
15%
16% Elements are Euler-Bernoulli (cubic Hermite bending, linear axial and
17% torsion interpolation) with consistent mass including rotary inertia.
18% Shear deformation is neglected, which leaves the highest retained modes a
19% few percent sharp for a rod this thick; see the README.
20%
21% Sign conventions: in-plane bending uses theta_z = +dv/dx, giving the
22% textbook stiffness block; out-of-plane bending uses theta_y = -dw/dx, so
23% the odd coupling terms flip sign (the similarity transform Sf below).
24%
25% Returns K and M, each 3 nn x 3 nn, symmetric; M is positive definite, so
26% chol(M) reduces the generalized eigenproblem to a symmetric standard one.
28nn = size(nodes, 2);
29ndof = 3 * nn;
30K = zeros(ndof);
31M = zeros(ndof);
33Aa = pi * d^2 / 4; % cross-section area
34Iz = pi * d^4 / 64; % bending second moment (circular)
35Jp = pi * d^4 / 32; % polar moment = torsion constant (circular)
36Gs = E / (2 * (1 + nu));
38for e = 1:nn - 1
39 dx = nodes(1, e + 1) - nodes(1, e);
40 dy = nodes(2, e + 1) - nodes(2, e);
41 l = sqrt(dx^2 + dy^2);
42 c = dx / l;
43 s = dy / l;
45 % bending block on the (transverse, rotation) pairs of both nodes
46 Kb = (E * Iz / l^3) * ...
47 [12, 6 * l, -12, 6 * l; ...
48 6 * l, 4 * l^2, -6 * l, 2 * l^2; ...
49 -12, -6 * l, 12, -6 * l; ...
50 6 * l, 2 * l^2, -6 * l, 4 * l^2];
51 Mb = (rho * Aa * l / 420) * ...
52 [156, 22 * l, 54, -13 * l; ...
53 22 * l, 4 * l^2, 13 * l, -3 * l^2; ...
54 54, 13 * l, 156, -22 * l; ...
55 -13 * l, -3 * l^2, -22 * l, 4 * l^2] ...
56 + (rho * Iz / (30 * l)) * ...
57 [36, 3 * l, -36, 3 * l; ...
58 3 * l, 4 * l^2, -3 * l, -l^2; ...
59 -36, -3 * l, 36, -3 * l; ...
60 3 * l, -l^2, -3 * l, 4 * l^2];
62 Ke = zeros(6, 6);
63 Me = zeros(6, 6);
64 if family == 0
65 % local order (u1, v1, t1, u2, v2, t2): axial on u, bending on (v, t)
66 Ka2 = (E * Aa / l) * [1, -1; -1, 1];
67 Ma2 = (rho * Aa * l / 6) * [2, 1; 1, 2];
68 ia = [1 4];
69 ib = [2 3 5 6];
70 Ke(ia, ia) = Ka2;
71 Me(ia, ia) = Ma2;
72 Ke(ib, ib) = Kb;
73 Me(ib, ib) = Mb;
74 T = [c, s, 0; -s, c, 0; 0, 0, 1];
75 else
76 % local order (w1, tx1, ty1, w2, tx2, ty2): torsion on tx (about the leg
77 % axis), bending on (w, ty) with theta_y = -dw/dx
78 Kt2 = (Gs * Jp / l) * [1, -1; -1, 1];
79 Mt2 = (rho * Jp * l / 6) * [2, 1; 1, 2];
80 Sf = diag([1, -1, 1, -1]);
81 ia = [2 5];
82 ib = [1 3 4 6];
83 Ke(ia, ia) = Kt2;
84 Me(ia, ia) = Mt2;
85 Ke(ib, ib) = Sf * Kb * Sf;
86 Me(ib, ib) = Sf * Mb * Sf;
87 T = [1, 0, 0; 0, c, s; 0, -s, c];
88 end
90 TT = blkdiag(T, T);
91 Kg = TT' * Ke * TT;
92 Mg = TT' * Me * TT;
93 dofs = [3 * e - 2, 3 * e - 1, 3 * e, 3 * e + 1, 3 * e + 2, 3 * e + 3];
94 K(dofs, dofs) = K(dofs, dofs) + Kg;
95 M(dofs, dofs) = M(dofs, dofs) + Mg;
96end
97end