Free edges for both finite-difference methods
The free-free plate has no ghost-layer form, so plate_matrices.m assembles the
stiffness from the bending energy (finite-difference energy method) with the
free-edge conditions as natural boundary conditions; its null space is exactly
the three rigid-body motions, which plate_rigid_modes.m removes after every
step. Edge curvatures come from the moment condition (plate_curvature_free.m).
Checked against Leissa's free square plate frequencies; the existing CFL bound
still holds.
18 changed files+920−147
README.mdmodified+2−2View file
@@ -21,7 +21,7 @@ The app and the tutorials are peers over the same core and components.
2121 Three methods solve the same problem specification, and the app can run any two and compare them directly:
2222
2323 - **Modal superposition** (`plate_modal.m`): closed-form modes of the simply supported rectangle, each integrated exactly in time for piecewise-linear forcing. No time-step limit; the approximation is the modal truncation.
24-- **Explicit finite differences** (`plate_fd.m`): second-order scheme with the 13-point biharmonic stencil, simply supported or clamped edges. Stable only below the CFL bound from `plate_fd_plan.m`; the slider goes past it on purpose.
24+- **Explicit finite differences** (`plate_fd.m`): second-order scheme with the 13-point biharmonic stencil; simply supported, clamped or free edges (free edges use the energy-method operators of `plate_matrices.m`). Stable only below the CFL bound from `plate_fd_plan.m`; the slider goes past it on purpose.
2525 - **Implicit finite differences** (`plate_implicit.m`): the same stencils assembled as sparse matrices (`plate_matrices.m`), stiffness averaged over three time levels (theta scheme, theta = 1/4 is the trapezoidal rule). One sparse Cholesky factorization, two triangular solves per step, no stability limit; the step is chosen for accuracy. Needs numbl 0.4.21 or later for its sparse direct solvers.
2626
2727 A method is a data record (`packages/core/src/methods/*.ts`) plus a `.m` file. Adding one means adding both; nothing else in core changes.
@@ -53,7 +53,7 @@ npx numbl run --path packages/core/src/m my_script.m
5353 ## Honest limits
5454
5555 - Kirchhoff-Love thin-plate theory: no shear deformation or rotary inertia, so thick plates and very high modes are wrong.
56-- Edge conditions are idealizations; free edges are not yet supported.
56+- Edge conditions are idealizations. Free edges are handled by a discrete bending energy whose natural boundary conditions are first-order accurate at the edge, and the rigid-body motion of a free plate is removed rather than modelled, as the supports of a real free-free rig would do.
5757 - Damping is a two-parameter model, chosen, not measured.
5858 - Stress is a second derivative of the displacement and converges more slowly than displacement in both methods; near the strike it depends on the load footprint.
5959 - Radiated sound is Rayleigh's integral in an infinite baffle: an estimate of the waveform, not a calibrated level.
apps/tutorials/src/tutorials/cfl.mdxmodified+1−1View file
@@ -59,7 +59,7 @@ the output of the planning function for four grids on the same plate.
5959 The lab shows an estimate of this cost before each finite-difference run. For a
6060 long simulated time on a fine grid the modal method is the practical choice on
6161 a simply supported plate; the finite-difference method earns its cost where no
62-closed-form modes exist, such as clamped edges.
62+closed-form modes exist, such as clamped or free edges.
6363
6464 ## Going further
6565
packages/components/src/params/ParameterPanel.tsxmodified+379−58View file
@@ -56,7 +56,8 @@ export interface ParameterPanelProps {
5656 export function ParameterPanel(props: ParameterPanelProps) {
5757 const { spec, onChange, sections = ALL_SECTIONS, method, issues = [], fdDtMax, open } = props;
5858 const set = (patch: Partial<ProblemSpec>) => onChange({ ...spec, ...patch });
59- const issue = (path: string) => issues.find((i) => i.path === path || i.path.startsWith(path + "."));
59+ const issue = (path: string) =>
60+ issues.find((i) => i.path === path || i.path.startsWith(path + "."));
6061 const invalid = (path: string) => issue(path)?.severity === "error";
6162 const show = (s: ParameterSection) => sections.includes(s);
6263 const isOpen = (s: ParameterSection) => (open ? open.includes(s) : true);
@@ -66,7 +67,8 @@ export function ParameterPanel(props: ParameterPanelProps) {
6667 (p) => p.material.name === spec.material.name && p.material.E === spec.material.E,
6768 )?.id;
6869 const platePresetId = PLATE_PRESETS.find(
69- (p) => p.plate.Lx === spec.plate.Lx && p.plate.Ly === spec.plate.Ly && p.plate.h === spec.plate.h,
70+ (p) =>
71+ p.plate.Lx === spec.plate.Lx && p.plate.Ly === spec.plate.Ly && p.plate.h === spec.plate.h,
7072 )?.id;
7173
7274 return (
@@ -91,10 +93,36 @@ export function ParameterPanel(props: ParameterPanelProps) {
9193 ))}
9294 </select>
9395 </div>
94- {platePresetId && <div className="pl-small pl-muted">{PLATE_PRESETS.find((p) => p.id === platePresetId)?.note}</div>}
95- <NumberField label="Length Lx" value={spec.plate.Lx} scale={1000} unit="mm" invalid={invalid("plate.Lx")} onChange={(v) => set({ plate: { ...spec.plate, Lx: v } })} />
96- <NumberField label="Width Ly" value={spec.plate.Ly} scale={1000} unit="mm" invalid={invalid("plate.Ly")} onChange={(v) => set({ plate: { ...spec.plate, Ly: v } })} />
97- <NumberField label="Thickness h" value={spec.plate.h} scale={1000} unit="mm" step={0.1} invalid={invalid("plate.h")} onChange={(v) => set({ plate: { ...spec.plate, h: v } })} />
96+ {platePresetId && (
97+ <div className="pl-small pl-muted">
98+ {PLATE_PRESETS.find((p) => p.id === platePresetId)?.note}
99+ </div>
100+ )}
101+ <NumberField
102+ label="Length Lx"
103+ value={spec.plate.Lx}
104+ scale={1000}
105+ unit="mm"
106+ invalid={invalid("plate.Lx")}
107+ onChange={(v) => set({ plate: { ...spec.plate, Lx: v } })}
108+ />
109+ <NumberField
110+ label="Width Ly"
111+ value={spec.plate.Ly}
112+ scale={1000}
113+ unit="mm"
114+ invalid={invalid("plate.Ly")}
115+ onChange={(v) => set({ plate: { ...spec.plate, Ly: v } })}
116+ />
117+ <NumberField
118+ label="Thickness h"
119+ value={spec.plate.h}
120+ scale={1000}
121+ unit="mm"
122+ step={0.1}
123+ invalid={invalid("plate.h")}
124+ onChange={(v) => set({ plate: { ...spec.plate, h: v } })}
125+ />
98126 </Group>
99127 )}
100128 {show("material") && (
@@ -118,33 +146,89 @@ export function ParameterPanel(props: ParameterPanelProps) {
118146 </select>
119147 </div>
120148 {materialPresetId && (
121- <div className="pl-small pl-muted">{MATERIAL_PRESETS.find((p) => p.id === materialPresetId)?.source} Handbook values, not a certificate.</div>
149+ <div className="pl-small pl-muted">
150+ {MATERIAL_PRESETS.find((p) => p.id === materialPresetId)?.source} Handbook values, not
151+ a certificate.
152+ </div>
122153 )}
123- <NumberField label="Young's modulus E" value={spec.material.E} scale={1e-9} unit="GPa" invalid={invalid("material.E")} onChange={(v) => set({ material: { ...spec.material, E: v, name: customName(spec.material.name) } })} />
124- <NumberField label="Poisson's ratio ν" value={spec.material.nu} step={0.01} invalid={invalid("material.nu")} onChange={(v) => set({ material: { ...spec.material, nu: v, name: customName(spec.material.name) } })} />
125- <NumberField label="Density ρ" value={spec.material.rho} unit="kg/m³" invalid={invalid("material.rho")} onChange={(v) => set({ material: { ...spec.material, rho: v, name: customName(spec.material.name) } })} />
154+ <NumberField
155+ label="Young's modulus E"
156+ value={spec.material.E}
157+ scale={1e-9}
158+ unit="GPa"
159+ invalid={invalid("material.E")}
160+ onChange={(v) =>
161+ set({ material: { ...spec.material, E: v, name: customName(spec.material.name) } })
162+ }
163+ />
164+ <NumberField
165+ label="Poisson's ratio ν"
166+ value={spec.material.nu}
167+ step={0.01}
168+ invalid={invalid("material.nu")}
169+ onChange={(v) =>
170+ set({ material: { ...spec.material, nu: v, name: customName(spec.material.name) } })
171+ }
172+ />
173+ <NumberField
174+ label="Density ρ"
175+ value={spec.material.rho}
176+ unit="kg/m³"
177+ invalid={invalid("material.rho")}
178+ onChange={(v) =>
179+ set({ material: { ...spec.material, rho: v, name: customName(spec.material.name) } })
180+ }
181+ />
126182 </Group>
127183 )}
128184 {show("edges") && (
129185 <Group title="Edges" issues={issues} prefix="bc" open={isOpen("edges")}>
130186 <div className="pl-field pl-field-wide">
131187 <label>Condition</label>
132- <select className="pl-select" value={spec.bc} onChange={(ev) => set({ bc: ev.target.value as BoundaryCondition })}>
188+ <select
189+ className="pl-select"
190+ value={spec.bc}
191+ onChange={(ev) => set({ bc: ev.target.value as BoundaryCondition })}
192+ >
133193 <option value="simply-supported">simply supported (w = 0, no moment)</option>
134194 <option value="clamped">clamped (w = 0, no slope)</option>
195+ <option value="free">free (no moment, no shear): the free-free plate</option>
135196 </select>
136197 </div>
137198 {method && !method.supports.bc.includes(spec.bc) && (
138- <div className="pl-issue">{method.name} has no {spec.bc} solution; only the finite-difference method can run this.</div>
199+ <div className="pl-issue">
200+ {method.name} has no {spec.bc} solution; the finite-difference methods can run this.
201+ </div>
139202 )}
140- <div className="pl-small pl-muted">Both are idealizations: a real clamp has finite stiffness and a real support some rotational restraint.</div>
203+ <div className="pl-small pl-muted">
204+ All three are idealizations: a real clamp has finite stiffness, a real support some
205+ rotational restraint, and a free plate rests on something that takes its rigid-body
206+ motion, which the solvers remove.
207+ </div>
141208 </Group>
142209 )}
143210 {show("damping") && (
144211 <Group title="Damping" issues={issues} prefix="damping" open={isOpen("damping")}>
145- <NumberField label="σ₀ (all frequencies)" value={spec.damping.sigma0} unit="1/s" step={0.1} invalid={invalid("damping")} onChange={(v) => set({ damping: { ...spec.damping, sigma0: v } })} />
146- <NumberField label="σ₁ (rises with frequency)" value={spec.damping.sigma1} unit="m²/s" step={0.001} invalid={invalid("damping")} onChange={(v) => set({ damping: { ...spec.damping, sigma1: v } })} />
147- <div className="pl-small pl-muted">A mode with Laplacian eigenvalue λ decays at rate σ₀ + σ₁λ. Illustrative values; real losses depend on the mounting.</div>
212+ <NumberField
213+ label="σ₀ (all frequencies)"
214+ value={spec.damping.sigma0}
215+ unit="1/s"
216+ step={0.1}
217+ invalid={invalid("damping")}
218+ onChange={(v) => set({ damping: { ...spec.damping, sigma0: v } })}
219+ />
220+ <NumberField
221+ label="σ₁ (rises with frequency)"
222+ value={spec.damping.sigma1}
223+ unit="m²/s"
224+ step={0.001}
225+ invalid={invalid("damping")}
226+ onChange={(v) => set({ damping: { ...spec.damping, sigma1: v } })}
227+ />
228+ <div className="pl-small pl-muted">
229+ A mode with Laplacian eigenvalue λ decays at rate σ₀ + σ₁λ. Illustrative values; real
230+ losses depend on the mounting.
231+ </div>
148232 </Group>
149233 )}
150234 {show("excitation") && (
@@ -159,7 +243,14 @@ export function ParameterPanel(props: ParameterPanelProps) {
159243 const next: Excitation =
160244 type === "impulse"
161245 ? { type, x: e.x, y: e.y, force: e.force, duration: 1e-3, radius: e.radius }
162- : { type, x: e.x, y: e.y, force: Math.min(e.force, 5), frequency: 500, radius: e.radius };
246+ : {
247+ type,
248+ x: e.x,
249+ y: e.y,
250+ force: Math.min(e.force, 5),
251+ frequency: 500,
252+ radius: e.radius,
253+ };
163254 set({ excitation: next });
164255 }}
165256 >
@@ -167,66 +258,276 @@ export function ParameterPanel(props: ParameterPanelProps) {
167258 <option value="harmonic">harmonic point force</option>
168259 </select>
169260 </div>
170- <NumberField label="Strike x" value={e.x} scale={1000} unit="mm" invalid={invalid("excitation")} onChange={(v) => set({ excitation: { ...e, x: v } })} />
171- <NumberField label="Strike y" value={e.y} scale={1000} unit="mm" invalid={invalid("excitation")} onChange={(v) => set({ excitation: { ...e, y: v } })} />
172- <NumberField label={e.type === "impulse" ? "Peak force" : "Force amplitude"} value={e.force} unit="N" invalid={invalid("excitation.force")} onChange={(v) => set({ excitation: { ...e, force: v } })} />
261+ <NumberField
262+ label="Strike x"
263+ value={e.x}
264+ scale={1000}
265+ unit="mm"
266+ invalid={invalid("excitation")}
267+ onChange={(v) => set({ excitation: { ...e, x: v } })}
268+ />
269+ <NumberField
270+ label="Strike y"
271+ value={e.y}
272+ scale={1000}
273+ unit="mm"
274+ invalid={invalid("excitation")}
275+ onChange={(v) => set({ excitation: { ...e, y: v } })}
276+ />
277+ <NumberField
278+ label={e.type === "impulse" ? "Peak force" : "Force amplitude"}
279+ value={e.force}
280+ unit="N"
281+ invalid={invalid("excitation.force")}
282+ onChange={(v) => set({ excitation: { ...e, force: v } })}
283+ />
173284 {e.type === "impulse" ? (
174- <NumberField label="Contact time" value={e.duration} scale={1000} unit="ms" step={0.1} invalid={invalid("excitation.duration")} onChange={(v) => set({ excitation: { ...e, duration: v } })} />
285+ <NumberField
286+ label="Contact time"
287+ value={e.duration}
288+ scale={1000}
289+ unit="ms"
290+ step={0.1}
291+ invalid={invalid("excitation.duration")}
292+ onChange={(v) => set({ excitation: { ...e, duration: v } })}
293+ />
175294 ) : (
176- <NumberField label="Frequency" value={e.frequency} unit="Hz" invalid={invalid("excitation.frequency")} onChange={(v) => set({ excitation: { ...e, frequency: v } })} />
295+ <NumberField
296+ label="Frequency"
297+ value={e.frequency}
298+ unit="Hz"
299+ invalid={invalid("excitation.frequency")}
300+ onChange={(v) => set({ excitation: { ...e, frequency: v } })}
301+ />
177302 )}
178- <NumberField label="Footprint radius" value={e.radius} scale={1000} unit="mm" step={1} invalid={invalid("excitation.radius")} onChange={(v) => set({ excitation: { ...e, radius: v } })} />
179- <div className="pl-small pl-muted">A finite footprint stands in for the contact patch; a point load would give infinite stress under the hammer.</div>
303+ <NumberField
304+ label="Footprint radius"
305+ value={e.radius}
306+ scale={1000}
307+ unit="mm"
308+ step={1}
309+ invalid={invalid("excitation.radius")}
310+ onChange={(v) => set({ excitation: { ...e, radius: v } })}
311+ />
312+ <div className="pl-small pl-muted">
313+ A finite footprint stands in for the contact patch; a point load would give infinite
314+ stress under the hammer.
315+ </div>
180316 </Group>
181317 )}
182318 {show("probe") && (
183319 <Group title="Probe" issues={issues} prefix="probe" open={isOpen("probe")}>
184- <NumberField label="x" value={spec.probe.x} scale={1000} unit="mm" invalid={invalid("probe")} onChange={(v) => set({ probe: { ...spec.probe, x: v } })} />
185- <NumberField label="y" value={spec.probe.y} scale={1000} unit="mm" invalid={invalid("probe")} onChange={(v) => set({ probe: { ...spec.probe, y: v } })} />
186- <div className="pl-small pl-muted">Time series and spectra are recorded here. Finite differences snap it to the nearest grid node.</div>
320+ <NumberField
321+ label="x"
322+ value={spec.probe.x}
323+ scale={1000}
324+ unit="mm"
325+ invalid={invalid("probe")}
326+ onChange={(v) => set({ probe: { ...spec.probe, x: v } })}
327+ />
328+ <NumberField
329+ label="y"
330+ value={spec.probe.y}
331+ scale={1000}
332+ unit="mm"
333+ invalid={invalid("probe")}
334+ onChange={(v) => set({ probe: { ...spec.probe, y: v } })}
335+ />
336+ <div className="pl-small pl-muted">
337+ Time series and spectra are recorded here. Finite differences snap it to the nearest
338+ grid node.
339+ </div>
187340 </Group>
188341 )}
189342 {show("listener") && (
190343 <Group title="Listener (sound)" issues={issues} prefix="listener" open={isOpen("listener")}>
191- <NumberField label="x" value={spec.listener.x} scale={1000} unit="mm" onChange={(v) => set({ listener: { ...spec.listener, x: v } })} />
192- <NumberField label="y" value={spec.listener.y} scale={1000} unit="mm" onChange={(v) => set({ listener: { ...spec.listener, y: v } })} />
193- <NumberField label="Height z" value={spec.listener.z} unit="m" step={0.1} onChange={(v) => set({ listener: { ...spec.listener, z: v } })} />
194- <div className="pl-small pl-muted">Rayleigh integral in an infinite baffle: an estimate of the waveform, not a calibrated level.</div>
344+ <NumberField
345+ label="x"
346+ value={spec.listener.x}
347+ scale={1000}
348+ unit="mm"
349+ onChange={(v) => set({ listener: { ...spec.listener, x: v } })}
350+ />
351+ <NumberField
352+ label="y"
353+ value={spec.listener.y}
354+ scale={1000}
355+ unit="mm"
356+ onChange={(v) => set({ listener: { ...spec.listener, y: v } })}
357+ />
358+ <NumberField
359+ label="Height z"
360+ value={spec.listener.z}
361+ unit="m"
362+ step={0.1}
363+ onChange={(v) => set({ listener: { ...spec.listener, z: v } })}
364+ />
365+ <div className="pl-small pl-muted">
366+ Rayleigh integral in an infinite baffle: an estimate of the waveform, not a calibrated
367+ level.
368+ </div>
195369 </Group>
196370 )}
197371 {show("output") && (
198372 <Group title="Output" issues={issues} prefix="output" open={isOpen("output")}>
199- <NumberField label="Sample rate" value={spec.output.fs} unit="Hz" integer invalid={invalid("output.fs")} onChange={(v) => set({ output: { ...spec.output, fs: v } })} />
200- <NumberField label="Duration (modal, audio)" value={spec.output.duration} unit="s" step={0.1} invalid={invalid("output.duration")} onChange={(v) => set({ output: { ...spec.output, duration: v } })} />
201- <NumberField label="Animation frames" value={spec.output.nframes} integer invalid={invalid("output.nframes")} onChange={(v) => set({ output: { ...spec.output, nframes: v } })} />
202- <NumberField label="Display grid nx" value={spec.output.grid[0]} integer onChange={(v) => set({ output: { ...spec.output, grid: [v, spec.output.grid[1]] } })} />
203- <NumberField label="Display grid ny" value={spec.output.grid[1]} integer onChange={(v) => set({ output: { ...spec.output, grid: [spec.output.grid[0], v] } })} />
373+ <NumberField
374+ label="Sample rate"
375+ value={spec.output.fs}
376+ unit="Hz"
377+ integer
378+ invalid={invalid("output.fs")}
379+ onChange={(v) => set({ output: { ...spec.output, fs: v } })}
380+ />
381+ <NumberField
382+ label="Duration (modal, audio)"
383+ value={spec.output.duration}
384+ unit="s"
385+ step={0.1}
386+ invalid={invalid("output.duration")}
387+ onChange={(v) => set({ output: { ...spec.output, duration: v } })}
388+ />
389+ <NumberField
390+ label="Animation frames"
391+ value={spec.output.nframes}
392+ integer
393+ invalid={invalid("output.nframes")}
394+ onChange={(v) => set({ output: { ...spec.output, nframes: v } })}
395+ />
396+ <NumberField
397+ label="Display grid nx"
398+ value={spec.output.grid[0]}
399+ integer
400+ onChange={(v) => set({ output: { ...spec.output, grid: [v, spec.output.grid[1]] } })}
401+ />
402+ <NumberField
403+ label="Display grid ny"
404+ value={spec.output.grid[1]}
405+ integer
406+ onChange={(v) => set({ output: { ...spec.output, grid: [spec.output.grid[0], v] } })}
407+ />
204408 </Group>
205409 )}
206410 {show("modal") && (
207- <Group title="Modal method" issues={issues} prefix="modal" open={isOpen("modal")} disabled={method ? method.settingsKey !== "modal" : false}>
208- <NumberField label="Max frequency" value={spec.modal.fmax} unit="Hz" invalid={invalid("modal.fmax")} onChange={(v) => set({ modal: { ...spec.modal, fmax: v } })} />
209- <NumberField label="Max modes" value={spec.modal.max_modes} integer onChange={(v) => set({ modal: { ...spec.modal, max_modes: v } })} />
210- <div className="pl-small pl-muted">Truncation of the modal sum: displacement converges fast, stress near the strike slowly.</div>
411+ <Group
412+ title="Modal method"
413+ issues={issues}
414+ prefix="modal"
415+ open={isOpen("modal")}
416+ disabled={method ? method.settingsKey !== "modal" : false}
417+ >
418+ <NumberField
419+ label="Max frequency"
420+ value={spec.modal.fmax}
421+ unit="Hz"
422+ invalid={invalid("modal.fmax")}
423+ onChange={(v) => set({ modal: { ...spec.modal, fmax: v } })}
424+ />
425+ <NumberField
426+ label="Max modes"
427+ value={spec.modal.max_modes}
428+ integer
429+ onChange={(v) => set({ modal: { ...spec.modal, max_modes: v } })}
430+ />
431+ <div className="pl-small pl-muted">
432+ Truncation of the modal sum: displacement converges fast, stress near the strike slowly.
433+ </div>
211434 </Group>
212435 )}
213436 {show("fd") && (
214- <Group title="Explicit finite differences" issues={issues} prefix="fd" open={isOpen("fd")} disabled={method ? method.settingsKey !== "fd" : false}>
215- <NumberField label="Grid nx" value={spec.fd.nx} integer min={5} invalid={invalid("fd")} onChange={(v) => set({ fd: { ...spec.fd, nx: v } })} />
216- <NumberField label="Grid ny" value={spec.fd.ny} integer min={5} invalid={invalid("fd")} onChange={(v) => set({ fd: { ...spec.fd, ny: v } })} />
217- <NumberField label="Duration" value={spec.fd.duration} scale={1000} unit="ms" step={0.5} onChange={(v) => set({ fd: { ...spec.fd, duration: v } })} />
218- <div className="pl-small" style={{ marginTop: 6 }}>Time step</div>
219- <CflSlider value={spec.fd.cfl} dtMax={fdDtMax} onChange={(v) => set({ fd: { ...spec.fd, cfl: v } })} />
437+ <Group
438+ title="Explicit finite differences"
439+ issues={issues}
440+ prefix="fd"
441+ open={isOpen("fd")}
442+ disabled={method ? method.settingsKey !== "fd" : false}
443+ >
444+ <NumberField
445+ label="Grid nx"
446+ value={spec.fd.nx}
447+ integer
448+ min={5}
449+ invalid={invalid("fd")}
450+ onChange={(v) => set({ fd: { ...spec.fd, nx: v } })}
451+ />
452+ <NumberField
453+ label="Grid ny"
454+ value={spec.fd.ny}
455+ integer
456+ min={5}
457+ invalid={invalid("fd")}
458+ onChange={(v) => set({ fd: { ...spec.fd, ny: v } })}
459+ />
460+ <NumberField
461+ label="Duration"
462+ value={spec.fd.duration}
463+ scale={1000}
464+ unit="ms"
465+ step={0.5}
466+ onChange={(v) => set({ fd: { ...spec.fd, duration: v } })}
467+ />
468+ <div className="pl-small" style={{ marginTop: 6 }}>
469+ Time step
470+ </div>
471+ <CflSlider
472+ value={spec.fd.cfl}
473+ dtMax={fdDtMax}
474+ onChange={(v) => set({ fd: { ...spec.fd, cfl: v } })}
475+ />
220476 </Group>
221477 )}
222478 {show("implicit") && (
223- <Group title="Implicit finite differences" issues={issues} prefix="implicit" open={isOpen("implicit")} disabled={method ? method.settingsKey !== "implicit" : false}>
224- <NumberField label="Grid nx" value={spec.implicit.nx} integer min={5} invalid={invalid("implicit")} onChange={(v) => set({ implicit: { ...spec.implicit, nx: v } })} />
225- <NumberField label="Grid ny" value={spec.implicit.ny} integer min={5} invalid={invalid("implicit")} onChange={(v) => set({ implicit: { ...spec.implicit, ny: v } })} />
226- <NumberField label="Duration" value={spec.implicit.duration} scale={1000} unit="ms" step={5} onChange={(v) => set({ implicit: { ...spec.implicit, duration: v } })} />
227- <NumberField label="Steps per sample" value={spec.implicit.steps_per_sample} step={0.5} min={0.25} invalid={invalid("implicit.steps_per_sample")} onChange={(v) => set({ implicit: { ...spec.implicit, steps_per_sample: v } })} />
228- <NumberField label="θ (time averaging)" value={spec.implicit.theta} step={0.05} min={0} max={1} invalid={invalid("implicit.theta")} onChange={(v) => set({ implicit: { ...spec.implicit, theta: v } })} />
229- <div className="pl-small pl-muted">No stability limit for θ ≥ 1/4; the step sets accuracy instead. More steps per sample tune the high modes better and cost proportionally.</div>
479+ <Group
480+ title="Implicit finite differences"
481+ issues={issues}
482+ prefix="implicit"
483+ open={isOpen("implicit")}
484+ disabled={method ? method.settingsKey !== "implicit" : false}
485+ >
486+ <NumberField
487+ label="Grid nx"
488+ value={spec.implicit.nx}
489+ integer
490+ min={5}
491+ invalid={invalid("implicit")}
492+ onChange={(v) => set({ implicit: { ...spec.implicit, nx: v } })}
493+ />
494+ <NumberField
495+ label="Grid ny"
496+ value={spec.implicit.ny}
497+ integer
498+ min={5}
499+ invalid={invalid("implicit")}
500+ onChange={(v) => set({ implicit: { ...spec.implicit, ny: v } })}
501+ />
502+ <NumberField
503+ label="Duration"
504+ value={spec.implicit.duration}
505+ scale={1000}
506+ unit="ms"
507+ step={5}
508+ onChange={(v) => set({ implicit: { ...spec.implicit, duration: v } })}
509+ />
510+ <NumberField
511+ label="Steps per sample"
512+ value={spec.implicit.steps_per_sample}
513+ step={0.5}
514+ min={0.25}
515+ invalid={invalid("implicit.steps_per_sample")}
516+ onChange={(v) => set({ implicit: { ...spec.implicit, steps_per_sample: v } })}
517+ />
518+ <NumberField
519+ label="θ (time averaging)"
520+ value={spec.implicit.theta}
521+ step={0.05}
522+ min={0}
523+ max={1}
524+ invalid={invalid("implicit.theta")}
525+ onChange={(v) => set({ implicit: { ...spec.implicit, theta: v } })}
526+ />
527+ <div className="pl-small pl-muted">
528+ No stability limit for θ ≥ 1/4; the step sets accuracy instead. More steps per sample
529+ tune the high modes better and cost proportionally.
530+ </div>
230531 </Group>
231532 )}
232533 </div>
@@ -237,18 +538,38 @@ function customName(name: string): string {
237538 return name.endsWith(" (edited)") ? name : `${name} (edited)`;
238539 }
239540
240-function Group(props: { title: string; prefix: string; issues: ValidationIssue[]; open: boolean; disabled?: boolean; children: ReactNode }) {
241- const related = props.issues.filter((i) => i.path === props.prefix || i.path.startsWith(props.prefix + "."));
541+function Group(props: {
542+ title: string;
543+ prefix: string;
544+ issues: ValidationIssue[];
545+ open: boolean;
546+ disabled?: boolean;
547+ children: ReactNode;
548+}) {
549+ const related = props.issues.filter(
550+ (i) => i.path === props.prefix || i.path.startsWith(props.prefix + "."),
551+ );
242552 return (
243553 <details className={`pl-group${props.disabled ? " pl-disabled" : ""}`} open={props.open}>
244554 <summary>
245555 {props.title}
246- {related.some((i) => i.severity === "error") && <span className="pl-badge pl-badge-error" style={{ marginLeft: 6 }}>check</span>}
247- {props.disabled && <span className="pl-badge pl-badge-info" style={{ marginLeft: 6 }}>not used by this method</span>}
556+ {related.some((i) => i.severity === "error") && (
557+ <span className="pl-badge pl-badge-error" style={{ marginLeft: 6 }}>
558+ check
559+ </span>
560+ )}
561+ {props.disabled && (
562+ <span className="pl-badge pl-badge-info" style={{ marginLeft: 6 }}>
563+ not used by this method
564+ </span>
565+ )}
248566 </summary>
249567 {props.children}
250568 {related.map((i) => (
251- <div key={i.path + i.message} className={`pl-issue${i.severity === "warning" ? " pl-issue-warning" : ""}`}>
569+ <div
570+ key={i.path + i.message}
571+ className={`pl-issue${i.severity === "warning" ? " pl-issue-warning" : ""}`}
572+ >
252573 {i.path}: {i.message}
253574 </div>
254575 ))}
packages/core/src/generate/library.tsmodified+82−14View file
@@ -15,6 +15,8 @@ import plate_compare from "../m/plate_compare.m?raw";
1515 import plate_implicit from "../m/plate_implicit.m?raw";
1616 import plate_implicit_plan from "../m/plate_implicit_plan.m?raw";
1717 import plate_matrices from "../m/plate_matrices.m?raw";
18+import plate_rigid_modes from "../m/plate_rigid_modes.m?raw";
19+import plate_curvature_free from "../m/plate_curvature_free.m?raw";
1820 import plate_result_schema from "../m/plate_result_schema.m?raw";
1921
2022 export interface SourceFile {
@@ -29,24 +31,90 @@ export interface LibraryFile extends SourceFile {
2931
3032 /** Every `.m` file the solvers can call. The runner ships all of them with each run. */
3133 export const LIBRARY_FILES: LibraryFile[] = [
32- { name: "plate_modal.m", source: plate_modal, summary: "Modal superposition on a simply supported rectangle" },
34+ {
35+ name: "plate_modal.m",
36+ source: plate_modal,
37+ summary: "Modal superposition on a simply supported rectangle",
38+ },
3339 { name: "plate_fd.m", source: plate_fd, summary: "Explicit finite-difference time stepping" },
34- { name: "plate_fd_plan.m", source: plate_fd_plan, summary: "Grid, stable time step and cost of the FD scheme" },
35- { name: "plate_implicit.m", source: plate_implicit, summary: "Implicit theta-scheme: sparse Cholesky once, two triangular solves per step" },
36- { name: "plate_implicit_plan.m", source: plate_implicit_plan, summary: "Grid, accuracy-chosen time step and size of the implicit scheme" },
37- { name: "plate_matrices.m", source: plate_matrices, summary: "Sparse Laplacian and biharmonic matrices on the interior nodes" },
38- { name: "plate_setup.m", source: plate_setup, summary: "Bending stiffness and other derived constants" },
39- { name: "plate_force_time.m", source: plate_force_time, summary: "Force history F(t): hammer pulse or harmonic drive" },
40- { name: "plate_force_footprint.m", source: plate_force_footprint, summary: "Raised-cosine spatial footprint of the load" },
41- { name: "plate_stress.m", source: plate_stress, summary: "Surface stresses, von Mises and principal stress from curvatures" },
42- { name: "plate_sdof.m", source: plate_sdof, summary: "Exact step response of one damped oscillator" },
43- { name: "plate_listener.m", source: plate_listener, summary: "Rayleigh-integral weights and delays to the listener" },
40+ {
41+ name: "plate_fd_plan.m",
42+ source: plate_fd_plan,
43+ summary: "Grid, stable time step and cost of the FD scheme",
44+ },
45+ {
46+ name: "plate_implicit.m",
47+ source: plate_implicit,
48+ summary: "Implicit theta-scheme: sparse Cholesky once, two triangular solves per step",
49+ },
50+ {
51+ name: "plate_implicit_plan.m",
52+ source: plate_implicit_plan,
53+ summary: "Grid, accuracy-chosen time step and size of the implicit scheme",
54+ },
55+ {
56+ name: "plate_matrices.m",
57+ source: plate_matrices,
58+ summary: "Sparse mass, Laplacian and stiffness matrices for every edge condition",
59+ },
60+ {
61+ name: "plate_rigid_modes.m",
62+ source: plate_rigid_modes,
63+ summary: "Rigid-body motions of a free plate and the projector that removes them",
64+ },
65+ {
66+ name: "plate_setup.m",
67+ source: plate_setup,
68+ summary: "Bending stiffness and other derived constants",
69+ },
70+ {
71+ name: "plate_force_time.m",
72+ source: plate_force_time,
73+ summary: "Force history F(t): hammer pulse or harmonic drive",
74+ },
75+ {
76+ name: "plate_force_footprint.m",
77+ source: plate_force_footprint,
78+ summary: "Raised-cosine spatial footprint of the load",
79+ },
80+ {
81+ name: "plate_stress.m",
82+ source: plate_stress,
83+ summary: "Surface stresses, von Mises and principal stress from curvatures",
84+ },
85+ {
86+ name: "plate_sdof.m",
87+ source: plate_sdof,
88+ summary: "Exact step response of one damped oscillator",
89+ },
90+ {
91+ name: "plate_listener.m",
92+ source: plate_listener,
93+ summary: "Rayleigh-integral weights and delays to the listener",
94+ },
4495 { name: "plate_spectrum.m", source: plate_spectrum, summary: "Windowed amplitude spectrum" },
45- { name: "plate_pad.m", source: plate_pad, summary: "Ghost layer implementing the edge conditions" },
96+ {
97+ name: "plate_pad.m",
98+ source: plate_pad,
99+ summary: "Ghost layer implementing the edge conditions",
100+ },
46101 { name: "plate_laplacian.m", source: plate_laplacian, summary: "Five-point Laplacian" },
47- { name: "plate_curvature.m", source: plate_curvature, summary: "Second derivatives by central differences" },
102+ {
103+ name: "plate_curvature.m",
104+ source: plate_curvature,
105+ summary: "Second derivatives by central differences",
106+ },
107+ {
108+ name: "plate_curvature_free.m",
109+ source: plate_curvature_free,
110+ summary: "Curvatures of a free plate, edges from the moment condition",
111+ },
48112 { name: "plate_compare.m", source: plate_compare, summary: "Error norms between two solutions" },
49- { name: "plate_result_schema.m", source: plate_result_schema, summary: "Documentation of the result struct" },
113+ {
114+ name: "plate_result_schema.m",
115+ source: plate_result_schema,
116+ summary: "Documentation of the result struct",
117+ },
50118 ];
51119
52120 export function libraryFile(name: string): LibraryFile {
packages/core/src/generate/params.tsmodified+1−1View file
@@ -17,7 +17,7 @@ export function generateParamsBlock(spec: ProblemSpec): string {
1717 "params = struct();",
1818 `params.plate = struct('Lx', ${n(spec.plate.Lx)}, 'Ly', ${n(spec.plate.Ly)}, 'h', ${n(spec.plate.h)}); % side lengths and thickness [m]`,
1919 `params.material = struct('name', ${s(spec.material.name)}, 'E', ${n(spec.material.E)}, 'nu', ${n(spec.material.nu)}, 'rho', ${n(spec.material.rho)}); % Young's modulus [Pa], Poisson's ratio, density [kg/m^3]`,
20- `params.bc = ${s(spec.bc)}; % 'simply-supported' or 'clamped' (clamped: finite differences only)`,
20+ `params.bc = ${s(spec.bc)}; % 'simply-supported', 'clamped' or 'free' (modal: simply-supported only)`,
2121 `params.damping = struct('sigma0', ${n(spec.damping.sigma0)}, 'sigma1', ${n(spec.damping.sigma1)}); % loss rates: frequency-independent [1/s], frequency-dependent [m^2/s]`,
2222 `params.excitation = ${excitation}`,
2323 `params.probe = struct('x', ${n(spec.probe.x)}, 'y', ${n(spec.probe.y)}); % where the time series are recorded [m]`,
packages/core/src/m/plate_curvature_free.madded+33−0View file
@@ -0,0 +1,33 @@
1+function [wxx, wyy, wxy] = plate_curvature_free(w, dx, dy, nu)
2+% PLATE_CURVATURE_FREE Second derivatives of the displacement of a plate with free edges.
3+%
4+% [wxx, wyy, wxy] = plate_curvature_free(w, dx, dy, nu) returns the curvatures
5+% at every node of the ny x nx field w (rows index y, columns index x).
6+% Inside they are the usual central differences. On the edges there is no
7+% ghost layer to lean on, so the free-edge conditions themselves supply the
8+% missing values:
9+%
10+% no bending moment on an edge w_nn = -nu w_tt (n normal, t tangent)
11+% no twisting moment at a corner w_xy = 0, hence also w_xx = w_yy = 0 there
12+%
13+% The tangential curvature along an edge is a central difference along the
14+% edge, and the twist on an edge is one-sided in the normal direction. Both
15+% direct stresses therefore vanish at a corner, and on an edge the normal
16+% stress vanishes while the tangential stress is -E z w_tt.
17+
18+[ny, nx] = size(w);
19+wxx = zeros(ny, nx);
20+wyy = zeros(ny, nx);
21+wxy = zeros(ny, nx);
22+wxx(:, 2:end - 1) = (w(:, 3:end) - 2 * w(:, 2:end - 1) + w(:, 1:end - 2)) / dx^2;
23+wyy(2:end - 1, :) = (w(3:end, :) - 2 * w(2:end - 1, :) + w(1:end - 2, :)) / dy^2;
24+% Edge columns and rows: the moment condition, using the tangential curvature just computed
25+wxx(2:end - 1, [1 end]) = -nu * wyy(2:end - 1, [1 end]);
26+wyy([1 end], 2:end - 1) = -nu * wxx([1 end], 2:end - 1);
27+% Twist: central inside, one-sided in the normal direction on the edges, zero at the corners
28+wxy(2:end - 1, 2:end - 1) = (w(3:end, 3:end) - w(3:end, 1:end - 2) - w(1:end - 2, 3:end) + w(1:end - 2, 1:end - 2)) / (4 * dx * dy);
29+wxy(2:end - 1, 1) = (w(3:end, 2) - w(3:end, 1) - w(1:end - 2, 2) + w(1:end - 2, 1)) / (2 * dx * dy);
30+wxy(2:end - 1, end) = (w(3:end, end) - w(3:end, end - 1) - w(1:end - 2, end) + w(1:end - 2, end - 1)) / (2 * dx * dy);
31+wxy(1, 2:end - 1) = (w(2, 3:end) - w(2, 1:end - 2) - w(1, 3:end) + w(1, 1:end - 2)) / (2 * dx * dy);
32+wxy(end, 2:end - 1) = (w(end, 3:end) - w(end, 1:end - 2) - w(end - 1, 3:end) + w(end - 1, 1:end - 2)) / (2 * dx * dy);
33+end
packages/core/src/m/plate_fd.mmodified+33−8View file
@@ -13,8 +13,11 @@ function result = plate_fd(params)
1313 % where w-, w, w+ are consecutive time levels and L is the five-point
1414 % Laplacian. Applying L twice gives the 13-point biharmonic stencil. The
1515 % scheme is second order in space and time (Bilbao, Numerical Sound
16-% Synthesis, 2009, ch. 12). Edge conditions enter through one layer of ghost
17-% nodes (plate_pad): simply supported or clamped.
16+% Synthesis, 2009, ch. 12). Simply supported and clamped edges enter through
17+% one layer of ghost nodes (plate_pad). Free edges have no ghost-layer form,
18+% so for them L and L(L .) are the sparse energy-method operators of
19+% plate_matrices applied to every node, and the three rigid-body motions are
20+% projected out of the displacement after each step (plate_rigid_modes).
1821 %
1922 % The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
2023 % is unstable for cfl > 1 and the run stops as soon as the field diverges,
@@ -65,6 +68,15 @@ wgt = L.weight(:);
6568 na = ceil(nsteps * dt * fs) + L.nbins + 1;
6669 Pbins = zeros(L.nbins, na);
6770
71+% Free edges: sparse operators on every node, and the rigid-body projector
72+free = strcmp(P.bc, 'free');
73+if free
74+ [Lm, Km, Mm] = plate_matrices(nx, ny, dx, dy, P.bc, P.material.nu);
75+ mvec = full(diag(Mm));
76+ [Rb, Qb] = plate_rigid_modes(X, Y, mvec);
77+ notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
78+end
79+
6880 % State
6981 w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
7082 w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
@@ -78,12 +90,21 @@ diverged_step = 0;
7890 kf = 1;
7991
8092 for nstep = 1:nsteps
81- Lw = plate_laplacian(plate_pad(w, P.bc), dx, dy);
82- B = zeros(ny, nx);
83- B(2:end - 1, 2:end - 1) = (Lw(3:end, 2:end - 1) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(1:end - 2, 2:end - 1)) / dy^2 ...
84- + (Lw(2:end - 1, 3:end) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(2:end - 1, 1:end - 2)) / dx^2;
93+ if free
94+ Lw = reshape((Lm * w(:)) ./ mvec, ny, nx);
95+ B = reshape((Km * w(:)) ./ mvec, ny, nx);
96+ else
97+ Lw = plate_laplacian(plate_pad(w, P.bc), dx, dy);
98+ B = zeros(ny, nx);
99+ B(2:end - 1, 2:end - 1) = (Lw(3:end, 2:end - 1) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(1:end - 2, 2:end - 1)) / dy^2 ...
100+ + (Lw(2:end - 1, 3:end) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(2:end - 1, 1:end - 2)) / dx^2;
101+ end
85102 wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
86- wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;
103+ if free
104+ wn = wn - reshape(Rb * (Qb' * wn(:)), ny, nx);
105+ else
106+ wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;
107+ end
87108
88109 % Probe displacement and curvatures from the 3 x 3 neighbourhood
89110 w_p(nstep + 1) = wn(ip, jp);
@@ -97,7 +118,11 @@ for nstep = 1:nsteps
97118 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [L.nbins 1]) * (dt * fs);
98119
99120 if kf <= nf && nstep == frame_steps(kf)
100- [wxx, wyy, wxy] = plate_curvature(plate_pad(wn, P.bc), dx, dy);
121+ if free
122+ [wxx, wyy, wxy] = plate_curvature_free(wn, dx, dy, P.material.nu);
123+ else
124+ [wxx, wyy, wxy] = plate_curvature(plate_pad(wn, P.bc), dx, dy);
125+ end
101126 S = plate_stress(P, wxx, wyy, wxy);
102127 Wf(:, :, kf) = wn; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
103128 kf = kf + 1;
packages/core/src/m/plate_fd_plan.mmodified+5−0View file
@@ -18,6 +18,11 @@ function plan = plate_fd_plan(params)
1818 % which for a square grid without sigma1 is the familiar dt <= dx^2 / (4 kappa).
1919 % Because dt scales with dx^2, halving the grid spacing costs eight times
2020 % as many cell updates for the same simulated time.
21+%
22+% The same bound serves the free plate: the largest eigenvalue of its
23+% energy-method operator M \ K (plate_matrices) lies just below lambda_max^2
24+% (about 0.97 to 0.99 of it on the grids used here), so the step is safe and
25+% slightly conservative there.
2126
2227 P = plate_setup(params);
2328 plan.nx = P.fd.nx;
packages/core/src/m/plate_implicit.mmodified+42−20View file
@@ -11,14 +11,19 @@ function result = plate_implicit(params)
1111 % (w+ - 2w + w-)/dt^2 + sigma0 (w+ - w-)/dt - sigma1 L (w+ - w-)/dt
1212 % + kappa^2 B (theta w+ + (1 - 2 theta) w + theta w-) = F s / mu,
1313 %
14-% with L and B the sparse Laplacian and biharmonic matrices on the interior
15-% nodes (plate_matrices). Multiplying by dt^2 and collecting terms gives one
16-% linear system per step with a constant matrix,
14+% with L and B the sparse Laplacian and biharmonic matrices of plate_matrices,
15+% which also supplies the diagonal mass matrix M: the identity on the interior
16+% nodes for simply supported or clamped edges, the trapezoidal node weights on
17+% every node for free edges. Multiplying by dt^2 and collecting terms gives
18+% one linear system per step with a constant matrix,
1719 %
18-% A w+ = C w - A_old w- + dt^2 F s / mu
19-% A = (1 + sigma0 dt) I - sigma1 dt L + theta kappa^2 dt^2 B
20-% C = 2 I - (1 - 2 theta) kappa^2 dt^2 B
21-% A_old = (1 - sigma0 dt) I + sigma1 dt L + theta kappa^2 dt^2 B.
20+% A w+ = C w - A_old w- + dt^2 F M s / mu
21+% A = (1 + sigma0 dt) M - sigma1 dt L + theta kappa^2 dt^2 B
22+% C = 2 M - (1 - 2 theta) kappa^2 dt^2 B
23+% A_old = (1 - sigma0 dt) M + sigma1 dt L + theta kappa^2 dt^2 B.
24+%
25+% For free edges the three rigid-body motions are projected out of w+ after
26+% every solve (plate_rigid_modes).
2227 %
2328 % A is symmetric positive definite, so it is factored once by sparse Cholesky
2429 % and each step costs two sparse triangular solves plus two sparse products.
@@ -46,14 +51,20 @@ if theta < 0.25
4651 notes{end + 1} = sprintf('theta = %g is below 1/4, so the scheme is only conditionally stable', theta);
4752 end
4853
49-% ---- Matrices on the interior nodes, factored once
50-[L, B, mx, my] = plate_matrices(nx, ny, dx, dy, P.bc);
51-m = mx * my;
52-I = speye(m);
54+% ---- Matrices on the unknowns, factored once
55+free = strcmp(P.bc, 'free');
56+[L, B, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, P.bc, P.material.nu);
57+m = numel(idx);
58+mvec = full(diag(M));
59+if free
60+ rr = 1:ny; cc = 1:nx; % every node is an unknown
61+else
62+ rr = 2:ny - 1; cc = 2:nx - 1; % the edge nodes carry w = 0
63+end
5364 k2dt2 = P.kappa^2 * dt^2;
54-A = (1 + P.damping.sigma0 * dt) * I - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
55-C = 2 * I - (1 - 2 * theta) * k2dt2 * B;
56-A_old = (1 - P.damping.sigma0 * dt) * I + P.damping.sigma1 * dt * L + theta * k2dt2 * B;
65+A = (1 + P.damping.sigma0 * dt) * M - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
66+C = 2 * M - (1 - 2 * theta) * k2dt2 * B;
67+A_old = (1 - P.damping.sigma0 * dt) * M + P.damping.sigma1 * dt * L + theta * k2dt2 * B;
5768 t_factor = tic;
5869 R = chol(A); % A = R' * R
5970 Rt = R'; % transposed once, not once per step
@@ -61,8 +72,12 @@ factor_s = toc(t_factor);
6172
6273 % ---- Load, probe, frames, listener (as in plate_fd)
6374 s = plate_force_footprint(P.excitation, X, Y, dx * dy);
64-s_int = s(2:end - 1, 2:end - 1);
65-load = (dt^2 / P.mu) * s_int(:);
75+s_unk = s(rr, cc);
76+load = (dt^2 / P.mu) * (mvec .* s_unk(:)); % M s on the unknowns
77+if free
78+ [Rb, Qb] = plate_rigid_modes(X(rr, cc), Y(rr, cc), mvec);
79+ notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
80+end
6681 tt = (0:nsteps) * dt;
6782 F = plate_force_time(P.excitation, tt);
6883
@@ -83,7 +98,7 @@ wgt = Lis.weight(:);
8398 na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
8499 Pbins = zeros(Lis.nbins, na);
85100
86-% ---- Time stepping on the interior vector
101+% ---- Time stepping on the vector of unknowns
87102 w = zeros(m, 1); wm = zeros(m, 1);
88103 wfull = zeros(ny, nx);
89104 w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
@@ -94,20 +109,27 @@ kf = 1;
94109 for nstep = 1:nsteps
95110 rhs = C * w - A_old * wm + F(nstep) * load;
96111 wn = R \ (Rt \ rhs); % forward then back substitution
112+ if free
113+ wn = wn - Rb * (Qb' * wn);
114+ end
97115
98- wfull(2:end - 1, 2:end - 1) = reshape(wn, my, mx);
116+ wfull(rr, cc) = reshape(wn, my, mx);
99117 w_p(nstep + 1) = wfull(ip, jp);
100118 wxx_p(nstep + 1) = (wfull(ip, jp + 1) - 2 * wfull(ip, jp) + wfull(ip, jp - 1)) / dx^2;
101119 wyy_p(nstep + 1) = (wfull(ip + 1, jp) - 2 * wfull(ip, jp) + wfull(ip - 1, jp)) / dy^2;
102120 wxy_p(nstep + 1) = (wfull(ip + 1, jp + 1) - wfull(ip + 1, jp - 1) - wfull(ip - 1, jp + 1) + wfull(ip - 1, jp - 1)) / (4 * dx * dy);
103121
104122 acc = zeros(ny, nx);
105- acc(2:end - 1, 2:end - 1) = reshape((wn - 2 * w + wm) / dt^2, my, mx);
123+ acc(rr, cc) = reshape((wn - 2 * w + wm) / dt^2, my, mx);
106124 ia = floor(nstep * dt * fs) + 1;
107125 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [Lis.nbins 1]) * (dt * fs);
108126
109127 if kf <= nf && nstep == frame_steps(kf)
110- [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
128+ if free
129+ [wxx, wyy, wxy] = plate_curvature_free(wfull, dx, dy, P.material.nu);
130+ else
131+ [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
132+ end
111133 S = plate_stress(P, wxx, wyy, wxy);
112134 Wf(:, :, kf) = wfull; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
113135 kf = kf + 1;
packages/core/src/m/plate_implicit_plan.mmodified+5−1View file
@@ -21,7 +21,11 @@ plan.dy = P.plate.Ly / (plan.ny - 1);
2121 plan.dt = 1 / (P.output.fs * P.implicit.steps_per_sample);
2222 plan.nsteps = ceil(P.implicit.duration / plan.dt);
2323 plan.ncells = plan.nx * plan.ny;
24-plan.nunknowns = (plan.nx - 2) * (plan.ny - 2);
24+if strcmp(P.bc, 'free')
25+ plan.nunknowns = plan.nx * plan.ny; % every node moves; see plate_matrices
26+else
27+ plan.nunknowns = (plan.nx - 2) * (plan.ny - 2);
28+end
2529 plan.cell_steps = plan.ncells * plan.nsteps;
2630 plan.theta = P.implicit.theta;
2731 end
packages/core/src/m/plate_matrices.mmodified+93−34View file
@@ -1,43 +1,102 @@
1-function [L, B, mx, my] = plate_matrices(nx, ny, dx, dy, bc)
2-% PLATE_MATRICES Sparse Laplacian and biharmonic matrices on the interior nodes.
1+function [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu)
2+% PLATE_MATRICES Sparse Laplacian, stiffness and mass matrices for the plate.
33 %
4-% [L, B, mx, my] = plate_matrices(nx, ny, dx, dy, bc) assembles, for the
5-% mx = nx-2 by my = ny-2 interior nodes of the grid (edge nodes carry w = 0),
6-% L the five-point Laplacian (5 nonzeros per row)
7-% B the thirteen-point biharmonic (13 nonzeros per row)
8-% as sparse matrices acting on w(:) with the interior field stored as an
9-% my x mx array, so the y index runs fastest, matching w(2:end-1, 2:end-1)(:).
4+% [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu) assembles
5+% the spatial operators of the finite-difference schemes as sparse matrices
6+% acting on the vector w(idx) of unknowns, where idx indexes the ny x nx field
7+% with y running fastest (column-major, as w(:) does) and the unknowns form an
8+% my x mx array. The semi-discrete equation is
109 %
11-% The stencils are the same as the ghost-node scheme in plate_fd: the second
12-% difference with w = 0 at the edge is D2; the fourth difference is D2*D2 for
13-% simply supported edges (odd reflection, w_nn = 0), and for clamped edges the
14-% even reflection w_ghost = w_1 changes the end coefficient of the fourth
15-% difference from 5/dx^4 to 7/dx^4. The mixed term 2 w_xxyy never reaches a
16-% ghost node. A test in the lab checks B against the explicit stencil.
10+% M w_tt + 2 sigma0 M w_t - 2 sigma1 L w_t + kappa^2 K w = F(t) M s / mu
1711 %
18-% Everything is built from Kronecker products so the matrices stay sparse.
19-
20-mx = nx - 2;
21-my = ny - 2;
22-ex = ones(mx, 1);
23-ey = ones(my, 1);
24-D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
25-D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
26-Ix = speye(mx);
27-Iy = speye(my);
12+% Simply supported and clamped edges (w = 0 on the edge nodes). The unknowns
13+% are the mx = nx-2 by my = ny-2 interior nodes, M is the identity, L the
14+% five-point Laplacian and K the thirteen-point biharmonic, with the edges
15+% handled exactly as in the ghost-node scheme of plate_fd / plate_pad: D2*D2
16+% encodes the odd reflection (w_nn = 0), and for clamped edges the even
17+% reflection changes the end coefficient of the fourth difference from 5/dx^4
18+% to 7/dx^4. The mixed term 2 w_xxyy never reaches a ghost node. nu is not
19+% used. A test in the lab checks K against the explicit stencil.
20+%
21+% Free edges. Every node is an unknown (mx = nx, my = ny). The free-edge
22+% conditions, no bending moment and no Kirchhoff shear on an edge and no
23+% twisting moment at a corner, involve third derivatives and do not fit in a
24+% ghost layer. K is therefore assembled from the bending energy
25+%
26+% U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
27+%
28+% using second differences at the nodes, the cross difference at the cell
29+% centres and trapezoidal quadrature weights (the finite-difference energy
30+% method; Bushnell, Computers & Structures 1973). The free-edge conditions
31+% are then natural boundary conditions of the discrete energy: the scheme is
32+% second order inside and first order at the boundary. M holds the trapezoidal
33+% node weights (1 inside, 1/2 on an edge, 1/4 at a corner) and L is the
34+% Laplacian energy operator with the natural zero-normal-derivative condition.
35+% K is symmetric positive semidefinite and its null space is exactly the three
36+% rigid-body motions 1, x and y, which plate_rigid_modes removes from the
37+% state. All three matrices are divided by dx*dy, so inside the plate M = 1
38+% and K and L reduce to the usual stencils, and the largest eigenvalue of
39+% M \ K stays below the interior bound (4/dx^2 + 4/dy^2)^2 that plate_fd_plan
40+% uses for the time step.
41+%
42+% Everything is built from Kronecker products so the matrices stay sparse;
43+% kron(Ax, By) acts on w(:) as By * W * Ax'.
2844
29-D4x = D2x * D2x;
30-D4y = D2y * D2y;
45+if nargin < 6
46+ nu = NaN;
47+end
3148 switch bc
32- case 'simply-supported'
33- % D2*D2 already encodes the odd reflection
34- case 'clamped'
35- D4x(1, 1) = 7 / dx^4; D4x(mx, mx) = 7 / dx^4;
36- D4y(1, 1) = 7 / dy^4; D4y(my, my) = 7 / dy^4;
49+ case {'simply-supported', 'clamped'}
50+ mx = nx - 2;
51+ my = ny - 2;
52+ ex = ones(mx, 1);
53+ ey = ones(my, 1);
54+ D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
55+ D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
56+ Ix = speye(mx);
57+ Iy = speye(my);
58+ D4x = D2x * D2x;
59+ D4y = D2y * D2y;
60+ if strcmp(bc, 'clamped')
61+ D4x(1, 1) = 7 / dx^4; D4x(mx, mx) = 7 / dx^4;
62+ D4y(1, 1) = 7 / dy^4; D4y(my, my) = 7 / dy^4;
63+ end
64+ L = kron(Ix, D2y) + kron(D2x, Iy);
65+ K = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
66+ M = speye(mx * my);
67+ mask = false(ny, nx);
68+ mask(2:end - 1, 2:end - 1) = true;
69+ idx = find(mask);
70+ case 'free'
71+ if ~isfinite(nu)
72+ error('plate_matrices: free edges need Poisson''s ratio nu');
73+ end
74+ mx = nx;
75+ my = ny;
76+ % One-dimensional pieces, as triplets so the layout is unambiguous:
77+ % second differences at the interior nodes ((n-2) x n), first differences
78+ % at the cell midpoints ((n-1) x n), selectors of the interior nodes, and
79+ % trapezoidal node weights.
80+ kx = (1:nx - 2)'; ky = (1:ny - 2)';
81+ d2x = sparse([kx; kx; kx], [kx; kx + 1; kx + 2], [ones(nx - 2, 1); -2 * ones(nx - 2, 1); ones(nx - 2, 1)], nx - 2, nx) / dx^2;
82+ d2y = sparse([ky; ky; ky], [ky; ky + 1; ky + 2], [ones(ny - 2, 1); -2 * ones(ny - 2, 1); ones(ny - 2, 1)], ny - 2, ny) / dy^2;
83+ jx = (1:nx - 1)'; jy = (1:ny - 1)';
84+ d1x = sparse([jx; jx], [jx; jx + 1], [-ones(nx - 1, 1); ones(nx - 1, 1)], nx - 1, nx) / dx;
85+ d1y = sparse([jy; jy], [jy; jy + 1], [-ones(ny - 1, 1); ones(ny - 1, 1)], ny - 1, ny) / dy;
86+ Sx = sparse(kx, kx + 1, ones(nx - 2, 1), nx - 2, nx);
87+ Sy = sparse(ky, ky + 1, ones(ny - 2, 1), ny - 2, ny);
88+ tx = ones(nx, 1); tx([1 end]) = 0.5; Tx = spdiags(tx, 0, nx, nx);
89+ ty = ones(ny, 1); ty([1 end]) = 0.5; Ty = spdiags(ty, 0, ny, ny);
90+ % The four terms of the energy density, each as D' * (weights) * D
91+ Kxx = kron(d2x' * d2x, Ty); % w_xx^2 at nodes 2..nx-1, every row
92+ Kyy = kron(Tx, d2y' * d2y); % w_yy^2 at nodes 2..ny-1, every column
93+ Kc = nu * (kron(d2x' * Sx, Sy' * d2y) + kron(Sx' * d2x, d2y' * Sy)); % 2 nu w_xx w_yy where both exist
94+ Kxy = 2 * (1 - nu) * kron(d1x' * d1x, d1y' * d1y); % 2 (1-nu) w_xy^2 at the cell centres
95+ K = Kxx + Kyy + Kc + Kxy;
96+ L = -(kron(d1x' * d1x, Ty) + kron(Tx, d1y' * d1y));
97+ M = kron(Tx, Ty);
98+ idx = (1:nx * ny)';
3799 otherwise
38100 error('plate_matrices: unknown boundary condition "%s"', bc);
39101 end
40-
41-L = kron(Ix, D2y) + kron(D2x, Iy);
42-B = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
43102 end
packages/core/src/m/plate_rigid_modes.madded+34−0View file
@@ -0,0 +1,34 @@
1+function [R, Q] = plate_rigid_modes(X, Y, mvec)
2+% PLATE_RIGID_MODES Basis of the rigid-body motions of a free plate and its M-orthogonal projector.
3+%
4+% [R, Q] = plate_rigid_modes(X, Y, mvec) takes the node coordinates X, Y and
5+% the diagonal mvec of the mass matrix (the trapezoidal weights of
6+% plate_matrices) and returns R = [1, x - xc, y - yc] evaluated at the nodes
7+% (n x 3) and Q = M R (R' M R)^-1 (n x 3), so that
8+%
9+% w - R * (Q' * w)
10+%
11+% is the part of w orthogonal to the rigid-body motions in the mass inner
12+% product: what remains when the mean displacement and the two mean tilts are
13+% taken out.
14+%
15+% Why remove them. A completely free plate has three zero-frequency modes,
16+% and a strike gives it net momentum and angular momentum: with only the loss
17+% sigma0 to stop it, the plate drifts to a displacement of order
18+% (impulse / mass) / (2 sigma0), millimetres against the micrometres of the
19+% bending motion. A real free-free measurement rests the plate on foam or
20+% hangs it from threads, which absorb exactly this motion and nothing else,
21+% so the solvers project it out after every step. The bending stiffness K
22+% annihilates R exactly, so this changes nothing about the flexural dynamics;
23+% the sigma1 loss operator does not, and the projection also keeps that term
24+% from feeding a rigid tilt back into flexure.
25+
26+x = X(:);
27+y = Y(:);
28+m = mvec(:);
29+xc = sum(m .* x) / sum(m);
30+yc = sum(m .* y) / sum(m);
31+R = [ones(numel(x), 1), x - xc, y - yc];
32+G = R' * (m .* R); % 3 x 3 Gram matrix in the mass inner product
33+Q = (m .* R) / G;
34+end
packages/core/src/m/plate_setup.mmodified+2−2View file
@@ -27,8 +27,8 @@ end
2727 if pl.h > 0.1 * min(pl.Lx, pl.Ly)
2828 warning('plate_setup: thickness is more than a tenth of the shorter side; Kirchhoff-Love thin-plate theory is not reliable here');
2929 end
30-if ~any(strcmp(params.bc, {'simply-supported', 'clamped'}))
31- error('plate_setup: unknown boundary condition "%s" (expected simply-supported or clamped)', params.bc);
30+if ~any(strcmp(params.bc, {'simply-supported', 'clamped', 'free'}))
31+ error('plate_setup: unknown boundary condition "%s" (expected simply-supported, clamped or free)', params.bc);
3232 end
3333
3434 P.D = mat.E * pl.h^3 / (12 * (1 - mat.nu^2));
packages/core/src/methods/common.tsmodified+14−0View file
@@ -19,6 +19,20 @@ export const GOVERNING_EQUATIONS: Equation[] = [
1919 },
2020 ];
2121
22+/** How the two finite-difference methods handle free edges, which have no ghost-layer form. */
23+export const FREE_EDGE_EQUATIONS: Equation[] = [
24+ {
25+ tex: String.raw`M_n = -D\left(w_{nn} + \nu\, w_{tt}\right) = 0, \qquad V_n = -D\left(w_{nnn} + (2-\nu)\, w_{ntt}\right) = 0, \qquad w_{nt} = 0 \text{ at a corner}`,
26+ caption:
27+ "Free edge: no bending moment and no Kirchhoff shear, n normal and t tangent to the edge; at a free corner the twisting moment vanishes as well. The shear condition holds a third derivative, so no single ghost layer can encode it.",
28+ },
29+ {
30+ tex: String.raw`U = \frac{D}{2}\iint \left[ w_{xx}^2 + w_{yy}^2 + 2\nu\, w_{xx} w_{yy} + 2(1-\nu)\, w_{xy}^2 \right] dA, \qquad \mathbf K = \sum_{q} \omega_q\, \mathsf D_q^{\!\top} \mathsf C\, \mathsf D_q`,
31+ caption:
32+ "For free edges the stiffness matrix K is assembled from the bending energy instead: finite differences of the curvatures at the nodes and cell centres with trapezoidal weights (the finite-difference energy method). K reproduces the 13-point stencil inside, satisfies the free-edge conditions as natural boundary conditions, and has exactly the three rigid-body motions in its null space; those are projected out after every step, as the soft supports of a real free-free rig would absorb them.",
33+ },
34+];
35+
2236 export const STRESS_EQUATIONS: Equation[] = [
2337 {
2438 tex: String.raw`\sigma_{xx} = -\frac{E z}{1-\nu^2}\left(w_{xx} + \nu w_{yy}\right),\quad \sigma_{yy} = -\frac{E z}{1-\nu^2}\left(w_{yy} + \nu w_{xx}\right),\quad \tau_{xy} = -\frac{E z}{1+\nu}\,w_{xy}`,
packages/core/src/methods/fd.tsmodified+13−3View file
@@ -1,5 +1,10 @@
11 import type { MethodDefinition } from "./types";
2-import { GOVERNING_EQUATIONS, RADIATION_EQUATIONS, STRESS_EQUATIONS } from "./common";
2+import {
3+ FREE_EDGE_EQUATIONS,
4+ GOVERNING_EQUATIONS,
5+ RADIATION_EQUATIONS,
6+ STRESS_EQUATIONS,
7+} from "./common";
38
49 export const FD_METHOD: MethodDefinition = {
510 id: "fd",
@@ -9,7 +14,7 @@ export const FD_METHOD: MethodDefinition = {
914 solverFunction: "plate_fd",
1015 planFunction: "plate_fd_plan",
1116 supports: {
12- bc: ["simply-supported", "clamped"],
17+ bc: ["simply-supported", "clamped", "free"],
1318 excitation: ["impulse", "harmonic"],
1419 audio: "clip",
1520 },
@@ -20,6 +25,9 @@ export const FD_METHOD: MethodDefinition = {
2025 "plate_pad.m",
2126 "plate_laplacian.m",
2227 "plate_curvature.m",
28+ "plate_matrices.m",
29+ "plate_rigid_modes.m",
30+ "plate_curvature_free.m",
2331 "plate_force_time.m",
2432 "plate_force_footprint.m",
2533 "plate_stress.m",
@@ -37,13 +45,14 @@ export const FD_METHOD: MethodDefinition = {
3745 {
3846 tex: String.raw`\mathsf L w_{ij} = \frac{w_{i,j+1} - 2w_{ij} + w_{i,j-1}}{\Delta x^2} + \frac{w_{i+1,j} - 2w_{ij} + w_{i-1,j}}{\Delta y^2}`,
3947 caption:
40- "Edge conditions enter through one ghost layer: odd reflection for simply supported (w = w_nn = 0), even reflection for clamped (w = w_n = 0).",
48+ "Supported and clamped edges enter through one ghost layer: odd reflection for simply supported (w = w_nn = 0), even reflection for clamped (w = w_n = 0). Free edges take the energy-method operators below.",
4149 },
4250 {
4351 tex: String.raw`\Delta t \le \frac{2}{\lambda_{\max}\left(\sigma_1 + \sqrt{\kappa^2 + \sigma_1^2}\right)}, \qquad \lambda_{\max} = \frac{4}{\Delta x^2} + \frac{4}{\Delta y^2}`,
4452 caption:
4553 "Stability bound from von Neumann analysis; for a square grid without sigma1 it reads dt <= dx^2 / (4 kappa). Because dt scales with dx^2, halving the spacing costs eight times more work.",
4654 },
55+ ...FREE_EDGE_EQUATIONS,
4756 ...STRESS_EQUATIONS,
4857 ...RADIATION_EQUATIONS,
4958 ],
@@ -53,6 +62,7 @@ export const FD_METHOD: MethodDefinition = {
5362 "Stress fields are evaluated only at the saved frames, so the envelope samples the true peak.",
5463 "The time-step limit makes long simulated times expensive: cost grows as (1/dx)^4 for a fixed duration.",
5564 "Audio is a short clip covering only the simulated window, box-averaged onto audio samples.",
65+ "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
5666 ],
5767 durationOf: (spec) => spec.fd.duration,
5868 gridOf: (spec) => [spec.fd.nx, spec.fd.ny],
packages/core/src/methods/implicit.tsmodified+11−2View file
@@ -1,5 +1,10 @@
11 import type { MethodDefinition } from "./types";
2-import { GOVERNING_EQUATIONS, RADIATION_EQUATIONS, STRESS_EQUATIONS } from "./common";
2+import {
3+ FREE_EDGE_EQUATIONS,
4+ GOVERNING_EQUATIONS,
5+ RADIATION_EQUATIONS,
6+ STRESS_EQUATIONS,
7+} from "./common";
38
49 export const IMPLICIT_METHOD: MethodDefinition = {
510 id: "implicit",
@@ -10,7 +15,7 @@ export const IMPLICIT_METHOD: MethodDefinition = {
1015 solverFunction: "plate_implicit",
1116 planFunction: "plate_implicit_plan",
1217 supports: {
13- bc: ["simply-supported", "clamped"],
18+ bc: ["simply-supported", "clamped", "free"],
1419 excitation: ["impulse", "harmonic"],
1520 audio: "clip",
1621 },
@@ -18,6 +23,8 @@ export const IMPLICIT_METHOD: MethodDefinition = {
1823 "plate_implicit.m",
1924 "plate_implicit_plan.m",
2025 "plate_matrices.m",
26+ "plate_rigid_modes.m",
27+ "plate_curvature_free.m",
2128 "plate_setup.m",
2229 "plate_force_time.m",
2330 "plate_force_footprint.m",
@@ -45,6 +52,7 @@ export const IMPLICIT_METHOD: MethodDefinition = {
4552 caption:
4653 "No stability limit, so the step is chosen for accuracy. The price is period elongation: modes near half the sample rate run slow by tens of percent unless several steps per sample are used.",
4754 },
55+ ...FREE_EDGE_EQUATIONS,
4856 ...STRESS_EQUATIONS,
4957 ...RADIATION_EQUATIONS,
5058 ],
@@ -54,6 +62,7 @@ export const IMPLICIT_METHOD: MethodDefinition = {
5462 "The probe is snapped to the nearest interior grid node; stress fields are evaluated only at the saved frames.",
5563 "Cost per step is a sparse triangular solve pair, about a millisecond for a 61 x 41 grid in the browser, independent of the plate stiffness.",
5664 "Audio covers only the simulated window.",
65+ "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
5766 ],
5867 durationOf: (spec) => spec.implicit.duration,
5968 gridOf: (spec) => [spec.implicit.nx, spec.implicit.ny],
packages/core/src/spec/types.tsmodified+2−1View file
@@ -7,7 +7,8 @@
77 * Units are SI throughout: metres, seconds, newtons, pascals, kilograms.
88 */
99
10-export type BoundaryCondition = "simply-supported" | "clamped";
10+/** The same condition on all four edges; "free" is the free-free plate of Chladni figures and modal testing. */
11+export type BoundaryCondition = "simply-supported" | "clamped" | "free";
1112
1213 export interface PlateGeometry {
1314 /** Side along x [m] */
packages/core/test/free.test.tsadded+168−0View file
@@ -0,0 +1,168 @@
1+import { describe, expect, it } from "vitest";
2+import { LIBRARY_FILES, buildCompareBundle, isPlainStruct, parseCompareResult } from "../src/index";
3+import { maxAbs, runner, solve, testSpec } from "./helpers";
4+
5+/** Run a MATLAB snippet against the library and return the struct it leaves in `out`. */
6+async function script(source: string) {
7+ const run = await runner.run({
8+ bundle: {
9+ files: [{ name: "c.m", source }, ...LIBRARY_FILES],
10+ entry: "c.m",
11+ resultVariable: "out",
12+ },
13+ });
14+ if (!run.ok) throw new Error(`${run.error}\n${run.output.join("")}`);
15+ if (!isPlainStruct(run.value)) throw new Error("unexpected shape");
16+ return run.value.fields as Record<string, number>;
17+}
18+
19+describe("free edges", () => {
20+ it("assembles a symmetric stiffness whose null space is exactly the rigid-body motions and whose interior is the 13-point stencil", async () => {
21+ const out = await script(`
22+ nx = 19; ny = 13; dx = 0.3/(nx-1); dy = 0.2/(ny-1); nu = 0.3;
23+ [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, 'free', nu);
24+ x = linspace(0, 0.3, nx); y = linspace(0, 0.2, ny); [X, Y] = meshgrid(x, y);
25+ out.symm = norm(full(K - K'), 1) / norm(full(K), 1);
26+ out.null1 = max(abs(K * ones(nx*ny, 1)));
27+ out.nullx = max(abs(K * X(:))) * dx^4;
28+ out.nully = max(abs(K * Y(:))) * dx^4;
29+ out.lnull = max(abs(L * ones(nx*ny, 1)));
30+ % a curved field: the interior rows must be the biharmonic stencil, per unit mass
31+ W = sin(3*X) .* cos(5*Y) + X.^2 .* Y;
32+ Kw = reshape(K * W(:), ny, nx); Mv = reshape(full(diag(M)), ny, nx);
33+ Lw = plate_laplacian(plate_pad(W, 'simply-supported'), dx, dy); % pad only feeds the outer ring, dropped below
34+ Bw = (Lw(3:end,2:end-1) - 2*Lw(2:end-1,2:end-1) + Lw(1:end-2,2:end-1))/dy^2 + (Lw(2:end-1,3:end) - 2*Lw(2:end-1,2:end-1) + Lw(2:end-1,1:end-2))/dx^2;
35+ inner = Kw(3:end-2, 3:end-2) ./ Mv(3:end-2, 3:end-2);
36+ ref = Bw(2:end-1, 2:end-1);
37+ out.interior = norm(inner(:) - ref(:)) / norm(ref(:));
38+ out.m_inner = Mv(5, 7); out.m_edge = Mv(1, 7); out.m_corner = Mv(1, 1);
39+ out.nunk = numel(idx); out.mx = mx; out.my = my;
40+ `);
41+ expect(out.symm).toBeLessThan(1e-12);
42+ expect(out.null1).toBeLessThan(1e-6);
43+ expect(out.nullx).toBeLessThan(1e-9);
44+ expect(out.nully).toBeLessThan(1e-9);
45+ expect(out.lnull).toBeLessThan(1e-9);
46+ expect(out.interior).toBeLessThan(1e-10);
47+ expect(out.m_inner).toBe(1);
48+ expect(out.m_edge).toBe(0.5);
49+ expect(out.m_corner).toBe(0.25);
50+ expect(out.nunk).toBe(19 * 13);
51+ expect(out.mx).toBe(19);
52+ expect(out.my).toBe(13);
53+ });
54+
55+ it("reproduces the frequencies of a completely free square plate (Leissa, nu = 0.3)", async () => {
56+ // omega a^2 sqrt(rho h / D) for the first three elastic modes of the free square plate,
57+ // Ritz values tabulated by Leissa (1969): 13.49, 19.79, 24.43. The energy method converges
58+ // to them from below at about first order in dx (the edge treatment), so a 17-node grid
59+ // should land within a couple of percent, and the three rigid modes must sit at zero.
60+ const out = await script(`
61+ n = 17; a = 1; h = a/(n-1); nu = 0.3;
62+ [~, K, M] = plate_matrices(n, n, h, h, 'free', nu);
63+ Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh;
64+ S = (S + S') / 2;
65+ ev = sort(eig(S));
66+ out.zero3 = max(abs(ev(1:3))) / ev(end);
67+ out.f1 = sqrt(ev(4)) * a^2; out.f2 = sqrt(ev(5)) * a^2; out.f3 = sqrt(ev(6)) * a^2;
68+ out.lam_ratio = ev(end) / (4/h^2 + 4/h^2)^2;
69+ `);
70+ expect(out.zero3).toBeLessThan(1e-9);
71+ expect(Math.abs(out.f1 - 13.49) / 13.49).toBeLessThan(0.02);
72+ // and it lies below the simply supported square's 2 pi^2 = 19.74: free edges are the softest support
73+ expect(out.f1).toBeLessThan(2 * Math.PI ** 2);
74+ expect(Math.abs(out.f2 - 19.79) / 19.79).toBeLessThan(0.03);
75+ expect(Math.abs(out.f3 - 24.43) / 24.43).toBeLessThan(0.03);
76+ // the explicit scheme's CFL bound, derived for the interior stencil, remains an upper bound
77+ expect(out.lam_ratio).toBeLessThan(1);
78+ expect(out.lam_ratio).toBeGreaterThan(0.9);
79+ });
80+
81+ it("runs the explicit scheme on a free plate: stable below the CFL limit, no rigid-body drift, corners stress-free", async () => {
82+ const spec = testSpec();
83+ spec.bc = "free";
84+ spec.fd.cfl = 0.98;
85+ const { result } = await solve(spec, "fd");
86+ expect(result.info.stable).toBe(true);
87+ expect(result.probe.w.every(Number.isFinite)).toBe(true);
88+ expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
89+ expect(result.info.notes.some((n) => /rigid-body/.test(n))).toBe(true);
90+ // the last frame carries no mean displacement and no mean tilt in the mass inner product
91+ // (trapezoidal node weights: 1 inside, 1/2 on an edge, 1/4 at a corner), which is the
92+ // inner product the projector uses
93+ const { nx, ny, nf, data } = result.frames.w;
94+ const k = nf - 1;
95+ let sum = 0,
96+ sx = 0,
97+ sy = 0,
98+ peak = 0,
99+ mass = 0;
100+ for (let j = 0; j < nx; j++)
101+ for (let i = 0; i < ny; i++) {
102+ const v = data[i + j * ny + k * nx * ny];
103+ const m = (j === 0 || j === nx - 1 ? 0.5 : 1) * (i === 0 || i === ny - 1 ? 0.5 : 1);
104+ mass += m;
105+ sum += m * v;
106+ sx += m * v * (j - (nx - 1) / 2);
107+ sy += m * v * (i - (ny - 1) / 2);
108+ peak = Math.max(peak, Math.abs(v));
109+ }
110+ expect(peak).toBeGreaterThan(0);
111+ expect(Math.abs(sum) / mass).toBeLessThan(1e-9 * peak);
112+ expect(Math.abs(sx) / (mass * nx)).toBeLessThan(1e-9 * peak);
113+ expect(Math.abs(sy) / (mass * ny)).toBeLessThan(1e-9 * peak);
114+ // at a free corner both curvatures and the twist vanish, so every stress there is zero
115+ const vm = result.frames.vm;
116+ for (const [i, j] of [
117+ [0, 0],
118+ [0, nx - 1],
119+ [ny - 1, 0],
120+ [ny - 1, nx - 1],
121+ ])
122+ expect(vm.data[i + j * ny + k * nx * ny]).toBe(0);
123+ // and the edges move, unlike a supported plate
124+ let edgePeak = 0;
125+ for (let j = 0; j < nx; j++)
126+ edgePeak = Math.max(edgePeak, Math.abs(data[0 + j * ny + k * nx * ny]));
127+ expect(edgePeak).toBeGreaterThan(0.05 * peak);
128+ });
129+
130+ it("diverges just above the CFL limit, as the supported plate does", async () => {
131+ const spec = testSpec();
132+ spec.bc = "free";
133+ spec.fd.cfl = 1.04;
134+ const { result } = await solve(spec, "fd");
135+ expect(result.info.stable).toBe(false);
136+ expect(result.probe.w.every(Number.isFinite)).toBe(true);
137+ });
138+
139+ it("agrees between the implicit and explicit schemes", async () => {
140+ const spec = testSpec();
141+ spec.bc = "free";
142+ spec.implicit = {
143+ nx: spec.fd.nx,
144+ ny: spec.fd.ny,
145+ duration: spec.fd.duration,
146+ steps_per_sample: 2,
147+ theta: 0.25,
148+ };
149+ const imp = await solve(spec, "implicit");
150+ const exp = await solve(spec, "fd");
151+ expect(imp.result.info.stable).toBe(true);
152+ expect(imp.result.info.extra.nunknowns).toBe(spec.fd.nx * spec.fd.ny);
153+ const cmp = await runner.run({
154+ bundle: buildCompareBundle(),
155+ inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
156+ });
157+ expect(cmp.ok).toBe(true);
158+ if (!cmp.ok) return;
159+ const c = parseCompareResult(cmp.value);
160+ expect(c.w_rel_l2).toBeLessThan(0.03);
161+ });
162+
163+ it("is refused by the modal method with a legible message", async () => {
164+ const spec = testSpec();
165+ spec.bc = "free";
166+ await expect(solve(spec, "modal")).rejects.toThrow(/simply supported/);
167+ });
168+});