Blame packages/core/src/m/plate_mount.m at main - proof-of-concept/plate-dynamics-lab

proof-of-concept / plate-dynamics-lab

80c0b6d Plate dynamics lab: plates in realizable mounts under gravity, explicit finite differences on WebGPUJeremy Magland 1 month ago1function mnt = plate_mount(P, dom)

2% PLATE_MOUNT How the plate is held, as data on the grid.

4% mnt = plate_mount(P, dom), with P from plate_setup and dom from

5% plate_domain, turns params.mount into what the solvers need: the

6% direction of gravity in the plate's coordinates, the nodes the mount holds

7% still, the springs, dashpots and one-sided contacts it puts under the

8% plate, the in-plane constraints for the self-weight stress solve

9% (plate_prestress), and the rigid-body basis the results are decomposed in

10% (plate_rigid_basis). Every mount is something that can be built on a

11% bench, and each carries an idealization that is stated in mnt.notes.

12%

13% Gravity. The plate lies in the (x, y) plane and w is along its normal z.

14% mnt.gx, mnt.gy, mnt.gz is the unit vector along which gravity pulls, in

15% these coordinates, and P.g its magnitude. The in-plane part (gx, gy) puts

16% the plate under its own weight in tension or compression (plate_prestress);

17% the normal part gz loads it transversely (plate_static).

18%

19% Mounts.

20%

21% 'thread' One thread through a hole at (x, y), of length L, from a fixed

22% hook. A plate hung from one point turns in its own plane until

23% its centre of mass is straight below the hole, which fixes

24% (gx, gy). For small out-of-plane motion the thread is a

25% transverse spring of stiffness T / L at the hole, T = m g its

26% tension (the linearized pendulum), and a torsional spring

27% G J / L about the vertical axis with J = pi d^4 / 32 from the

28% thread's diameter d and shear modulus G (mount.thread).

29% 'threads' Two threads of the same length L from hooks at the same

30% height, through holes 1 and 2. The line of the holes is then

31% horizontal, and the tensions come from statics: the plate is

32% refused when its centre of mass is not between the holes,

33% because one thread would go slack.

34% 'pads' A horizontal plate resting on pads (foam, felt or a knife

35% edge): each pad is a rectangle of half-widths (rx, ry) about

36% (x, y), with a stiffness k [N/m^3] and damping c [N s/m^3] per

37% unit area. A pad pushes only while compressed, so the plate

38% can lift off and land again.

39% 'cymbal' A horizontal plate on a stand: a rod of radius r_rod through

40% the centre (no transverse effect) between two felt washers of

41% outer radius r_felt, each of stiffness k and damping c per

42% unit area and both compressed by the wing nut's preload

43% [m]. Each felt is a one-sided contact.

44% 'vise' One edge gripped in a vise to a depth [m]: the nodes within

45% the jaws are held at w = 0 and the clamped condition follows

46% from the bending energy (plate_bending). The angle tilt [deg]

47% sets the orientation: 0 hangs straight down from the jaw, 90 is

48% a horizontal cantilever, 180 stands up on the jaw with the

49% self-weight in compression.

50% 'pinched' Held between finger and thumb over a patch of radius r about

51% (x, y): a two-sided spring k [N/m] and dashpot c [N s/m] spread

52% over the patch, the patch also fixed in the plane. tilt [deg]

53% runs from 0 (hanging below the hand) to 90 (held flat).

54%

55% Fields.

56% type the mount

57% gx, gy, gz unit direction of gravity in plate coordinates

58% held ny x nx logical: nodes held at w = 0 (vise)

59% ks, cs ny x nx: linear transverse spring [N/m] and

60% dashpot [N s/m] at each node (threads, pinch)

61% kb, cb, ob bottom contacts: stiffness [N/m], dashpot

62% [N s/m] and surface offset [m] per node; the

63% contact pushes up while ob - w > 0

64% kt, ct, ot top contacts: push down while ot + w > 0

65% ktor torsional stiffness about the vertical axis

66% through the pivot [N m/rad] (one thread)

67% fix in-plane constraints for plate_prestress: rows

68% [node, dir_x, dir_y] meaning the displacement

69% of node along (dir_x, dir_y) is zero

70% free_rotation the in-plane rotation about the first fixed

71% node is free (one thread) and must be

72% regularized, which needs the load to have no

73% moment about it

74% pivot [x y] the rigid-body basis is taken about

