1function result = plate_fd(params)
2% PLATE_FD Explicit finite differences for the damped Kirchhoff-Love plate.
3%
4% result = plate_fd(params) advances
5%
6% w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
7%
8% on a uniform nx x ny grid (nodes include the edges) with the scheme
9%
10% (1 + sigma0 dt) w+ = 2 w - (1 - sigma0 dt) w- - dt^2 kappa^2 L(L w)
11% + 2 sigma1 dt (L w - L w-) + dt^2 F s / mu
12%
13% where w-, w, w+ are consecutive time levels and L is the five-point
14% Laplacian. Applying L twice gives the 13-point biharmonic stencil. The
15% 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.
18%
19% The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
20% is unstable for cfl > 1 and the run stops as soon as the field diverges,
21% reporting the step at which it did (info.diverged_step).
22%
23% Approximations. Second-order truncation error in dx, dy, dt; the probe is
24% snapped to the nearest grid node; stress fields are evaluated only at the
25% saved frames, so the envelope is a sample of the true peak; the audio is
26% the Rayleigh integral box-averaged onto audio samples.
27%
28% See plate_result_schema for the fields of result.
30t_start = tic;
31P = plate_setup(params);
32plan = plate_fd_plan(params);
33nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
34x = linspace(0, P.plate.Lx, nx);
35y = linspace(0, P.plate.Ly, ny);
36[X, Y] = meshgrid(x, y); % w(i, j) lives at (x(j), y(i))
37fs = P.output.fs;
38notes = {};
39if P.fd.cfl > 1
40 notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
41end
43s = plate_force_footprint(P.excitation, X, Y, dx * dy);
44tt = (0:nsteps) * dt;
45F = plate_force_time(P.excitation, tt);
47% Probe: nearest interior node
48[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
49[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), ny - 1);
50if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
51 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
52end
54% Frames at the steps nearest to nframes equally spaced times in [0, duration]
55frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));
56nf = numel(frame_steps);
57Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
58tf = zeros(1, nf);
60% Listener: every grid node is a patch; the acceleration of a step is added into
61% the audio sample that contains it, delayed by the node's travel time.
62L = plate_listener(P.listener, X, Y, dx * dy, fs);
63bins = L.delay(:) + 1;
64wgt = L.weight(:);
65na = ceil(nsteps * dt * fs) + L.nbins + 1;
66Pbins = zeros(L.nbins, na);
68% State
69w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
70w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
71c_old = 1 - P.damping.sigma0 * dt;
72c_new = 1 / (1 + P.damping.sigma0 * dt);
73k2dt2 = P.kappa^2 * dt^2;
74s1dt2 = 2 * P.damping.sigma1 * dt;
75fdt2 = dt^2 / P.mu;
76blowup = 1e3 * P.plate.h; % deflection beyond which the run is declared unstable
77diverged_step = 0;
78kf = 1;
80for 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;
85 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;
88 % Probe displacement and curvatures from the 3 x 3 neighbourhood
89 w_p(nstep + 1) = wn(ip, jp);
90 wxx_p(nstep + 1) = (wn(ip, jp + 1) - 2 * wn(ip, jp) + wn(ip, jp - 1)) / dx^2;
91 wyy_p(nstep + 1) = (wn(ip + 1, jp) - 2 * wn(ip, jp) + wn(ip - 1, jp)) / dy^2;
92 wxy_p(nstep + 1) = (wn(ip + 1, jp + 1) - wn(ip + 1, jp - 1) - wn(ip - 1, jp + 1) + wn(ip - 1, jp - 1)) / (4 * dx * dy);
94 % Radiated pressure: acceleration at this step, weighted and delayed per node
95 acc = (wn - 2 * w + wm) / dt^2;
96 ia = floor(nstep * dt * fs) + 1;
97 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [L.nbins 1]) * (dt * fs);
99 if kf <= nf && nstep == frame_steps(kf)
100 [wxx, wyy, wxy] = plate_curvature(plate_pad(wn, P.bc), dx, dy);
101 S = plate_stress(P, wxx, wyy, wxy);
102 Wf(:, :, kf) = wn; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
103 kf = kf + 1;
104 end
106 peak = max(abs(wn(:)));
107 if ~(peak < blowup) % also catches NaN and Inf
108 diverged_step = nstep;
109 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
110 break
111 end
112 wm = w; w = wn; Lw_m = Lw;
113end
115% Truncate after a blow-up so the outputs stay finite
116if diverged_step > 0
117 nkeep = diverged_step; % keep the last finite level
118 tt = tt(1:nkeep); w_p = w_p(1:nkeep); wxx_p = wxx_p(1:nkeep); wyy_p = wyy_p(1:nkeep); wxy_p = wxy_p(1:nkeep);
119 keepf = frame_steps < diverged_step;
120 Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
121end
123% Audio: sum the delay bins, then trim to the simulated duration
124p_audio = zeros(1, na);
125for d = 1:L.nbins
126 p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
127end
128n_audio = min(na, ceil(numel(tt) * dt * fs));
129p_audio = p_audio(1:n_audio);
131Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
132[fw, Aw] = plate_spectrum(tt, w_p);
133[~, Avm] = plate_spectrum(tt, Sp.vm);
135result.method = 'fd';
136result.grid = struct('x', x, 'y', y);
137result.t = tt;
138result.probe = struct('x', x(jp), 'y', y(ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
139result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
140if isempty(tf)
141 result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
142else
143 result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
144end
145result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
146result.audio = struct('fs', fs, 'p', p_audio, 't0', L.t0);
147result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', plan.dt_max, 'cfl', P.fd.cfl, ...
148 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
149 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
150end