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, 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";
11const steel = { E: 200e9, nu: 0.3, rho: 7850 };
13function rnd(seed: number): () => number {
14 let s = seed >>> 0;
15 return () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296 - 0.5; };
16}
17function randomFields(pl: Plate, amp: number, seed: number): Fields {
18 const g = pl.grid, r = rnd(seed), f = fields(g.n);
19 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
20 if (!insideNode(g, i, j)) continue;
21 const k = j * g.nx + i;
22 f.u[k] = 0.02 * amp * r(); f.v[k] = 0.02 * amp * r(); f.w[k] = amp * r();
23 }
24 return f;
25}
26function bowed(g: Grid, sag: number): Float64Array {
27 const w0 = new Float64Array(g.n);
28 const Lx = (g.nx - 1) * g.hx;
29 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);
30 return w0;
31}
32function withStand(pl: Plate): Plate {
33 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 } };
34}
36test("the linear stiffness is symmetric with exactly three rigid modes per block (uniform, tapered, disc with hole)", () => {
37 const cases: { name: string; pl: Plate; fields: 1 | 3 }[] = [];
38 const g1 = makeGrid(0.12, 0.09, 0.015);
39 cases.push({ name: "uniform rectangle", pl: uniformPlate(g1, steel, 1e-3), fields: 3 });
40 const tapered = uniformPlate(g1, steel, 1e-3);
41 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);
42 cases.push({ name: "tapered rectangle", pl: tapered, fields: 3 });
43 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);
44 cases.push({ name: "disc with hole", pl: uniformPlate(gd, steel, 1e-3), fields: 1 });
45 for (const c of cases) {
46 const sys = assemble(c.pl, c.fields);
47 // Symmetry is built into the banded storage; check the rigid modes: K x = 0 for translations and rotations.
48 const g = c.pl.grid;
49 const nullVectors: Float64Array[] = [];
50 const mk = (fn: (i: number, j: number, comp: number) => number) => {
51 const x = new Float64Array(sys.n);
52 for (let d = 0; d < sys.n; d++) {
53 const k = sys.dofNode[d];
54 const comp = sys.fields === 1 ? 2 : d % 3;
55 x[d] = fn(k % g.nx, Math.floor(k / g.nx), comp);
56 }
57 return x;
58 };
59 nullVectors.push(mk((_i, _j, c) => (c === 2 ? 1 : 0)));
60 nullVectors.push(mk((i, _j, c) => (c === 2 ? nodeX(g, i) : 0)));
61 nullVectors.push(mk((_i, j, c) => (c === 2 ? nodeY(g, j) : 0)));
62 if (c.fields === 3) {
63 nullVectors.push(mk((_i, _j, c) => (c === 0 ? 1 : 0)));
64 nullVectors.push(mk((_i, _j, c) => (c === 1 ? 1 : 0)));
65 nullVectors.push(mk((i, j, c) => (c === 0 ? -nodeY(g, j) : c === 1 ? nodeX(g, i) : 0)));
66 }
67 const y = new Float64Array(sys.n);
68 for (const x of nullVectors) {
69 const { stiffnessMatvec } = require_assemble();
70 stiffnessMatvec(sys, x, y);
71 let ny = 0, nx = 0;
72 for (let d = 0; d < sys.n; d++) { ny += y[d] * y[d]; nx += x[d] * x[d]; }
73 const scale = Math.abs(stiffnessEntry(sys, 0, 0)) * Math.sqrt(nx);
74 assert.ok(Math.sqrt(ny) < 1e-9 * scale, `${c.name}: rigid mode residual ${Math.sqrt(ny) / scale}`);
75 }
76 // A non-rigid vector must not be in the null space.
77 const z = mk((i, j, c) => (c === 2 ? Math.sin(i * 0.9) * Math.cos(j * 1.1) : 0));
78 const { stiffnessMatvec } = require_assemble();
79 stiffnessMatvec(sys, z, y);
80 assert.ok(y.some((v) => Math.abs(v) > 1e-3), `${c.name}: bending is not trivial`);
81 }
82});
84import * as assembleModule from "../src/model/assemble.ts";
85function require_assemble(): typeof assembleModule { return assembleModule; }
87test("forces are the exact gradient of the discrete energy (flat and bowed, tapered, free and on a stand)", () => {
88 const g = makeGrid(0.1, 0.08, 0.01);
89 const variants: { name: string; pl: Plate }[] = [];
90 variants.push({ name: "flat free", pl: uniformPlate(g, steel, 1e-3) });
91 variants.push({ name: "bowed free", pl: uniformPlate(g, steel, 1e-3, bowed(g, 3e-3)) });
92 const tapered = uniformPlate(g, steel, 1e-3, bowed(g, 2e-3));
93 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);
94 variants.push({ name: "bowed tapered", pl: tapered });
95 variants.push({ name: "bowed on a stand", pl: withStand(uniformPlate(g, steel, 1e-3, bowed(g, 3e-3))) });
96 for (const v of variants) {
97 const pl = v.pl;
98 const sr = standRates(pl);
99 const f = randomFields(pl, 2.5e-3, 7);
100 const masses = nodeMasses(pl);
101 const acc = fields(g.n);
102 accelerations(pl, f, acc, false, sr);
103 const dir = randomFields(pl, 1, 11);
104 // Directional derivative of U along dir versus -F . dir.
105 const eps = 1e-7;
106 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]) });
107 const dU = (potentialEnergy(pl, shift(eps), sr) - potentialEnergy(pl, shift(-eps), sr)) / (2 * eps);
108 let fd = 0;
109 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]);
110 assert.ok(Math.abs(dU - fd) < 1e-6 * Math.abs(dU), `${v.name}: dU ${dU} vs -F.dir ${fd}`);
111 }
112});
114test("the loss terms and the pad dashpots only remove energy", () => {
115 const g = makeGrid(0.1, 0.08, 0.01);
116 const pl = withStand(uniformPlate(g, steel, 1e-3));
117 const sr = standRates(pl);
118 const loss = fullLoss({ sigma0: 1, sigma1: 5e-4, kappa: 1 }, sr, 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, sr);
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("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
131 const g = makeGrid(0.1, 0.1, 0.01);
132 const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
133 // The largest frequency is set by the in-plane waves, independent of h.
134 assert.ok(Math.abs(w1 / w2 - 1) < 1e-6, `in-plane omega_max ${w1} vs ${w2}`);
135 const { stiffnessEntry: ke } = assembleModule;
136 const s1 = assemble(uniformPlate(g, steel, 1e-3)), s2 = assemble(uniformPlate(g, steel, 2e-3));
137 const c = Math.floor(s1.n / 2);
138 // Bending stiffness scales as h^3 and the mass as h: omega^2 scales as h^2.
139 assert.ok(Math.abs(ke(s2, c, c) / ke(s1, c, c) - 8) < 1e-9);
140 assert.ok(Math.abs(s2.mass[c] / s1.mass[c] - 2) < 1e-12);
141 assert.ok(Math.abs(bendingStiffness(steel, 2e-3) / bendingStiffness(steel, 1e-3) - 8) < 1e-12);
142});
144test("the stand resists a rigid translation with its total stiffness and a tilt with the pads' second moment", () => {
145 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);
146 const pl = withStand(uniformPlate(g, steel, 1e-3));
147 const sr = standRates(pl)!;
148 const masses = nodeMasses(pl);
149 const f = fields(g.n), acc = fields(g.n);
150 const dz = 1e-4;
151 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? dz : 0;
152 accelerations(pl, f, acc, false, sr);
153 let F = 0;
154 for (let k = 0; k < g.n; k++) F += masses[k] * acc.w[k];
155 assert.ok(Math.abs(-F / dz / sr.totalKw - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${sr.totalKw}`);
156 // Tilt about y: w = theta x; restoring moment = theta * sum_k kw_k x_k^2 (bending of a plane is
157 // zero; the von Karman membrane strain of a rotation is theta^2 / 2, so a tiny tilt keeps it negligible).
158 const th = 1e-7;
159 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;
160 accelerations(pl, f, acc, false, sr);
161 let M = 0, I2 = 0;
162 for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
163 assert.ok(Math.abs(-M / th / I2 - 1) < 1e-6, `tilt stiffness ${-M / th} vs ${I2}`);
164 // With lift-off, a small motion within the preload is the same linear spring; a large upward motion unloads the bottom pad.
165 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 1e-5 : 0;
166 const lin = fields(g.n), lo = fields(g.n);
167 accelerations(pl, f, lin, false, sr);
168 accelerations(pl, f, lo, true, sr);
169 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])));
170 for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 0.01 : 0;
171 accelerations(pl, f, lo, true, sr);
172 let Flo = 0;
173 for (let k = 0; k < g.n; k++) Flo += masses[k] * lo.w[k];
174 // Bottom pad unloaded: only the top pad (compressed by preload + 10 mm) pushes down, minus the rest balance.
175 const kb = sr.totalKw / 2;
176 const expected = -(kb * (sr.preloadTop + 0.01) - kb * sr.preloadTop) - kb * sr.preloadBottom;
177 assert.ok(Math.abs(Flo / expected - 1) < 1e-6, `lift-off force ${Flo} vs ${expected}`);
178});