75% ed, eh in-plane unit vectors: down and horizontal

76% nr number of rigid-body coordinates (0 or 3)

77% rigid_names {1 x nr} names of the coordinates

78% supports struct array: name, unit, nodes (index list),

79% kind ('spring', 'bottom' or 'top': which force

80% the series reports at those nodes)

81% mass, xc, yc plate mass [kg] and centre of mass, from the

82% node weights

83% notes remarks about the idealization

85mount = P.mount;

86nx = dom.nx; ny = dom.ny; dx = dom.dx; dy = dom.dy;

87X = dom.X; Y = dom.Y;

88inside = dom.inside;

89a = dom.m * dx * dy; % node areas [m^2]

90Mi = P.mu * a; % node masses [kg]

91mass = sum(Mi(:));

92xc = sum(Mi(:) .* X(:)) / mass;

93yc = sum(Mi(:) .* Y(:)) / mass;

94zero = zeros(ny, nx);

95mnt.type = mount.type;

96mnt.held = false(ny, nx);

97mnt.ks = zero; mnt.cs = zero;

98mnt.kb = zero; mnt.cb = zero; mnt.ob = zero;

99mnt.kt = zero; mnt.ct = zero; mnt.ot = zero;

100mnt.ktor = 0;

101mnt.fix = zeros(0, 3);

102mnt.free_rotation = false;

103mnt.mass = mass; mnt.xc = xc; mnt.yc = yc;

104supports = struct('name', {}, 'unit', {}, 'nodes', {}, 'kind', {});

105notes = {};

106

107switch mount.type

108 case 'thread'

109 kA = nearest_node([mount.x mount.y], dom);

110 A = [X(kA) Y(kA)]; % the plate hangs from the node the hole snaps to

111 d = [xc yc] - A;

112 if norm(d) < 1e-3 * max(P.plate.Lx, P.plate.Ly)

113 error('plate_mount: the hole is at the centre of mass; a plate hung there has no preferred orientation');

114 end

115 ed = d / norm(d);

116 eh = [-ed(2) ed(1)];

117 mnt.gx = ed(1); mnt.gy = ed(2); mnt.gz = 0;

118 T = mass * P.g;

119 mnt.ks(kA) = T / mount.length;

120 mnt.fix = [kA 1 0; kA 0 1];

121 mnt.free_rotation = true;

122 mnt.ktor = torsion(mount);

123 mnt.pivot = [X(kA) Y(kA)];

124 mnt.ed = ed; mnt.eh = eh;

125 mnt.nr = 3;

126 mnt.rigid_names = {'sway of the hole', 'swing about the hole', 'spin about the thread'};

127 supports(end + 1) = struct('name', 'thread', 'unit', 'N', 'nodes', kA, 'kind', 'spring');

128 ang = atan2(ed(1), -ed(2)) * 180 / pi;

129 notes{end + 1} = sprintf('hung from the node (%.4g, %.4g) m; the plate turns %.3g degrees in its plane so the centre of mass hangs below the hole; tension %.3g N gives a transverse spring of %.3g N/m at the hole', X(kA), Y(kA), ang, T, T / mount.length);

130 if mnt.ktor > 0

131 notes{end + 1} = sprintf('spin about the thread is restrained only by its torsion, %.3g N m/rad', mnt.ktor);

132 else

133 notes{end + 1} = 'spin about the thread is free: nothing but the sigma0 loss slows it';

134 end

135 case 'threads'

136 A1 = [mount.x1 mount.y1]; A2 = [mount.x2 mount.y2];

137 k1 = nearest_node(A1, dom); k2 = nearest_node(A2, dom);

138 if k1 == k2

139 error('plate_mount: the two holes fall on the same grid node; move them apart or refine the grid');

140 end

141 A1 = [X(k1) Y(k1)]; A2 = [X(k2) Y(k2)];

142 e = (A2 - A1) / norm(A2 - A1);

143 ed = [-e(2) e(1)]; % perpendicular to the line of the holes

144 if dot([xc yc] - A1, ed) < 0

145 ed = -ed; % the plate hangs below its threads

146 end

147 eh = [-ed(2) ed(1)];

148 s2 = norm(A2 - A1);

149 sc = dot([xc yc] - A1, eh) / dot(A2 - A1, eh); % centre of mass along the line, 0 at hole 1, 1 at hole 2

150 if ~(sc > 0 && sc < 1)

