Plate: clamped edges (clamped 13-point stencil, computed modes) and a general membrane tension T for held edges; mallet unchanged
5 changed files+100−45
src/index.htmlmodified+1−1View file
@@ -34,7 +34,7 @@
3434 <li><a href="./double-pendulum/"><strong>Double pendulum</strong><div class="muted">Strike it, watch and hear the motion turn regular as the energy decays; Lyapunov exponent and Poincaré section.</div></a></li>
3535 <li><a href="./fput-chain/"><strong>FPUT chain</strong><div class="muted">A struck chain of N masses with cubic springs; energy spreading up the modes, Lyapunov exponent and Poincaré sections for N = 2 and 3.</div></a></li>
3636 <li><a href="./beam/"><strong>Beam</strong><div class="muted">A struck beam of N masses with bending stiffness and stretching: clamped, pinned or cantilever; inharmonic, cymbal-like overtones.</div></a></li>
37- <li><a href="./plate/"><strong>Plate</strong><div class="muted">A struck square plate of N × N masses with bending stiffness and von Kármán stretching, simply supported or free; modes (m, n) at m² + n², pitch glides when struck hard.</div></a></li>
37+ <li><a href="./plate/"><strong>Plate</strong><div class="muted">A struck square plate of N × N masses with bending stiffness and von Kármán stretching, simply supported, clamped or free, optional membrane tension (drum-like); modes (m, n) at m² + n², pitch glides when struck hard.</div></a></li>
3838 </ul>
3939 </main>
4040 </body>
src/plate/README.mdmodified+7−1View file
@@ -1,6 +1,6 @@
11 # plate
22
3-Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram; step 2 (2026-10-09), von Kármán stretching; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices).
3+Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram; step 2 (2026-10-09), von Kármán stretching; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices); step 5 (2026-10-09), clamped edges and membrane tension.
44
55 Current state: N × N masses (N per side, 2 to 32, default 8) on vertical rails at the nodes of a square grid, spacing a, moving up and down only. Bending κ on the five-point Laplacian at each node, energy κ (L u)² / 2, force −κ L(L u). The edge nodes are supports at u = 0 with no bending term (mirror ghost), which makes the edge simply supported. Stretching (von Kármán), below, set by the thickness h (default 5 cm; 0 = linear). Damping as in the beam: drag (τ, all modes alike) and a frequency-dependent loss (dashpots alongside the grid springs, force −c L v, set by τ top at the top of the band). Parameters are given for the 8 × 8 plate (m 1 kg, κ 20 N/m, h 0.05 m, a 0.5 m, so the side is 4.5 m; τ 540 s, τ top 200 s); Masses per side refines the same plate (`refine`: m with the cell area, κ inversely, the loss coefficient c unchanged). Strike: a kick (pyramid velocity bump of half-width 0.5 m) or the shared Hertz mallet (`../common/mallet.ts`, M 1 kg, K 1000) over the same bump, at (x, y) m from the lower left corner (default (1.3, 1.1)); clicking the plate moves the strike point. The card shows the plate from above as a height map (red up, blue down, scaled to the recent largest |u|), with the strike point (cross) and the listening point (circle). The sound panel is the beam's, reduced: velocity, displacement or acceleration at a listening point (default (3.1, 1.7), bilinear between masses), pitch of the lowest mode, the whole decay as a note through the shared renderer and anti-alias filter, spectrogram with dashed lines at the distinct mode pitches, energy chart, Replay and Save WAV. Not yet: Lyapunov exponent, ghost, mode bars and histories, cross-section.
66
@@ -19,3 +19,9 @@ Check against plate theory (Leissa, free square plate, ν = 0.3, ω L² √(ρh/
1919 Stretching with free edges (step 4, 2026-10-09): the edge is free in its own plane too (no normal or shear stress), F = ∂F/∂n = 0, so ∇⁴F has clamped edge conditions. s = l(u, u) at the interior nodes; φ = B⁻¹ s with B the clamped biharmonic on them, F = 0 on the edge nodes and an even ghost beyond (the standard 13-point stencil: 20 inside, 21 next to an edge, 22 at a corner): B = L² + 2C, C counting the edges a node touches. B differs from L² only on the 4(N − 3) nodes next to the edge, so B⁻¹ comes from two sine solves and a small dense capacitance matrix (Woodbury, computed once per N); against a dense solve 1e-13. Energy (ks/8) Σ s φ, force its exact gradient (checked against central differences, 1e-10), invariant under rigid motion (no net force or torque, 1e-16); a cylindrical bend (developable) has exactly zero stretching energy. Independent check: a spherical cap u = c r² (c = 0.01 /m) on a free 4 m square, against a continuum Galerkin solution of the clamped biharmonic (Legendre-based polynomials, degree 27, Python): 0.0027848 J; grid +20.3%, +5.5%, +1.4%, +0.36% at N = 9, 17, 33, 65 (second order). With the cell-centred layout the same test gave −44%, −28%, −16%, −8.7% at N = 8, 16, 32, 64 (in-plane edge on the outer masses), or +70% … +10% with it one ring outside; that is why the layout was changed. Energy balance at 100 J (stretching up to 40% of the energy): 5e-4. Nonlinear scale (lowest mode, h 5 cm, N = 8): 0.140 J. Cost per step at N = 8: 1.7 µs linear, 4.4 µs with stretching, 6.0 µs with the loss too (simply supported: 4.8 µs); the default note renders in about 6.5 s in Node, 5.5 s in the browser (simply supported: 4.3 s).
2020
2121 Sound (2026-10-09, free, N = 8, h 5 cm, mallet at (1.3, 1.1), velocity at (3.1, 1.7), lowest mode at 100 Hz): at 0.01 J linear; at 1 J the energy already spreads between modes within the first second; at 10 J a broadband attack; at 100 J a wash for about the first second, the low modes emerging as it decays. Spectral centroid at 0.2 s: 880, 1210, 1270, 1290, 1500 Hz for 0.01, 1, 2.51, 10, 100 J.
22+
23+Clamped edges and tension (step 5, 2026-10-09). Edges select: simply supported, clamped, free. Clamped: the simply supported layout with the slope held at zero (even ghost beyond the supports); the support nodes then carry a bending term (L u = −2 u next to them, weight 1/2), so the energy is κ/2 (|L u|² + 2 Σ c u²), c the edges a node touches: the standard clamped 13-point stencil (21 next to an edge, 22 in a corner), the same operator as the free plate's stress function. Modes computed by the symmetry blocks (no (m, n)); stretching as simply supported (same in-plane edge). Check against Leissa's clamped square plate (ω L² √(ρh/D): 35.99, 73.39 (×2), 108.22, 131.58, 132.20, 165.15 (×2)): 33.62, 64.96, 94.51, 107.02 at N = 7; 35.33, 70.98, 104.10, 124.28 at N = 15; 35.82, 72.77, 107.12, 129.66 at N = 31 (lowest 6.6%, 1.8%, 0.47% low: second order). Time domain (N = 8, undamped kick): all 48 distinct frequencies within 3 cents (median 0.22), energy error 1.1e-3. Nonlinear scale (h 5 cm, N = 8): 0.44 J (simply supported 0.13 J).
24+
25+Tension T (N/m, Parameters; default 0; held edges only, disabled for free): a membrane tension, springs between neighbours and to the supports, energy T/2 uᵀLu, force −T L u, the same T at every N (T a²/a² per spring). Simply supported, the modes stay exact sines with ω² = (T l + κ l²)/m; the derived values compare the ratios with the continuum's, from ω² ∝ T k² + κ k⁴. With tension dominant (small κ) a held plate becomes a drumhead: ratios approach the square membrane's √(m² + n²) (clamped, κ ≈ 0: 1.580, 1.998, 2.230, 2.543 at N = 32 against 1.581, 2, 2.236, 2.550). Stretching then adds to the tension (tension modulation). A drum-like setting needs a thinner plate for its stretching to matter (the tension raises the linear stiffness): e.g. κ 2, T 200, h 1 mm gives a nonlinear scale near 1 J (simply supported, N = 8). Forces checked against central differences (1e-10, with tension, both held edges). A "drum" edge option was tried first and dropped for the general tension parameter.
26+
27+Known cost issue (2026-10-09): the 4x oversampled note forces at least 4 steps per audio sample, and the step stability needs is longer than that whenever the speed-up is small, that is when the lowest mode is close in frequency to the top of the band: clamped (lowest mode higher) and especially with tension (top/lowest only 7.5 at N = 8, so about 8x more steps than needed). Clamped default note: 12 s to render in the browser. Fix planned: read the signal between steps by interpolation (shared note renderer).
src/plate/index.htmlmodified+5−3View file
@@ -35,6 +35,7 @@
3535 <span class="muted" style="margin-left: 12px">Edges</span>
3636 <select id="edges">
3737 <option value="simply" selected title="held at the edge, free to turn there">simply supported</option>
38+ <option value="clamped" title="held at the edge and kept level there">clamped</option>
3839 <option value="free" title="nothing holds the plate (rigid motion removed)">free</option>
3940 </select>
4041 <span class="muted" style="margin-left: 12px">Speed</span>
@@ -55,6 +56,7 @@
5556 <div id="param-form">
5657 <label>m <input id="p-m" type="number" step="any" /> <span class="unit">kg</span></label>
5758 <label>κ <input id="p-kappa" type="number" step="any" /> <span class="unit">N/m (bending)</span></label>
59+ <label>T <input id="p-tension" type="number" step="any" /> <span class="unit">N/m (membrane tension, held edges; 0 = none)</span></label>
5860 <label>h <input id="p-thick" type="number" step="any" /> <span class="unit">m (thickness, sets the stretching; 0 = linear)</span></label>
5961 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
6062 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
@@ -102,15 +104,15 @@
102104
103105 <template id="help-plate">
104106 <h2>The plate</h2>
105- <p>A square plate as a grid of equal masses on vertical rails, moving up and down only, seen from above: red is up, blue down, white at rest (scaled to fit). Each mass resists bending (κ, on how far it sits off the mean of its four neighbours). The edges are simply supported (held at zero height but free to turn there) or free. A free plate would also slide and tilt as a whole; that motion is removed (the plate is followed in its own frame), and the readout says how much of a strike's energy it took. A bent plate also stretches in its own plane (von Kármán), and the in-plane stress stiffens it: this is the nonlinearity. It matters once the plate moves by about its thickness h, so weak strikes ring as independent normal modes, while strong ones start sharp and glide down in pitch as they decay, and the strongest spread their energy over many modes.</p>
107+ <p>A square plate as a grid of equal masses on vertical rails, moving up and down only, seen from above: red is up, blue down, white at rest (scaled to fit). Each mass resists bending (κ, on how far it sits off the mean of its four neighbours). The edges are simply supported (held at zero height but free to turn there), clamped (held and kept level), or free. A held plate can also carry a membrane tension T (Parameters), as a drumhead does. A free plate would also slide and tilt as a whole; that motion is removed (the plate is followed in its own frame), and the readout says how much of a strike's energy it took. A bent plate also stretches in its own plane (von Kármán), and the in-plane stress stiffens it: this is the nonlinearity. It matters once the plate moves by about its thickness h, so weak strikes ring as independent normal modes, while strong ones start sharp and glide down in pitch as they decay, and the strongest spread their energy over many modes.</p>
106108 <p>Strike gives the energy on the slider as a hammer blow (a short pyramid bump of velocity, one spacing of the 8 × 8 plate wide on each side) or with a mallet (mass M, contact stiffness K in Parameters) that presses into the plate over the same bump through a Hertz contact (force K δ^1.5) and lifts off. The point is in metres from the lower left corner; click the plate to move it there. The cross marks the strike point, the circle the listening point. The space bar pauses.</p>
107109 </template>
108110 <template id="help-params">
109111 <h2>Parameters</h2>
110- <p>Mass, bending κ, the thickness h (it sets the stretching, as the amplitude at which it matters; 0 for a linear plate), the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
112+ <p>Mass, bending κ, the membrane tension T (force per unit length, for simply supported or clamped edges; the same at every N), the thickness h (it sets the stretching, as the amplitude at which it matters; 0 for a linear plate), the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
111113 <p>The values are for the plate with 8 masses per side. Masses per side refines the same plate: same size, strike energies and decay, with each mass scaled with the area of its cell and κ inversely (the thickness stays), so the low modes stay put and more masses add higher ones.</p>
112114 <p>The stretching follows the von Kármán plate equations: the in-plane stress comes from a stress function, solved on the grid at every step. A simply supported edge may move outward but not along itself; a free edge is free in its own plane too. The nonlinear scale below is the energy of the lowest mode at which its stretching energy equals its bending energy.</p>
113- <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. For a simply supported plate, the normal modes are products of sines, numbered (m, n) by their half waves along x and y. Their frequencies go as m² + n² for the low ones (the continuum plate), so (1, 2) and (2, 1) share a pitch, and the ratios to the lowest, 2.5, 4, 5, 6.5, 8.5, ..., are not harmonic. The derived values list the lowest modes and how close the grid is to the continuum.</p>
115+ <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. Clamped edges hold the slope at zero too; the modes are computed (no (m, n)), and the plate is stiffer than a simply supported one, with overtones about as widely spread (plate theory: 1, 2.04, 3.01, 3.66, ...). With tension T and little bending (small κ), a held plate becomes a drumhead: its frequencies approach the ideal square membrane's, √(m² + n²), ratios 1, 1.58, 2, 2.24, 2.55, ..., also inharmonic but closer together than a plate's; stretching adds to the tension as it moves, so a hard strike starts sharp and settles, as on a tom-tom. For a simply supported plate, the normal modes are products of sines, numbered (m, n) by their half waves along x and y. Their frequencies go as m² + n² for the low ones (the continuum plate), so (1, 2) and (2, 1) share a pitch, and the ratios to the lowest, 2.5, 4, 5, 6.5, 8.5, ..., are not harmonic. The derived values list the lowest modes and how close the grid is to the continuum.</p>
114116 </template>
115117 <template id="help-sound">
116118 <h2>Sped up · listen</h2>
src/plate/main.tsmodified+20−11View file
@@ -36,10 +36,10 @@ const pauseButton = document.getElementById("pause") as HTMLButtonElement;
3636 // Parameters, as the beam's defaults: bending kappa = 20 on unit masses,
3737 // drag tau 540 s, the frequency-dependent loss tauTop 200 s; thickness 5 cm
3838 // for the stretching (2026-10-09: the default strike glides down by about
39-// half an octave, 100 J is a broadband wash). base: the page's
39+// half an octave, 100 J is a broadband wash); no tension. base: the page's
4040 // parameters, for the plate at REF_N masses per side; P: the same plate
4141 // refined to the chosen N (see refine).
42-let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
42+let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, tension: 0, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
4343 let P: Params = base;
4444 // The mallet: a physical object, not refined with N.
4545 let mallet: MalletParams = { M: 1, K: 1000 };
@@ -152,6 +152,7 @@ const relErr = () => (E0 > 0 ? (energy(P, s) + s.W - E0) / E0 : 0);
152152 const paramInputs = {
153153 m: document.getElementById("p-m") as HTMLInputElement,
154154 kappa: document.getElementById("p-kappa") as HTMLInputElement,
155+ tension: document.getElementById("p-tension") as HTMLInputElement,
155156 thick: document.getElementById("p-thick") as HTMLInputElement,
156157 a: document.getElementById("p-a") as HTMLInputElement,
157158 tau: document.getElementById("p-tau") as HTMLInputElement,
@@ -179,17 +180,24 @@ function distinctModes(): { w: number; count: number; k: number }[] {
179180 }
180181 return out;
181182 }
182-// Values that follow from the parameters: the lowest modes, with their (m,
183-// n) and the continuum's ratio (m^2 + n^2) / 2 when simply supported, and
184-// the decay times.
183+// Values that follow from the parameters: the lowest modes, with, when
184+// simply supported, their (m, n) and the continuum's ratio (from omega^2 ~ T
185+// k^2 + kappa k^4 with k^2 ~ m^2 + n^2: (m^2 + n^2) / 2 without tension,
186+// sqrt((m^2 + n^2) / 2) with tension alone), and the decay times.
185187 function showDerived() {
186188 const w1 = modes.omega[0], hz = (w: number) => w / (2 * Math.PI);
187- const free = P.edges === "free";
189+ const free = P.edges !== "simply"; // no (m, n)
190+ const kk = (Math.PI / (P.N + 1)) ** 2, T = P.edges === "free" ? 0 : P.tension;
191+ const cont = (m: number, n: number) => {
192+ const l = kk * (m * m + n * n);
193+ return Math.sqrt(T * l + P.kappa * l * l);
194+ };
188195 const rows = distinctModes().slice(0, 8).map(({ w, count, k }) => {
189196 const m = modes.m[k], n = modes.n[k];
190197 return (
191198 ` ${free ? "" : `(${m}, ${n})`}${count > 1 ? ` ×${count} ` : " "}${hz(w).toPrecision(3).padStart(7)} Hz ` +
192- `ratio ${(w / w1).toFixed(3)}` + (free ? "" : ` continuum ${((m * m + n * n) / 2).toFixed(1)}`)
199+ `ratio ${(w / w1).toFixed(3)}` +
200+ (free ? "" : ` continuum ${(cont(m, n) / cont(1, 1)).toFixed(3)}`)
193201 );
194202 });
195203 const g = (x: number) => (x === 0 ? "0" : x.toPrecision(4));
@@ -202,7 +210,7 @@ function showDerived() {
202210 (P.N === REF_N ? "" : `at ${P.N} × ${P.N}: m ${g(P.m)} kg, κ ${g(P.kappa)} N/m, a ${g(P.a)} m\n`) +
203211 `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
204212 `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
205- `lowest modes${free ? " (free edges, rigid motions removed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
213+ `lowest modes${P.edges === "free" ? " (free edges, rigid motions removed)" : P.edges === "clamped" ? " (clamped, computed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
206214 `nonlinear scale ${Number.isFinite(nonlinearScale(P, modes)) ? fmtEnergy(nonlinearScale(P, modes)) : "none (linear)"} (${free ? "lowest mode" : "mode 1, 1"}: stretching = bending energy)\n` +
207215 `decay time ${decay}\n` +
208216 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
@@ -306,7 +314,7 @@ function renderNote() {
306314 maxLength: 2 * predictedNoteLength(), signal: sig.key, x, y,
307315 };
308316 noteInfo = `${sig.label} (${x}, ${y}) m`;
309- noteFile = `plate-ss-N${P.N}-${strikeLabel}-${energy(P, lastStrike).toPrecision(3)}J-${sig.key}x${x}y${y}m-${pitch().toFixed(0)}Hz.wav`;
317+ noteFile = `plate-${P.edges}-N${P.N}-${strikeLabel}-${energy(P, lastStrike).toPrecision(3)}J-${sig.key}x${x}y${y}m-${pitch().toFixed(0)}Hz.wav`;
310318 soundStatus.textContent = "computing…";
311319 soundHint.hidden = true;
312320 player.stop();
@@ -343,6 +351,7 @@ function setParams(next: Params, N: number) {
343351 base = next;
344352 P = refine(base, next.edges === "free" ? Math.max(3, N) : N); // free: at least 3 per side
345353 edgesSelect.value = P.edges;
354+ paramInputs.tension.disabled = P.edges === "free";
346355 modes = new Modes(P);
347356 H = stepFor(P, modes, 0);
348357 s = newState(P.N);
@@ -376,10 +385,10 @@ edgesSelect.addEventListener("change", () => setParams({ ...base, edges: edgesSe
376385 document.getElementById("param-apply")!.addEventListener("click", () => {
377386 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
378387 const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
379- const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick"].includes(key) && x === 0))) &&
388+ const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick", "tension"].includes(key) && x === 0))) &&
380389 Object.values(mv).every((x) => Number.isFinite(x) && x > 0);
381390 if (!ok) {
382- paramStatus.textContent = "parameters must be positive (h, τ, τ top may be 0)";
391+ paramStatus.textContent = "parameters must be positive (T, h, τ, τ top may be 0)";
383392 return;
384393 }
385394 paramStatus.textContent = "";
src/plate/plate.tsmodified+67−29View file
@@ -1,10 +1,11 @@
11 // A square plate as N x N equal masses on vertical rails at the nodes of a
22 // grid, spacing a, moving up and down only, by u (m), stored row by row:
3-// node (i, j) at index i + N j. SI units. Two kinds of edge (see "Free
4-// edges" below for the second):
3+// node (i, j) at index i + N j. SI units. Three kinds of edge (see
4+// "Clamped" and "Free edges" below):
55 // - "simply" (simply supported): node (i, j) at x = (i + 1) a, y = (j + 1)
66 // a; the edge nodes (i or j = -1 or N) are supports at u = 0 (held, free
77 // to turn), the plate's side (N + 1) a.
8+// - "clamped": the same layout, the edge held and kept level.
89 // - "free": nothing holds the plate; node (i, j) at x = i a, y = j a, the
910 // edge nodes on the plate's edge, the side (N - 1) a. Each node carries
1011 // the plate's mass over its share of the surrounding cells: m inside, m /
@@ -37,6 +38,28 @@
3738 // would need another solver. The stretching is set by the thickness h:
3839 // ks / kappa = 12 (1 - nu^2) / h^2 (from D = Eh^3 / 12 (1 - nu^2), nu = 0.3),
3940 // so it matters once the plate moves by about h.
41+// Clamped: the supports also hold the slope at zero, an even ghost beyond
42+// them (u(-2, j) = u(0, j)). The support nodes then have a bending term,
43+// (L u) = -2 u(0, j) there, weighted 1/2 (the support's half cell): the
44+// energy becomes kappa / 2 (|L u|^2 + 2 sum c u^2), c the number of edges a
45+// node touches, i.e. the standard clamped 13-point stencil (21 at the
46+// centre next to an edge, 22 in a corner; the same operator as the free
47+// plate's stress function below, checked there against a continuum
48+// solution at second order). The Gaussian term integrates to zero on a
49+// held edge, so (Lap w)^2 is the whole bending energy, as when simply
50+// supported. The modes are computed (by symmetry blocks, as for the free
51+// plate); stretching as for the simply supported plate (same in-plane
52+// edge).
53+// Tension: a membrane tension T (N/m, force per unit length) on a held
54+// plate (simply supported or clamped), as on a drumhead: springs between
55+// neighbours and to the supports, energy T / 2 (u^T L u) (from T / 2 |grad
56+// w|^2; T a^2 / a^2 per spring, so T needs no scaling with the spacing),
57+// force -T L u, on top of the bending. Simply supported, L is diagonal in
58+// the same sine modes: omega_mn^2 = (T l + kappa l^2) / m, l = l_m + l_n.
59+// With the tension dominant the frequencies go as sqrt(m^2 + n^2), a square
60+// membrane; the bending stiffens the upper modes, as a real drumhead's does.
61+// Stretching then adds to the tension as the plate moves (tension
62+// modulation: a hard strike starts sharp). A free plate has no tension.
4063 // Free edges: the bending energy then needs its full form, D / 2 ((Lap w)^2
4164 // - 2 (1 - nu)(w_xx w_yy - w_xy^2)), whose second (Gaussian) term only
4265 // matters at a free edge and sets its natural conditions (zero moment and
@@ -77,13 +100,14 @@
77100
78101 import { symmetricEigen } from "../common/eigen";
79102
80-export type Edges = "simply" | "free";
103+export type Edges = "simply" | "clamped" | "free";
81104
82105 export interface Params {
83106 N: number; // masses per side (N^2 in all)
84107 edges: Edges;
85108 m: number; // kg, each mass
86109 kappa: number; // N/m, bending
110+ tension: number; // N/m, membrane tension (held edges only); 0 = none
87111 thick: number; // m, thickness, sets the stretching (von Karman); 0 = linear
88112 a: number; // m, rest spacing (for positions and drawing; it does not enter the equations)
89113 tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
@@ -114,6 +138,11 @@ export const isDamped = (P: Params) => P.tau > 0 || P.tauTop > 0;
114138 // Eigenvalues of the one-dimensional second difference, l_m, m = 1..N.
115139 const lam1 = (N: number, m: number) => 4 * Math.sin((m * Math.PI) / (2 * (N + 1))) ** 2;
116140 const isFree = (P: Params) => P.edges === "free";
141+const isClamped = (P: Params) => P.edges === "clamped";
142+// Membrane tension (N/m): 0 for a free plate.
143+export const tensionOf = (P: Params) => (isFree(P) ? 0 : P.tension);
144+// Edges a node touches (0, 1 or 2), for the clamped term.
145+const touched = (N: number, i: number, j: number) => (i === 0 ? 1 : 0) + (i === N - 1 ? 1 : 0) + (j === 0 ? 1 : 0) + (j === N - 1 ? 1 : 0);
117146 // Largest eigenvalue of L: 2 l_N; free, 8 (of m M^-1 L, exactly: see laplacian).
118147 const laplacianMax = (P: Params) => (isFree(P) ? 8 : 2 * lam1(P.N, P.N));
119148
@@ -504,7 +533,12 @@ export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Ar
504533 }
505534 laplacian(P.N, u, c);
506535 laplacian(P.N, c, F);
507- for (let k = 0; k < F.length; k++) F[k] *= -P.kappa;
536+ const T = tensionOf(P);
537+ for (let k = 0; k < F.length; k++) F[k] = -P.kappa * F[k] - T * c[k];
538+ if (isClamped(P)) {
539+ const N = P.N;
540+ for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) F[k] -= 2 * P.kappa * touched(N, i, j) * u[k];
541+ }
508542 if (!linear) stretching(P, u, F);
509543 }
510544
@@ -514,17 +548,21 @@ export function kinetic(P: Params, s: State): number {
514548 for (let k = 0; k < s.p.length; k++) T += (s.p[k] * s.p[k]) / M[k];
515549 return T / 2;
516550 }
517-// Bending energy alone (the linear part) at u.
518-export function bendingEnergy(P: Params, u: Float64Array): number {
551+// Linear energy at u: bending, plus the tension.
552+export function linearEnergy(P: Params, u: Float64Array): number {
519553 if (isFree(P)) return freeBending(P, u);
520- const c = new Float64Array(u.length);
521- laplacian(P.N, u, c);
522- let V = 0;
523- for (const x of c) V += x * x;
524- return (P.kappa * V) / 2;
554+ const N = P.N, c = new Float64Array(u.length);
555+ laplacian(N, u, c);
556+ let V = 0, uLu = 0;
557+ for (let k = 0; k < u.length; k++) {
558+ V += c[k] * c[k];
559+ uLu += u[k] * c[k];
560+ }
561+ if (isClamped(P)) for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) V += 2 * touched(N, i, j) * u[k] * u[k];
562+ return (P.kappa * V + tensionOf(P) * uLu) / 2;
525563 }
526564 export const stretchEnergy = (P: Params, u: Float64Array) => stretching(P, u);
527-export const potential = (P: Params, s: State) => bendingEnergy(P, s.u) + stretching(P, s.u);
565+export const potential = (P: Params, s: State) => linearEnergy(P, s.u) + stretching(P, s.u);
528566 export const energy = (P: Params, s: State) => kinetic(P, s) + potential(P, s);
529567
530568 // Velocity Verlet, with drag and loss in half steps around it (Strang), as
@@ -605,8 +643,8 @@ export class Verlet {
605643
606644 // The normal modes, sorted by frequency: omega (rad/s), and for a simply
607645 // supported plate the mode numbers (m, n) (ties: by m, then n); for a free
608-// plate the displacement shapes S (normalized as S^T M S = 1, M the node
609-// masses), computed, and m = n = 0.
646+// or clamped plate the displacement shapes S (normalized as S^T M S = 1, M
647+// the node masses), computed, and m = n = 0.
610648 export class Modes {
611649 N: number;
612650 count: number;
@@ -617,8 +655,8 @@ export class Modes {
617655 constructor(P: Params) {
618656 const N = P.N;
619657 this.N = N;
620- if (isFree(P)) {
621- const { lambda, S } = freeModes(P);
658+ if (isFree(P) || isClamped(P)) {
659+ const { lambda, S } = computedModes(P);
622660 this.count = lambda.length;
623661 this.omega = Float64Array.from(lambda, (l) => Math.sqrt(Math.max(l, 0)));
624662 this.S = S;
@@ -626,9 +664,9 @@ export class Modes {
626664 this.n = new Int32Array(this.count);
627665 return;
628666 }
629- const w0 = Math.sqrt(P.kappa / P.m);
667+ const T = tensionOf(P), w = (l: number) => Math.sqrt((T * l + P.kappa * l * l) / P.m);
630668 const all: [number, number, number][] = [];
631- for (let m = 1; m <= N; m++) for (let n = 1; n <= N; n++) all.push([w0 * (lam1(N, m) + lam1(N, n)), m, n]);
669+ for (let m = 1; m <= N; m++) for (let n = 1; n <= N; n++) all.push([w(lam1(N, m) + lam1(N, n)), m, n]);
632670 all.sort((x, y) => x[0] - y[0] || x[1] - y[1] || x[2] - y[2]);
633671 this.count = all.length;
634672 this.omega = Float64Array.from(all, (x) => x[0]);
@@ -640,13 +678,13 @@ export class Modes {
640678 }
641679 }
642680
643-// Modes of a free plate: eigenvalues lambda = omega^2 of M^-1/2 K M^-1/2 (K
644-// the bending stiffness, M the node masses) and displacement shapes, without
645-// the three rigid motions (lambda = 0). The square is symmetric under x ->
681+// Computed modes (free or clamped plate): eigenvalues lambda = omega^2 of M^-1/2
682+// K M^-1/2 (K the linear stiffness, M the node masses) and displacement
683+// shapes, without a free plate's three rigid motions (lambda = 0). The square is symmetric under x ->
646684 // -x and y -> -y (masses included), so the matrix splits into four blocks,
647685 // even or odd in x and in y (built on mirror pairs of nodes), each a quarter
648686 // the size: 1/16 of the work of one full eigensolve.
649-function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
687+function computedModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
650688 const N = P.N, n = N * N;
651689 // One-dimensional symmetric (even) or antisymmetric (odd) unit vectors on mirror pairs.
652690 const basis1 = (odd: boolean): [number, number][][] => {
@@ -672,7 +710,7 @@ function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
672710 for (let q = 0; q < d; q++) {
673711 e.fill(0);
674712 for (const [k, w] of vecs[q]) e[k] = w * rs[k];
675- freeBending(P, e, F); // F = -K M^-1/2 b_q
713+ forces(P, e, F, undefined, true); // F = -K M^-1/2 b_q
676714 for (let r = 0; r < d; r++) {
677715 let acc = 0;
678716 for (const [k, w] of vecs[r]) acc -= w * rs[k] * F[k];
@@ -691,22 +729,22 @@ function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
691729 }
692730 }
693731 found.sort((x, y) => x.l - y.l);
694- const kept = found.slice(3); // the rigid motions, lambda = 0 to round-off
732+ const kept = isFree(P) ? found.slice(3) : found; // a free plate's rigid motions, lambda = 0 to round-off
695733 return { lambda: Float64Array.from(kept, (x) => x.l), S: kept.map((x) => x.v) };
696734 }
697735
698736 // Amplitude decay time (s) of each mode, in Modes order (Infinity if
699737 // undamped): 1 / (1 / tau + c q / 2), q = S^T L S for the mode's shape
700738 // normalized by S^T M S = 1: (l_m + l_n) / m when simply supported (exact:
701-// L is diagonal in the modes); free, from the computed shape (L is not
702-// diagonal in the free modes, so this is approximate, as for the beam).
739+// L is diagonal in the modes); free or clamped, from the computed shape (L
740+// is not diagonal in those modes, so this is approximate, as for the beam).
703741 export function decayTimes(P: Params, modes: Modes): Float64Array {
704742 const c = lossCoefficient(P), Ls = new Float64Array(P.N * P.N);
705743 return Float64Array.from(modes.omega, (_, k) => {
706744 let q = 0;
707745 if (modes.S) {
708746 const v = modes.S[k];
709- laplacian(P.N, v, Ls, true);
747+ laplacian(P.N, v, Ls, isFree(P));
710748 for (let i = 0; i < v.length; i++) q += v[i] * Ls[i];
711749 } else q = (lam1(P.N, modes.m[k]) + lam1(P.N, modes.n[k])) / P.m;
712750 const rate = (P.tau > 0 ? 1 / P.tau : 0) + (c * q) / 2;
@@ -715,7 +753,7 @@ export function decayTimes(P: Params, modes: Modes): Float64Array {
715753 }
716754
717755 // Nonlinear scale (J): the energy of the lowest mode ((1, 1) when simply
718-// supported) at which its stretching energy equals its bending energy
756+// supported) at which its stretching energy equals its linear energy
719757 // (Infinity when linear). Bending goes as the amplitude squared, stretching
720758 // as its fourth power.
721759 export function nonlinearScale(P: Params, modes: Modes): number {
@@ -724,7 +762,7 @@ export function nonlinearScale(P: Params, modes: Modes): number {
724762 let u = new Float64Array(N * N);
725763 if (modes.S) u = Float64Array.from(modes.S[0]);
726764 else for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) u[i + N * j] = Math.sin((Math.PI * (i + 1)) / (N + 1)) * Math.sin((Math.PI * (j + 1)) / (N + 1));
727- const V2 = bendingEnergy(P, u), V4 = stretching(P, u);
765+ const V2 = linearEnergy(P, u), V4 = stretching(P, u);
728766 return V4 > 0 ? (V2 * V2) / V4 : Infinity;
729767 }
730768