Initial-shape imperfection (Marguerre shallow-shell strains) - proof-of-concept/steelsim21

proof-of-concept / steelsim21

Initial-shape imperfection (Marguerre shallow-shell strains)

Real sheet metal has a bow of a few mm; a bowed sheet is a shallow shell with an asymmetric stiffness, which is the likely source of the slow wobble of a real sheet. Adds a stress-free initial shape (across width, across height, or dish) with configurable sag, applied through the Marguerre form of the membrane strains, a test of the energy gradient with the imperfection, and a 3D view option to draw the total shape.

Jeremy Magland <jmagland@flatironinstitute.org> committed 1 month ago commit 8dcc11fc223c parent b9a5a7e Browse files

10 changed files+177−29

README.mdmodified+2−1View file

@@ -15,7 +15,8 @@ Requires Node 20 or later and a browser with WebGPU (Chrome or Edge 113+; Firefo

1515

1616 Press Start. The sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles, and then released by ramping the finger force to zero over the release time; the recording of audio and movie frames starts at the release. The view shows the transverse displacement as a color map (2D) or a lit height field (3D, drag to orbit and wheel to zoom), the spectrogram shows the listener signal, and the audio can be played or downloaded as a 16-bit WAV at any time.

1717

18-- Plate, excitation, and numerics parameters take effect at the next Start. Damping, listener position, and the radiation model can be changed while a simulation runs.

18+- The sheet has a stress-free initial bow (default 3 mm across the width; set the shape to Flat for the ideal plate). A negative push amount pushes against the bow.

19+- Plate, initial shape, excitation, and numerics parameters take effect at the next Start. Damping, listener position, and the radiation model can be changed while a simulation runs.

1920 - After the simulation completes, drag the time slider to explore the movie; Play starts audio from the slider position and moves the movie and the spectrogram cursor with it. Live returns to the running view.

2021 - To extend a finished simulation, raise the duration and press Continue; the state is kept on the GPU.

2122 - Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1`.

THEORY.mdmodified+13−1View file

@@ -27,6 +27,18 @@ with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. The quadratic terms in the

2727

2828 The loss terms follow the two-parameter model that is standard in plate sound synthesis. A mode of angular frequency $\omega$ decays at the rate $\sigma_0 + \sigma_1 k^2$ with $k^2 = \omega/\kappa$ and $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster. The defaults ($\sigma_0 = 1.15\ \mathrm{s^{-1}}$, $\sigma_1 = 1.08\times10^{-4}\ \mathrm{m^2/s}$) give a 60 dB decay time of about 6 s at low frequencies and about 1 s near 10 kHz. These are not measured values; both parameters are exposed in the interface.

2929

30+### Initial shape (imperfection)

31+

32+Real sheet metal is never flat; rolling and handling leave a bow of a few millimetres. A bowed sheet is a shallow shell, and its behaviour differs qualitatively from a flat plate: the flat plate has a symmetric, purely hardening nonlinearity, whereas a shell has an asymmetric one, may snap between two curved states, and oscillates more slowly at large amplitude as it passes through the soft flat configuration. This is the mechanism behind oil-canning and the wobble board, and it appears to be what produces the slow wobble of a real sheet. We therefore allow a stress-free initial shape $w_0(x,y)$ in the Marguerre form of the von Karman equations: the membrane strains use the total slope minus the stress-free slope,

33+

34+$$

35+\varepsilon_{xx} = u_x + \tfrac12\left[(w_0 + w)_x^2 - w_{0,x}^2\right], \qquad

36+\varepsilon_{yy} = v_y + \tfrac12\left[(w_0 + w)_y^2 - w_{0,y}^2\right], \qquad

37+\gamma_{xy} = u_y + v_x + (w_0 + w)_x (w_0 + w)_y - w_{0,x} w_{0,y},

38+$$

39+

40+the membrane forces act through the total slope $(w_0 + w)$, and bending acts on the deviation $w$ alone. The shape $w_0 = 0$, $w = 0$ carries no stress. Three shapes are offered, all parabolic with a prescribed sag of the centre relative to the edges: curved across the width, curved across the height, and a dish curved both ways. The default is a 3 mm bow across the width. A positive push moves the corner in the direction of the bow and a negative push against it; the two cases sound different, as they do on a real sheet.

41+

3042 ## 2. Boundary conditions and suspension

3143

3244 All four edges are free: the bending moment and the Kirchhoff shear vanish, and so does the in-plane traction. Rather than imposing these conditions explicitly, the discretization of Section 4 is derived from a discrete energy, and free edges are then the natural boundary conditions of that energy. This avoids the intricate ghost-point formulas that free edges otherwise require, and it guarantees a symmetric stiffness operator.

@@ -73,4 +85,4 @@ The pressure is reported in pascals, and the peak sound pressure level is displa

7385

7486 ## 7. Limitations

7587

76-The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; the corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The loss model is empirical and its parameters are not calibrated against a real sheet. The free-edge treatment near the edges is only first-order accurate. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The radiation model is approximate, as discussed above. Finally, single precision limits the dynamic range of the tail of the sound to roughly 100 dB below the peak.

88+The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; the corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The Marguerre treatment of the initial shape is likewise a shallow-shell approximation, adequate for a bow of a few millimetres, and the real shape of a given sheet is of course unknown; the three offered shapes are idealizations. The loss model is empirical and its parameters are not calibrated against a real sheet. The free-edge treatment near the edges is only first-order accurate. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The radiation model is approximate, as discussed above. Finally, single precision limits the dynamic range of the tail of the sound to roughly 100 dB below the peak.

src/main.tsmodified+39−4View file

@@ -1,6 +1,6 @@

11 // UI wiring: parameters, run loop, live and movie views, playback, downloads.

22

3-import { defaultParams, makeGrid, makeMaterial, t60, type PhysParams } from "./params.ts";

3+import { defaultParams, makeGrid, makeMaterial, t60, type BowShape, type PhysParams } from "./params.ts";

44 import { Simulation } from "./sim/simulation.ts";

55 import { Renderer, type ViewMode } from "./render/renderer.ts";

66 import { Spectrogram } from "./audio/spectrogram.ts";

@@ -47,6 +47,11 @@ const groups: { title: string; open: boolean; fields: FieldSpec[] }[] = [

4747 { key: "holeInset", label: "Hole inset", unit: "mm", scale: 1e-3, step: 1 },

4848 ],

4949 },

50+ {

51+ title: "Initial shape",

52+ open: true,

53+ fields: [{ key: "bowAmplitude", label: "Bow (center vs edges)", unit: "mm", scale: 1e-3, step: 0.5 }],

54+ },

5055 {

5156 title: "Damping",

5257 open: true,

@@ -123,6 +128,27 @@ function buildForm(): void {

123128 det.appendChild(r2);

124129 inputs.set("radiation", sel);

125130 }

131+ if (f.key === "bowAmplitude") {

132+ const r2 = document.createElement("div");

133+ r2.className = "field";

134+ const l2 = document.createElement("label");

135+ l2.textContent = "Bow shape";

136+ const sel = document.createElement("select");

137+ for (const [v, t] of [["width", "Curved across the width"], ["height", "Curved across the height"], ["dish", "Dish (doubly curved)"], ["none", "Flat"]]) {

138+ const o = document.createElement("option");

139+ o.value = v; o.textContent = t;

140+ sel.appendChild(o);

141+ }

142+ sel.value = params.bowShape;

143+ sel.addEventListener("change", () => { params.bowShape = sel.value as BowShape; });

144+ r2.append(l2, sel, document.createElement("span"));

145+ det.appendChild(r2);

146+ inputs.set("bowShape", sel);

147+ const note = document.createElement("p");

148+ note.className = "sub";

149+ note.textContent = "Real sheet metal is never flat. A bow of a few mm makes the sheet a shallow shell with an asymmetric stiffness; push with a negative amount to push against the bow.";

150+ det.appendChild(note);

151+ }

126152 if (f.key === "sigma1") {

127153 const note = document.createElement("p");

128154 note.className = "sub";

@@ -165,6 +191,7 @@ function updateDampingNote(): void {

165191

166192 function setFormEnabled(enabled: boolean): void {

167193 for (const g of groups) for (const f of g.fields) if (!f.live) inputs.get(f.key)!.disabled = !enabled;

194+ inputs.get("bowShape")!.disabled = !enabled;

168195 }

169196

170197 // ---------- State ----------

@@ -255,6 +282,7 @@ const btnWav = mkButton("Download .wav");

255282 const levelSel = document.createElement("select");

256283 let exagInput: HTMLInputElement;

257284 let rangeInput: HTMLInputElement;

285+let shapeCheck: HTMLInputElement;

258286

259287 function buildViewToolbar(): void {

260288 const tb = $("viewToolbar");

@@ -270,7 +298,9 @@ function buildViewToolbar(): void {

270298 hint.className = "sub";

271299 hint.style.margin = "0";

272300 hint.textContent = "3D: drag to orbit, wheel to zoom";

273- tb.append(seg, labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), hint);

301+ shapeCheck = document.createElement("input");

302+ shapeCheck.type = "checkbox";

303+ tb.append(seg, labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), labeled("3D: add initial shape", shapeCheck), hint);

274304 }

275305

276306 function buildMovieToolbar(): void {

@@ -371,7 +401,7 @@ async function start(): Promise<void> {

371401 btnPause.textContent = "Pause";

372402 params.duration = duration;

373403 sim = new Simulation(device, params);

374- renderer = new Renderer(device, $("view"), sim.grid, params);

404+ renderer = new Renderer(device, $("view"), sim.grid, params, sim.w0);

375405 movieBuf?.destroy();

376406 movieBuf = device.createBuffer({ size: Math.max(16, sim.grid.n * 4), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });

377407 movieBufFrame = -1;

@@ -427,6 +457,7 @@ function frameBody(): void {

427457 renderer.mode = viewMode.value;

428458 renderer.exaggeration = parseFloat(exagInput.value) || 1;

429459 renderer.colorRange = (parseFloat(rangeInput.value) || 0) * 1e-3;

460+ renderer.showShape = shapeCheck.checked;

430461 const pos = player.position();

431462 if (pos !== null) {

432463 scrubFrame = Math.min(sim.frames.length - 1, Math.floor(pos * sim.fps));

@@ -509,7 +540,7 @@ async function init(): Promise<void> {

509540 function applyUrlParams(): void {

510541 const q = new URLSearchParams(location.search);

511542 for (const [k, v] of q) {

512- if (k in params && k !== "radiation") {

543+ if (k in params && k !== "radiation" && k !== "bowShape") {

513544 const num = parseFloat(v);

514545 if (isFinite(num)) {

515546 (params as unknown as Record<string, number>)[k] = num;

@@ -519,6 +550,10 @@ function applyUrlParams(): void {

519550 }

520551 }

521552 if (k === "radiation" && (v === "rayleigh" || v === "dipole")) params.radiation = v;

553+ if (k === "bowShape" && ["none", "width", "height", "dish"].includes(v)) {

554+ params.bowShape = v as BowShape;

555+ (inputs.get("bowShape") as HTMLSelectElement).value = v;

556+ }

522557 if (k === "duration") {

523558 duration = params.duration;

524559 durInput.value = String(duration);

src/params.tsmodified+24−0View file

@@ -25,8 +25,12 @@ export interface PhysParams {

2525 listenerAzimuth: number; // [deg], rotation about the plate's vertical axis (0 = on the normal)

2626 listenerElevation: number; // [deg]

2727 radiation: "rayleigh" | "dipole";

28+ bowShape: BowShape; // stress-free initial shape of the sheet

29+ bowAmplitude: number; // sag of the initial shape, center relative to the edges [m]

2830 }

2931

32+export type BowShape = "none" | "width" | "height" | "dish";

33+

3034 export const defaultParams: PhysParams = {

3135 Lx: 18 * 0.0254,

3236 Ly: 12 * 0.0254,

@@ -52,8 +56,28 @@ export const defaultParams: PhysParams = {

5256 listenerAzimuth: 0,

5357 listenerElevation: 0,

5458 radiation: "rayleigh",

59+ bowShape: "width",

60+ bowAmplitude: 3e-3,

5561 };

5662

63+/** Stress-free initial shape w0 at every node. */

64+export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {

65+ const w0 = new Float32Array(g.n);

66+ if (p.bowShape === "none" || p.bowAmplitude === 0) return w0;

67+ for (let j = 0; j < g.ny; j++) {

68+ for (let i = 0; i < g.nx; i++) {

69+ const X = (2 * nodeX(g, i)) / p.Lx; // -1..1

70+ const Y = (2 * nodeY(g, j)) / p.Ly;

71+ let f = 0;

72+ if (p.bowShape === "width") f = 1 - X * X;

73+ else if (p.bowShape === "height") f = 1 - Y * Y;

74+ else f = 1 - 0.5 * (X * X + Y * Y);

75+ w0[j * g.nx + i] = p.bowAmplitude * f;

76+ }

77+ }

78+ return w0;

79+}

80+

5781 export const SPEED_OF_SOUND = 343; // [m/s]

5882 export const RHO_AIR = 1.2; // [kg/m^3]

5983 export const HIST_LEN = 256; // ring buffer length for the retarded-time audio sum [samples]

src/render/renderer.tsmodified+11−3View file

@@ -12,6 +12,8 @@ export class Renderer {

1212 exaggeration = 4;

1313 /** Color range |w| [m]; 0 means automatic (from the current maximum). */

1414 colorRange = 0;

15+ /** Add the stress-free initial shape to the 3D height field. */

16+ showShape = false;

1517 camera = { yaw: 0.6, pitch: 0.5, dist: 1.0 };

1618

1719 private device: GPUDevice;

@@ -29,20 +31,24 @@ export class Renderer {

2931 private params: PhysParams;

3032 private bindGroups = new Map<GPUBuffer, GPUBindGroup>();

3133

32- constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams) {

34+ private w0: GPUBuffer;

35+

36+ constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams, w0: GPUBuffer) {

3337 this.device = device;

38+ this.w0 = w0;

3439 this.canvas = canvas;

3540 this.grid = grid;

3641 this.params = params;

3742 this.ctx = canvas.getContext("webgpu") as GPUCanvasContext;

3843 this.format = navigator.gpu.getPreferredCanvasFormat();

3944 this.ctx.configure({ device, format: this.format, alphaMode: "opaque" });

40- this.uniform = device.createBuffer({ size: 128, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });

45+ this.uniform = device.createBuffer({ size: 144, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });

4146 const module = device.createShaderModule({ code: shaderCode });

4247 this.bgLayout = device.createBindGroupLayout({

4348 entries: [

4449 { binding: 0, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "uniform" } },

4550 { binding: 1, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },

51+ { binding: 2, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },

4652 ],

4753 });

4854 const layout = device.createPipelineLayout({ bindGroupLayouts: [this.bgLayout] });

@@ -114,7 +120,7 @@ export class Renderer {

114120 this.resize();

115121 const aspect = this.canvas.width / this.canvas.height;

116122 const range = this.colorRange > 0 ? this.colorRange : Math.max(wmax, 1e-6);

117- const u = new Float32Array(32);

123+ const u = new Float32Array(36);

118124 const ui = new Uint32Array(u.buffer);

119125 const eye = [

120126 this.camera.dist * Math.cos(this.camera.pitch) * Math.sin(this.camera.yaw),

@@ -132,6 +138,7 @@ export class Renderer {

132138 u[26] = p.Lx / 2; u[27] = -p.Ly / 2;

133139 u[28] = p.fingerRadius;

134140 u[29] = 0.4; u[30] = 0.5; u[31] = 1.0;

141+ ui[32] = this.showShape ? 1 : 0;

135142 this.device.queue.writeBuffer(this.uniform, 0, u);

136143

137144 let bg = this.bindGroups.get(wBuffer);

@@ -141,6 +148,7 @@ export class Renderer {

141148 entries: [

142149 { binding: 0, resource: { buffer: this.uniform } },

143150 { binding: 1, resource: { buffer: wBuffer } },

151+ { binding: 2, resource: { buffer: this.w0 } },

144152 ],

145153 });

146154 this.bindGroups.set(wBuffer, bg);

src/render/shaders.wgslmodified+16−3View file

@@ -16,10 +16,15 @@ struct RU {

1616 lightX: f32,

1717 lightY: f32,

1818 lightZ: f32,

19+ showShape: u32, // 3D: add the stress-free initial shape to the height field

20+ pad0: u32,

21+ pad1: u32,

22+ pad2: u32,

1923 };

2024

2125 @group(0) @binding(0) var<uniform> R: RU;

2226 @group(0) @binding(1) var<storage, read> w: array<f32>;

27+@group(0) @binding(2) var<storage, read> w0: array<f32>;

2328

2429 fn colormap(t: f32) -> vec3f {

2530 let s = clamp(t, -1.0, 1.0);

@@ -35,6 +40,14 @@ fn wAt(i: i32, j: i32) -> f32 {

3540 let jj = clamp(j, 0, i32(R.ny) - 1);

3641 return w[u32(jj) * R.nx + u32(ii)];

3742 }

43+// Height used by the 3D view: w, plus the initial shape when requested.

44+fn zAt(i: i32, j: i32) -> f32 {

45+ let ii = clamp(i, 0, i32(R.nx) - 1);

46+ let jj = clamp(j, 0, i32(R.ny) - 1);

47+ var z = w[u32(jj) * R.nx + u32(ii)];

48+ if (R.showShape != 0u) { z += w0[u32(jj) * R.nx + u32(ii)]; }

49+ return z;

50+}

3851

3952 // ---------- 2D view: full-screen triangle, bilinear lookup in the fragment shader ----------

4053

@@ -87,9 +100,9 @@ fn vs3d(@builtin(vertex_index) vi: u32) -> V3 {

87100 let x = (f32(i) - 0.5 * f32(R.nx - 1u)) * R.hx;

88101 let y = (f32(j) - 0.5 * f32(R.ny - 1u)) * R.hy;

89102 let val = wAt(i, j);

90- let z = val * R.wScale;

91- let dzdx = (wAt(i + 1, j) - wAt(i - 1, j)) * R.wScale / (2.0 * R.hx);

92- let dzdy = (wAt(i, j + 1) - wAt(i, j - 1)) * R.wScale / (2.0 * R.hy);

103+ let z = zAt(i, j) * R.wScale;

104+ let dzdx = (zAt(i + 1, j) - zAt(i - 1, j)) * R.wScale / (2.0 * R.hx);

105+ let dzdy = (zAt(i, j + 1) - zAt(i, j - 1)) * R.wScale / (2.0 * R.hy);

93106 var o: V3;

94107 o.pos = R.mvp * vec4f(x, y, z, 1.0);

95108 o.val = val;

src/sim/cpuRef.tsmodified+18−12View file

@@ -18,7 +18,7 @@ const G3 = 1 / Math.sqrt(3);

1818 * (natural free-edge conditions); membrane forces come from bilinear (Q4)

1919 * elements with 2x2 Gauss integration and the von Karman strains.

2020 */

21-export function accelerations(g: Grid, m: Material, f: Fields, out: Fields): void {

21+export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?: Float64Array): void {

2222 const { nx, ny, hx, hy } = g;

2323 const { u, v, w } = f;

2424 const hx2 = hx * hx;

@@ -51,6 +51,7 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields): voi

5151 const U = new Float64Array(4);

5252 const V = new Float64Array(4);

5353 const Wn = new Float64Array(4);

54+ const W0 = new Float64Array(4);

5455 const dNdx = new Float64Array(4);

5556 const dNdy = new Float64Array(4);

5657

@@ -61,6 +62,7 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields): voi

6162 U[0] = u[n0]; U[1] = u[n1]; U[2] = u[n2]; U[3] = u[n3];

6263 V[0] = v[n0]; V[1] = v[n1]; V[2] = v[n2]; V[3] = v[n3];

6364 Wn[0] = w[n0]; Wn[1] = w[n1]; Wn[2] = w[n2]; Wn[3] = w[n3];

65+ if (w0) { W0[0] = w0[n0]; W0[1] = w0[n1]; W0[2] = w0[n2]; W0[3] = w0[n3]; }

6466 for (let gp = 0; gp < 4; gp++) {

6567 const xi = gp & 1 ? G3 : -G3;

6668 const eta = gp & 2 ? G3 : -G3;

@@ -72,15 +74,17 @@ export function accelerations(g: Grid, m: Material, f: Fields, out: Fields): voi

7274 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);

7375 dNdy[2] = ((1 + xi) / 4) * (2 / hy);

7476 dNdy[3] = ((1 - xi) / 4) * (2 / hy);

75- let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0;

77+ // Marguerre strains: total slope (w0 + w) minus the stress-free slope of w0.

78+ let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;

7679 for (let b = 0; b < 4; b++) {

7780 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];

7881 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];

79- wx += dNdx[b] * Wn[b]; wy += dNdy[b] * Wn[b];

82+ wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);

83+ w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];

8084 }

81- const exx = ux + 0.5 * wx * wx;

82- const eyy = vy + 0.5 * wy * wy;

83- const gxy = uy + vx + wx * wy;

85+ const exx = ux + 0.5 * (wx * wx - w0x * w0x);

86+ const eyy = vy + 0.5 * (wy * wy - w0y * w0y);

87+ const gxy = uy + vx + wx * wy - w0x * w0y;

8488 const Nxx = m.C * (exx + nu * eyy);

8589 const Nyy = m.C * (eyy + nu * exx);

8690 const Nxy = m.G * gxy;

@@ -150,7 +154,7 @@ export function maxAngularFrequency(g: Grid, m: Material, nSmall = 12, iters = 3

150154 }

151155

152156 /** Discrete potential energy (bending + membrane) matching `accelerations`. */

153-export function potentialEnergy(g: Grid, m: Material, f: Fields): number {

157+export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Array): number {

154158 const { nx, ny, hx, hy } = g;

155159 const { u, v, w } = f;

156160 const nu = m.nu;

@@ -182,15 +186,17 @@ export function potentialEnergy(g: Grid, m: Material, f: Fields): number {

182186 const eta = gp & 2 ? G3 : -G3;

183187 const dNdx = [-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)].map((s) => (s / 4) * (2 / hx));

184188 const dNdy = [-(1 - xi), -(1 + xi), 1 + xi, 1 - xi].map((s) => (s / 4) * (2 / hy));

185- let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0;

189+ let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;

186190 for (let b = 0; b < 4; b++) {

191+ const z0 = w0 ? w0[nn[b]] : 0;

187192 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];

188193 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];

189- wx += dNdx[b] * w[nn[b]]; wy += dNdy[b] * w[nn[b]];

194+ wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);

195+ w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;

190196 }

191- const exx = ux + 0.5 * wx * wx;

192- const eyy = vy + 0.5 * wy * wy;

193- const gxy = uy + vx + wx * wy;

197+ const exx = ux + 0.5 * (wx * wx - w0x * w0x);

198+ const eyy = vy + 0.5 * (wy * wy - w0y * w0y);

199+ const gxy = uy + vx + wx * wy - w0x * w0y;

194200 Um += q * 0.5 * (m.C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + m.G * gxy * gxy);

195201 }

196202 }

src/sim/shaders/step.wgslmodified+9−4View file

@@ -14,6 +14,7 @@

1414 @group(0) @binding(5) var<storage, read_write> velOut: array<f32>;

1515 @group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node] acceleration of w

1616 @group(0) @binding(7) var<storage, read_write> react: array<f32>; // constraint reaction (acceleration of w) per node

17+@group(0) @binding(8) var<storage, read> w0: array<f32>; // stress-free initial shape

1718

1819 const TILE = 16;

1920 const HALO = 2;

@@ -37,6 +38,7 @@ fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }

3738 fn U(i: i32, j: i32) -> f32 { return posIn[id(i, j)]; }

3839 fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }

3940 fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }

41+fn W0(i: i32, j: i32) -> f32 { return w0[id(i, j)]; }

4042

4143 // Weighted bending "moments" at nodes; zero where the second difference is undefined.

4244 fn Ma(i: i32, j: i32) -> f32 {

@@ -71,7 +73,8 @@ fn cellForces(ci: i32, cj: i32, lc: u32) {

7173 }

7274 let Un = vec4f(U(ci, cj), U(ci + 1, cj), U(ci + 1, cj + 1), U(ci, cj + 1));

7375 let Vn = vec4f(V(ci, cj), V(ci + 1, cj), V(ci + 1, cj + 1), V(ci, cj + 1));

74- let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1));

76+ let W0n = vec4f(W0(ci, cj), W0(ci + 1, cj), W0(ci + 1, cj + 1), W0(ci, cj + 1));

77+ let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1)) + W0n; // total shape

7578 let q = 0.25 * P.hx * P.hy;

7679 let ihx = 1.0 / P.hx; let ihy = 1.0 / P.hy;

7780 let g3 = 0.57735026918962576;

@@ -84,9 +87,11 @@ fn cellForces(ci: i32, cj: i32, lc: u32) {

8487 let ux = dot(dNdx, Un); let uy = dot(dNdy, Un);

8588 let vx = dot(dNdx, Vn); let vy = dot(dNdy, Vn);

8689 let wx = dot(dNdx, Wn); let wy = dot(dNdy, Wn);

87- let exx = ux + 0.5 * wx * wx;

88- let eyy = vy + 0.5 * wy * wy;

89- let gxy = uy + vx + wx * wy;

90+ let w0x = dot(dNdx, W0n); let w0y = dot(dNdy, W0n);

91+ // Marguerre strains: total slope minus the stress-free slope.

92+ let exx = ux + 0.5 * (wx * wx - w0x * w0x);

93+ let eyy = vy + 0.5 * (wy * wy - w0y * w0y);

94+ let gxy = uy + vx + wx * wy - w0x * w0y;

9095 let Nxx = P.C * (exx + P.nu * eyy);

9196 let Nyy = P.C * (eyy + P.nu * exx);

9297 let Nxy = P.G * gxy;

src/sim/simulation.tsmodified+8−1View file

@@ -12,6 +12,7 @@ import {

1212 type Material,

1313 makeGrid,

1414 makeMaterial,

15+ initialShape,

1516 nodeX,

1617 nodeY,

1718 trapWeight,

@@ -82,6 +83,8 @@ export class Simulation {

8283 private vel: GPUBuffer[];

8384 private hist: GPUBuffer;

8485 private react: GPUBuffer;

86+ /** Stress-free initial shape (f32 per node), also used by the renderer. */

87+ readonly w0: GPUBuffer;

8588 private lw: GPUBuffer;

8689 private audioOut: GPUBuffer;

8790 private frameBuf: GPUBuffer;

@@ -131,6 +134,8 @@ export class Simulation {

131134 this.vel = [f32(3 * n), f32(3 * n)]; // [du | dv | dw]

132135 this.hist = f32(HIST_LEN * n);

133136 this.react = f32(n);

137+ this.w0 = f32(n);

138+ device.queue.writeBuffer(this.w0, 0, initialShape(params, g));

134139 this.lw = f32(2 * n);

135140 this.audioOut = f32(this.samplesPerBatch);

136141 this.frameBuf = f32(n);

@@ -186,6 +191,7 @@ export class Simulation {

186191 entry(5, "storage"),

187192 entry(6, "storage"),

188193 entry(7, "storage"),

194+ entry(8, "read-only-storage"),

189195 ],

190196 });

191197 this.stepPipe = mk(stepWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [stepLayout] }));

@@ -230,6 +236,7 @@ export class Simulation {

230236 { binding: 5, resource: { buffer: this.vel[1 - p] } },

231237 { binding: 6, resource: { buffer: this.hist } },

232238 { binding: 7, resource: { buffer: this.react } },

239+ { binding: 8, resource: { buffer: this.w0 } },

233240 ],

234241 }),

235242 );

@@ -564,7 +571,7 @@ export class Simulation {

564571 }

565572

566573 destroy(): void {

567- for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();

574+ for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();

568575 }

569576

570577 get fs(): number {

tests/cpuRef.test.tsmodified+37−0View file

@@ -148,3 +148,40 @@ test("maximum frequency estimate is close to the in-plane CFL estimate", () => {

148148 const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)

149149 assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);

150150 });

151+

152+test("forces are the energy gradient with an initial-shape imperfection", () => {

153+ const nx = 8, ny = 7;

154+ const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };

155+ const f = fields(g.n), acc = fields(g.n);

156+ const w0 = new Float64Array(g.n);

157+ let seed = 11;

158+ const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);

159+ for (let k = 0; k < g.n; k++) {

160+ const x = (k % nx) / (nx - 1) - 0.5;

161+ w0[k] = 0.003 * (1 - 4 * x * x); // cylindrical bow, 3 mm

162+ f.w[k] = 0.01 * rnd();

163+ f.u[k] = 1e-4 * rnd();

164+ f.v[k] = 1e-4 * rnd();

165+ }

166+ accelerations(g, mat, f, acc, w0);

167+ const delta = 1e-7;

168+ let worst = 0;

169+ for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {

170+ for (let k = 0; k < g.n; k += 5) {

171+ const mass = mat.rhoh * g.hx * g.hy * trapWeight(g, k % nx, Math.floor(k / nx));

172+ const orig = arr[k];

173+ arr[k] = orig + delta;

174+ const Ep = potentialEnergy(g, mat, f, w0);

175+ arr[k] = orig - delta;

176+ const Em = potentialEnergy(g, mat, f, w0);

177+ arr[k] = orig;

178+ const Ffd = -(Ep - Em) / (2 * delta);

179+ const Fan = accArr[k] * mass;

180+ worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));

181+ }

182+ }

183+ assert.ok(worst < 1e-4, `gradient mismatch ${worst}`);

184+ // The stress-free shape itself carries no membrane energy.

185+ const rest = fields(g.n);

186+ assert.ok(potentialEnergy(g, mat, rest, w0) < 1e-12, "w0 is stress free");

187+});