151 error('plate_mount: the centre of mass is not between the two holes, so one thread would go slack');

152 end

153 T1 = mass * P.g * (1 - sc);

154 T2 = mass * P.g * sc;

155 mnt.gx = ed(1); mnt.gy = ed(2); mnt.gz = 0;

156 mnt.ks(k1) = T1 / mount.length;

157 mnt.ks(k2) = T2 / mount.length;

158 mnt.fix = [k1 ed; k2 ed; k1 eh]; % threads inextensible, sway fixed at one hole

159 mnt.pivot = (A1 + A2) / 2;

160 mnt.ed = ed; mnt.eh = eh;

161 mnt.nr = 3;

162 mnt.rigid_names = {'sway of the threads', 'swing about the holes', 'yaw about the vertical'};

163 supports(end + 1) = struct('name', 'thread 1', 'unit', 'N', 'nodes', k1, 'kind', 'spring');

164 supports(end + 1) = struct('name', 'thread 2', 'unit', 'N', 'nodes', k2, 'kind', 'spring');

165 notes{end + 1} = sprintf('hung from the nodes (%.4g, %.4g) and (%.4g, %.4g) m on threads of %.3g m; tensions %.3g and %.3g N; the line of the holes is horizontal', A1(1), A1(2), A2(1), A2(2), mount.length, T1, T2);

166 notes{end + 1} = sprintf('span %.3g m between the holes', s2);

167 case 'pads'

168 mnt.gx = 0; mnt.gy = 0; mnt.gz = -1;

169 if ~isfield(mount, 'pads') || isempty(mount.pads)

170 error('plate_mount: the pads mount needs at least one pad');

171 end

172 npads = numel(mount.pads);

173 for k = 1:npads

174 pd = mount.pads(k);

175 in = inside & abs(X - pd.x) <= pd.rx + 1e-12 & abs(Y - pd.y) <= pd.ry + 1e-12;

176 if ~any(in(:))

177 kn = nearest_node([pd.x pd.y], dom);

178 in = false(ny, nx); in(kn) = true;

179 notes{end + 1} = sprintf('pad %d is smaller than a grid cell, so it holds the one node nearest its centre', k);

180 end

181 mnt.kb(in) = mnt.kb(in) + mount.k * a(in);

182 mnt.cb(in) = mnt.cb(in) + mount.c * a(in);

183 supports(end + 1) = struct('name', sprintf('pad %d', k), 'unit', 'N', 'nodes', find(in), 'kind', 'bottom');

184 end

185 mnt.pivot = [xc yc];

186 mnt.ed = [1 0]; mnt.eh = [0 1];

187 mnt.nr = 3;

188 mnt.rigid_names = {'heave', 'tilt along x', 'tilt along y'};

189 notes{end + 1} = sprintf('resting on %d pad(s); each pushes only while compressed, so the plate can lift off', npads);

190 case 'cymbal'

191 mnt.gx = 0; mnt.gy = 0; mnt.gz = -1;

192 r = sqrt((X - P.plate.Lx / 2).^2 + (Y - P.plate.Ly / 2).^2);

193 in = inside & r <= mount.r_felt + 1e-12 & r > mount.r_rod;

194 if ~any(in(:))

195 error('plate_mount: the felt washer covers no node; enlarge r_felt or refine the grid');

196 end

197 mnt.kb(in) = mount.k * a(in); mnt.cb(in) = mount.c * a(in); mnt.ob(in) = mount.preload;

198 mnt.kt(in) = mount.k * a(in); mnt.ct(in) = mount.c * a(in); mnt.ot(in) = mount.preload;

199 supports(end + 1) = struct('name', 'lower felt', 'unit', 'N', 'nodes', find(in), 'kind', 'bottom');

200 supports(end + 1) = struct('name', 'upper felt', 'unit', 'N', 'nodes', find(in), 'kind', 'top');

201 mnt.pivot = [P.plate.Lx / 2, P.plate.Ly / 2];

202 mnt.ed = [1 0]; mnt.eh = [0 1];

203 mnt.nr = 3;

204 mnt.rigid_names = {'heave', 'tilt along x', 'tilt along y'};

205 notes{end + 1} = sprintf('on a stand: felts over %d node(s) between radii %.3g and %.3g m, each preloaded by %.3g mm; the rod and the hole it passes through are not modelled', nnz(in), mount.r_rod, mount.r_felt, mount.preload * 1e3);

