d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 1// Float64 reference implementation of the spatial discretization. It mirrors
2// step.wgsl and is used for the time-step estimate, the tests, and the
3// GPU-versus-CPU self-test. It is not used in the time loop.
5import { type Grid, type Material, insideCell, insideNode, nodeWeight } from "../params.ts";
7export interface Fields {
8 u: Float64Array;
9 v: Float64Array;
10 w: Float64Array;
11}
13export function fields(n: number): Fields {
14 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
15}
17const G3 = 1 / Math.sqrt(3);
19/** Shape function derivatives of the bilinear cell at Gauss point gp (0..3), scaled by 2/h. */
20function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float64Array, dNdy: Float64Array): void {
21 const xi = gp & 1 ? G3 : -G3;
22 const eta = gp & 2 ? G3 : -G3;
23 dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
24 dNdx[1] = ((1 - eta) / 4) * (2 / hx);
25 dNdx[2] = ((1 + eta) / 4) * (2 / hx);
26 dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
27 dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
28 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
29 dNdy[2] = ((1 + xi) / 4) * (2 / hy);
30 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
31}
33/**
34 * Undamped accelerations (out.u, out.v, out.w) of every node for the
35 * displacement fields f and the stress-free shape w0. Bending forces are the
36 * gradient of the discrete Kirchhoff energy (natural free-edge conditions);
37 * membrane forces come from bilinear cells with 2x2 Gauss integration and
38 * the Marguerre form of the von Karman strains. Nodes outside the sheet get zero.
39 */
40export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?: Float64Array): void {
41 const { nx, ny, hx, hy } = g;
42 const { u, v, w } = f;
43 const hx2 = hx * hx, hy2 = hy * hy;
44 const nu = m.nu;
45 const id = (i: number, j: number) => j * nx + i;
47 // A second difference exists where the node and both of its neighbors belong to the sheet.
48 const hasA = (i: number, j: number) => insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j);
49 const hasB = (i: number, j: number) => insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1);
50 const dxx = (i: number, j: number) => (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
51 const dyy = (i: number, j: number) => (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
52 // Weighted bending moments at nodes; zero where the second difference is undefined.
53 const Ma = (i: number, j: number): number => {
54 if (!hasA(i, j)) return 0;
55 return nodeWeight(g, i, j) * (dxx(i, j) + nu * (hasB(i, j) ? dyy(i, j) : 0));
56 };
57 const Mb = (i: number, j: number): number => {
58 if (!hasB(i, j)) return 0;
59 return nodeWeight(g, i, j) * (dyy(i, j) + nu * (hasA(i, j) ? dxx(i, j) : 0));
60 };
61 // Mixed second difference at a cell center; zero for cells outside the sheet.
62 const Cc = (ci: number, cj: number): number => {
63 if (!insideCell(g, ci, cj)) return 0;
64 return (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
65 };
67 const q = (hx * hy) / 4; // area per Gauss point
68 const U = new Float64Array(4), V = new Float64Array(4), Wn = new Float64Array(4), W0 = new Float64Array(4);
69 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
71 // Adds the membrane force of cell (ci, cj) on its local node a to F = [Fu, Fv, Fw].
72 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
73 if (!insideCell(g, ci, cj)) return;
74 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
75 for (let b = 0; b < 4; b++) {
76 U[b] = u[nn[b]]; V[b] = v[nn[b]]; Wn[b] = w[nn[b]]; W0[b] = w0 ? w0[nn[b]] : 0;
77 }
78 for (let gp = 0; gp < 4; gp++) {
79 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
80 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
81 for (let b = 0; b < 4; b++) {
82 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
83 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
84 wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);
85 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
86 }
87 // Marguerre strains: total slope (w0 + w) minus the stress-free slope of w0.
88 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
89 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
90 const gxy = uy + vx + wx * wy - w0x * w0y;
91 const Nxx = m.C * (exx + nu * eyy);
92 const Nyy = m.C * (eyy + nu * exx);
93 const Nxy = m.G * gxy;
94 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
95 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
96 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
97 }
98 };
100 const F = new Float64Array(3);
101 for (let j = 0; j < ny; j++) {
102 for (let i = 0; i < nx; i++) {
103 const k = id(i, j);
104 const weight = nodeWeight(g, i, j);
105 if (weight === 0) {
106 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
107 continue;
108 }
109 const bend =
110 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
111 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
112 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
113 F[0] = 0; F[1] = 0; F[2] = -m.D * hx * hy * bend;
114 cellForce(i, j, 0, F);
115 cellForce(i - 1, j, 1, F);
116 cellForce(i - 1, j - 1, 2, F);
117 cellForce(i, j - 1, 3, F);
118 const mass = m.rhoh * hx * hy * weight;
119 out.u[k] = F[0] / mass;
120 out.v[k] = F[1] / mass;
121 out.w[k] = F[2] / mass;
122 }
123 }
124}
126/**
127 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
128 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
129 * transverse component, and -2 sigmaIn v on the in-plane components.
130 */
131export function addDamping(g: Grid, sigma0: number, sigma1: number, sigmaIn: number, vel: Fields, out: Fields): void {
132 const { nx, ny, hx, hy } = g;
133 const id = (i: number, j: number) => j * nx + i;
134 for (let j = 0; j < ny; j++) {
135 for (let i = 0; i < nx; i++) {
136 const k = id(i, j);
137 if (nodeWeight(g, i, j) === 0) continue;
138 const c = vel.w[k];
139 let lap = 0;
140 if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
141 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
142 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
143 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
144 out.u[k] -= 2 * (sigma0 + sigmaIn) * vel.u[k];
145 out.v[k] -= 2 * (sigma0 + sigmaIn) * vel.v[k];
146 out.w[k] += -2 * sigma0 * c + 2 * sigma1 * lap;
147 }
148 }
149}
151/** Discrete potential energy (bending + membrane) whose negative gradient is the force in `accelerations`. */
152export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Array): number {
153 const { nx, ny, hx, hy } = g;
154 const { u, v, w } = f;
155 const nu = m.nu;
156 const id = (i: number, j: number) => j * nx + i;
157 let Ub = 0;
158 for (let j = 0; j < ny; j++) {
159 for (let i = 0; i < nx; i++) {
160 let a = 0, b = 0;
161 if (insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
162 if (insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
163 Ub += nodeWeight(g, i, j) * (a * a + b * b + 2 * nu * a * b);
164 }
165 }
166 for (let cj = 0; cj < ny - 1; cj++) {
167 for (let ci = 0; ci < nx - 1; ci++) {
168 if (!insideCell(g, ci, cj)) continue;
169 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
170 Ub += 2 * (1 - nu) * c * c;
171 }
172 }
173 Ub *= 0.5 * m.D * hx * hy;
175 let Um = 0;
176 const q = (hx * hy) / 4;
177 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
178 for (let cj = 0; cj < ny - 1; cj++) {
179 for (let ci = 0; ci < nx - 1; ci++) {
180 if (!insideCell(g, ci, cj)) continue;
181 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
182 for (let gp = 0; gp < 4; gp++) {
183 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
184 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
185 for (let b = 0; b < 4; b++) {
186 const z0 = w0 ? w0[nn[b]] : 0;
187 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
188 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
189 wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);
190 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
191 }
192 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
193 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
194 const gxy = uy + vx + wx * wy - w0x * w0y;
195 Um += q * 0.5 * (m.C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + m.G * gxy * gxy);
196 }
197 }
198 }
199 return Ub + Um;
200}
202/**
203 * Largest angular frequency of the linearized discrete system, by power
204 * iteration on a small full grid with the same spacing. The extreme mode is a
205 * grid-scale in-plane mode, so neither the sheet size nor its outline matters.
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 206 * `inPlaneMass` is the factor by which the in-plane inertia is scaled.
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 208export function maxAngularFrequency(g: Grid, m: Material, inPlaneMass = 1, nSmall = 12, iters = 300): number {
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 209 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
210 const n = gs.n;
211 const x = fields(n), y = fields(n);
212 let seed = 12345;
213 const rnd = () => {
214 seed = (seed * 1664525 + 1013904223) >>> 0;
215 return seed / 4294967296 - 0.5;
216 };
217 for (let k = 0; k < n; k++) {
218 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
219 }
220 const eps = 1e-9; // keeps the nonlinear terms negligible
221 let lambda = 0;
222 for (let it = 0; it < iters; it++) {
223 let norm = 0;
224 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
225 norm = Math.sqrt(norm);
226 for (let k = 0; k < n; k++) {
227 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
228 }
229 accelerations(gs, m, x, y);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 230 if (inPlaneMass !== 1) {
231 for (let k = 0; k < n; k++) { y.u[k] /= inPlaneMass; y.v[k] /= inPlaneMass; }
232 }
234 for (let k = 0; k < n; k++) dot += x.u[k] * y.u[k] + x.v[k] * y.v[k] + x.w[k] * y.w[k];
235 lambda = -dot / (eps * eps); // Rayleigh quotient of -M^{-1} K
236 for (let k = 0; k < n; k++) {
237 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
238 }
239 }
240 return Math.sqrt(Math.max(lambda, 0));
241}