1// The float64 reference discretization: symmetry and rigid modes of the
2// linear stiffness, forces as the gradient of the energy, dissipativity,
3// thickness scaling, the stand, and the free-plate frequencies.
4import { test } from "node:test";
5import assert from "node:assert/strict";
6import { makeGrid, withMask, type Grid, type Plate, nodeMasses, standRates, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
7import { uniformPlate, accelerations, addDamping, addPadMemory, padState, stepPadMemory, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
8import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
9import { fullLoss } from "../src/model/reference.ts";
10import { feltModulus } from "../src/felt.ts";
12const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
14function rnd(seed: number): () => number {
15 let s = seed >>> 0;
16 return () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296 - 0.5; };
17}
18function randomFields(pl: Plate, amp: number, seed: number): Fields {
19 const g = pl.grid, r = rnd(seed), f = fields(g.n);
20 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
21 if (!insideNode(g, i, j)) continue;
22 const k = j * g.nx + i;
23 f.u[k] = 0.02 * amp * r(); f.v[k] = 0.02 * amp * r(); f.w[k] = amp * r();
24 }
25 return f;
26}
27function bowed(g: Grid, sag: number): Float64Array {
28 const w0 = new Float64Array(g.n);
29 const Lx = (g.nx - 1) * g.hx;
30 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = sag * (1 - (2 * nodeX(g, i) / Lx) ** 2) + 0.3 * sag * Math.sin(3 * nodeY(g, j) / Lx);
31 return w0;
32}
33function withStand(pl: Plate): Plate {
34 return { ...pl, stand: { x: 0, y: 0, pad: { innerRadius: 0.005, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
35}
37test("the linear stiffness is symmetric with exactly three rigid modes per block (uniform, tapered, disc with hole)", () => {
38 const cases: { name: string; pl: Plate; fields: 1 | 3 }[] = [];
39 const g1 = makeGrid(0.12, 0.09, 0.015);
40 cases.push({ name: "uniform rectangle", pl: uniformPlate(g1, steel, 1e-3), fields: 3 });
41 const tapered = uniformPlate(g1, steel, 1e-3);
42 for (let k = 0; k < g1.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.5 * (k % g1.nx) / g1.nx);
43 cases.push({ name: "tapered rectangle", pl: tapered, fields: 3 });
44 const gd = withMask(makeGrid(0.16, 0.16, 0.01), (x, y) => x * x + y * y <= 0.08 ** 2 && x * x + y * y >= 0.015 ** 2);
45 cases.push({ name: "disc with hole", pl: uniformPlate(gd, steel, 1e-3), fields: 1 });
46 for (const c of cases) {
47 const sys = assemble(c.pl, c.fields);
48 // Symmetry is built into the banded storage; check the rigid modes: K x = 0 for translations and rotations.
49 const g = c.pl.grid;
50 const nullVectors: Float64Array[] = [];
51 const mk = (fn: (i: number, j: number, comp: number) => number) => {
52 const x = new Float64Array(sys.n);
53 for (let d = 0; d < sys.n; d++) {
54 const k = sys.dofNode[d];
55 const comp = sys.fields === 1 ? 2 : d % 3;
56 x[d] = fn(k % g.nx, Math.floor(k / g.nx), comp);
57 }
58 return x;
59 };
60 nullVectors.push(mk((_i, _j, c) => (c === 2 ? 1 : 0)));
61 nullVectors.push(mk((i, _j, c) => (c === 2 ? nodeX(g, i) : 0)));
62 nullVectors.push(mk((_i, j, c) => (c === 2 ? nodeY(g, j) : 0)));
63 if (c.fields === 3) {
64 nullVectors.push(mk((_i, _j, c) => (c === 0 ? 1 : 0)));
65 nullVectors.push(mk((_i, _j, c) => (c === 1 ? 1 : 0)));
66 nullVectors.push(mk((i, j, c) => (c === 0 ? -nodeY(g, j) : c === 1 ? nodeX(g, i) : 0)));
67 }
68 const y = new Float64Array(sys.n);
69 for (const x of nullVectors) {
70 const { stiffnessMatvec } = require_assemble();
71 stiffnessMatvec(sys, x, y);
72 let ny = 0, nx = 0;
73 for (let d = 0; d < sys.n; d++) { ny += y[d] * y[d]; nx += x[d] * x[d]; }
74 const scale = Math.abs(stiffnessEntry(sys, 0, 0)) * Math.sqrt(nx);
75 assert.ok(Math.sqrt(ny) < 1e-9 * scale, `${c.name}: rigid mode residual ${Math.sqrt(ny) / scale}`);
76 }
77 // A non-rigid vector must not be in the null space.
78 const z = mk((i, j, c) => (c === 2 ? Math.sin(i * 0.9) * Math.cos(j * 1.1) : 0));
79 const { stiffnessMatvec } = require_assemble();
80 stiffnessMatvec(sys, z, y);
81 assert.ok(y.some((v) => Math.abs(v) > 1e-3), `${c.name}: bending is not trivial`);
82 }
83});
85import * as assembleModule from "../src/model/assemble.ts";
86function require_assemble(): typeof assembleModule { return assembleModule; }
88test("forces are the exact gradient of the discrete energy (flat and bowed, tapered, free and on a stand)", () => {
89 const g = makeGrid(0.1, 0.08, 0.01);
90 const variants: { name: string; pl: Plate }[] = [];
91 variants.push({ name: "flat free", pl: uniformPlate(g, steel, 1e-3) });
92 variants.push({ name: "bowed free", pl: uniformPlate(g, steel, 1e-3, bowed(g, 3e-3)) });
93 const tapered = uniformPlate(g, steel, 1e-3, bowed(g, 2e-3));
94 for (let k = 0; k < g.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.4 * (k % g.nx) / g.nx);
95 variants.push({ name: "bowed tapered", pl: tapered });
96 variants.push({ name: "bowed on a stand", pl: withStand(uniformPlate(g, steel, 1e-3, bowed(g, 3e-3))) });
97 for (const v of variants) {
98 const pl = v.pl;
99 const sr = standRates(pl);
100 const f = randomFields(pl, 2.5e-3, 7);
101 const masses = nodeMasses(pl);
102 const acc = fields(g.n);
103 accelerations(pl, f, acc, false, sr);
104 const dir = randomFields(pl, 1, 11);
105 // Directional derivative of U along dir versus -F . dir.
106 const eps = 1e-7;
107 const shift = (s: number): Fields => ({ u: f.u.map((x, k) => x + s * dir.u[k]), v: f.v.map((x, k) => x + s * dir.v[k]), w: f.w.map((x, k) => x + s * dir.w[k]) });
108 const dU = (potentialEnergy(pl, shift(eps), sr) - potentialEnergy(pl, shift(-eps), sr)) / (2 * eps);
109 let fd = 0;
110 for (let k = 0; k < g.n; k++) fd -= masses[k] * (acc.u[k] * dir.u[k] + acc.v[k] * dir.v[k] + acc.w[k] * dir.w[k]);
111 assert.ok(Math.abs(dU - fd) < 1e-6 * Math.abs(dU), `${v.name}: dU ${dU} vs -F.dir ${fd}`);
112 }
113});
115test("the loss terms only remove energy", () => {
116 const g = makeGrid(0.1, 0.08, 0.01);
117 const pl = uniformPlate(g, steel, 1e-3);
118 const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, 20);
119 const masses = nodeMasses(pl);
120 for (const seed of [1, 2, 3]) {
121 const vel = randomFields(pl, 1, seed);
122 const acc = fields(g.n);
123 addDamping(pl, loss, vel, acc);
124 let power = 0;
125 for (let k = 0; k < g.n; k++) power += masses[k] * (acc.u[k] * vel.u[k] + acc.v[k] * vel.v[k] + acc.w[k] * vel.w[k]);
126 assert.ok(power < 0, `power ${power}`);
127 }
128});
130test("the felt's chain, stepped in time, dissipates the work of its loss stiffness per cycle", () => {
131 const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
132 const pl = withStand(uniformPlate(g, steel, 1e-3));
133 const sr = standRates(pl)!;
134 const masses = nodeMasses(pl);
135 // A rigid translation w = A sin(omega t) of the plate, stepped finely; the work done
136 // against the arms over a cycle, once the transient of the arms has passed, is
137 // pi L(omega) K A^2 with L the loss stiffness per unit of K (eta at the bounce).
138 for (const f of [sr.bounceHz, 2000]) {
139 const omega = 2 * Math.PI * f, A = 1e-4;
140 const steps = 400, dt = 1 / (f * steps);
141 const pad = padState(pl, sr);
142 const x = fields(g.n), dw = fields(g.n), acc = fields(g.n);
143 const setW = (t: number, out: Fields) => { for (let k = 0; k < g.n; k++) out.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? A * Math.sin(omega * t) : 0; };
144 let work = 0;
145 const cycles = 4;
146 for (let s = 0; s < cycles * steps; s++) {
147 const t0 = s * dt, t1 = (s + 1) * dt;
148 setW(t0, x); setW(t1, dw);
149 for (let k = 0; k < g.n; k++) dw.w[k] -= x.w[k];
150 acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
151 addPadMemory(pl, sr, x, pad, false, acc);
152 // Work of the plate against the arms over the step, with the force at mid-step (trapezoid).
153 let F0 = 0;
154 for (let k = 0; k < g.n; k++) F0 += masses[k] * acc.w[k];
155 stepPadMemory(pl, sr, dt, dw, x, pad, false);
156 setW(t1, x);
157 acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
158 addPadMemory(pl, sr, x, pad, false, acc);
159 let F1 = 0;
160 for (let k = 0; k < g.n; k++) F1 += masses[k] * acc.w[k];
161 if (s >= (cycles - 1) * steps) work -= 0.5 * (F0 + F1) * dw.w[sr.nodes[0]];
162 }
163 const L = feltModulus(sr.chain, omega).loss;
164 const expected = Math.PI * L * sr.totalKw * A * A;
165 assert.ok(work > 0);
166 assert.ok(Math.abs(work / expected - 1) < 0.02, `at ${f.toFixed(0)} Hz: work per cycle ${work} vs ${expected}`);
167 if (f === sr.bounceHz) assert.ok(Math.abs(L / sr.eta - 1) < 0.02, `loss stiffness at the bounce ${L} vs eta ${sr.eta}`);
168 }
169});
171test("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
172 const g = makeGrid(0.1, 0.1, 0.01);
173 const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
174 // The largest frequency is set by the in-plane waves, independent of h.
175 assert.ok(Math.abs(w1 / w2 - 1) < 1e-6, `in-plane omega_max ${w1} vs ${w2}`);
176 const { stiffnessEntry: ke } = assembleModule;
177 const s1 = assemble(uniformPlate(g, steel, 1e-3)), s2 = assemble(uniformPlate(g, steel, 2e-3));
178 const c = Math.floor(s1.n / 2);
179 // Bending stiffness scales as h^3 and the mass as h: omega^2 scales as h^2.
180 assert.ok(Math.abs(ke(s2, c, c) / ke(s1, c, c) - 8) < 1e-9);
181 assert.ok(Math.abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12);
182 assert.ok(Math.abs(bendingStiffness(steel, 2e-3) / bendingStiffness(steel, 1e-3) - 8) < 1e-12);
183});
185test("the stand resists a rigid translation with its total stiffness and a tilt with the pads' second moment", () => {
186 const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
187 const pl = withStand(uniformPlate(g, steel, 1e-3));
188 const sr = standRates(pl)!;
189 const masses = nodeMasses(pl);
190 const f = fields(g.n), acc = fields(g.n);
191 const dz = 1e-4;
192 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? dz : 0;
193 accelerations(pl, f, acc, false, sr);
194 let F = 0;
195 for (let k = 0; k < g.n; k++) F += masses[k] * acc.w[k];
196 // The relaxed felt holds kInf of the stiffness at the bounce; the arms of the chain hold the rest.
197 const kInf = sr.chain.kInf;
198 assert.ok(kInf > 0.2 && kInf < 0.8, `relaxed fraction ${kInf}`);
199 assert.ok(Math.abs(-F / dz / (kInf * sr.totalKw) - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${kInf * sr.totalKw}`);
200 // Tilt about y: w = theta x; restoring moment = theta * sum_k kw_k x_k^2 (bending of a plane is
201 // zero; the von Karman membrane strain of a rotation is theta^2 / 2, so a tiny tilt keeps it negligible).
202 const th = 1e-7;
203 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? th * nodeX(g, k % g.nx) : 0;
204 accelerations(pl, f, acc, false, sr);
205 let M = 0, I2 = 0;
206 for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += kInf * sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
207 assert.ok(Math.abs(-M / th / I2 - 1) < 1e-6, `tilt stiffness ${-M / th} vs ${I2}`);
208 // With lift-off, a small motion within the preload is the same linear spring; a large upward motion unloads the bottom pad.
209 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 1e-5 : 0;
210 const lin = fields(g.n), lo = fields(g.n);
211 accelerations(pl, f, lin, false, sr);
212 accelerations(pl, f, lo, true, sr);
213 for (let k = 0; k < g.n; k++) assert.ok(Math.abs(lin.w[k] - lo.w[k]) < 1e-9 * Math.max(1, Math.abs(lin.w[k])));
214 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 0.01 : 0;
215 accelerations(pl, f, lo, true, sr);
216 let Flo = 0;
217 for (let k = 0; k < g.n; k++) Flo += masses[k] * lo.w[k];
218 // Bottom pad unloaded: only the top pad (compressed by preload + 10 mm) pushes down, minus the rest balance.
219 const kb = (kInf * sr.totalKw) / 2;
220 const expected = -(kb * (sr.preloadTop + 0.01) - kb * sr.preloadTop) - kb * sr.preloadBottom;
221 assert.ok(Math.abs(Flo / expected - 1) < 1e-6, `lift-off force ${Flo} vs ${expected}`);
222});