206 case 'vise'

207 depth = mount.depth;

208 switch mount.edge

209 case 'x0'

210 coord = X; edir = [1 0];

211 case 'x1'

212 coord = P.plate.Lx - X; edir = [-1 0];

213 case 'y0'

214 coord = Y; edir = [0 1];

215 case 'y1'

216 coord = P.plate.Ly - Y; edir = [0 -1];

217 otherwise

218 error('plate_mount: unknown vise edge "%s" (expected x0, x1, y0 or y1)', mount.edge);

219 end

220 held = inside & coord <= depth + 1e-12;

221 % the clamped condition needs a held node behind every held edge node

222 h = min(dx, dy);

223 if nnz(held) < 2 || max(coord(held)) < h - 1e-12

224 held = inside & coord <= h + 1e-12;

225 notes{end + 1} = sprintf('jaw depth widened to %.3g mm, two rows of nodes, which the clamped condition needs', h * 1e3);

226 end

227 if ~any(held(:))

228 error('plate_mount: the vise grips no node of the plate');

229 end

230 if nnz(inside & ~held) < 9

231 error('plate_mount: the vise leaves almost nothing of the plate free');

232 end

233 mnt.held = held;

234 tilt = mount.tilt * pi / 180;

235 mnt.gx = cos(tilt) * edir(1); mnt.gy = cos(tilt) * edir(2); mnt.gz = -sin(tilt);

236 kh = find(held);

237 mnt.fix = [kh ones(numel(kh), 1) zeros(numel(kh), 1); kh zeros(numel(kh), 1) ones(numel(kh), 1)];

238 mnt.pivot = [xc yc];

239 mnt.ed = edir; mnt.eh = [-edir(2) edir(1)];

240 mnt.nr = 0;

241 mnt.rigid_names = {};

242 notes{end + 1} = sprintf('%d nodes within %.3g mm of the %s edge are held in the vise; tilt %g degrees', nnz(held), depth * 1e3, mount.edge, mount.tilt);

243 if mount.tilt > 90

244 notes{end + 1} = 'standing on the jaw: the self-weight is compressive, which softens the lowest mode and can buckle a long thin plate';

245 end

246 case 'pinched'

247 A = [mount.x mount.y];

248 r = sqrt((X - A(1)).^2 + (Y - A(2)).^2);

249 in = inside & r <= mount.radius + 1e-12;

250 if ~any(in(:))

251 kA = nearest_node(A, dom);

252 in = false(ny, nx); in(kA) = true;

253 notes{end + 1} = 'the pinch patch is smaller than a grid cell, so it holds one node';

254 end

255 atot = sum(a(in));

256 mnt.ks(in) = mount.k * a(in) / atot;

257 mnt.cs(in) = mount.c * a(in) / atot;

258 kin = find(in);

259 mnt.fix = [kin ones(numel(kin), 1) zeros(numel(kin), 1); kin zeros(numel(kin), 1) ones(numel(kin), 1)];

260 d = [xc yc] - A;

261 if norm(d) < 1e-9

262 d = [0 -1];

263 end

264 ed = d / norm(d);

265 tilt = mount.tilt * pi / 180;

266 mnt.gx = cos(tilt) * ed(1); mnt.gy = cos(tilt) * ed(2); mnt.gz = -sin(tilt);

267 mnt.pivot = A;

268 mnt.ed = ed; mnt.eh = [-ed(2) ed(1)];

269 mnt.nr = 3;

270 mnt.rigid_names = {'heave at the hand', 'swing about the hand', 'turn about the hand'};

271 supports(end + 1) = struct('name', 'fingers', 'unit', 'N', 'nodes', kin, 'kind', 'spring');

272 notes{end + 1} = sprintf('pinched over %d node(s) about (%.4g, %.4g) m: a two-sided spring of %.3g N/m with %.3g N s/m of damping, tilt %g degrees', nnz(in), A(1), A(2), mount.k, mount.c, mount.tilt);

273end

274mnt.supports = supports;

275mnt.notes = notes;

276end

277

278function k = nearest_node(A, dom)

279d2 = (dom.X - A(1)).^2 + (dom.Y - A(2)).^2;

280d2(~dom.inside) = Inf;

281[dmin, k] = min(d2(:));

282if ~isfinite(dmin)

283 error('plate_mount: no node of the plate near (%.4g, %.4g) m', A(1), A(2));

284end