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). 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).
21%
22% The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
23% is unstable for cfl > 1 and the run stops as soon as the field diverges,
24% reporting the step at which it did (info.diverged_step).
25%
26% Approximations. Second-order truncation error in dx, dy, dt; the probe is
27% snapped to the nearest grid node; stress fields are evaluated only at the
28% saved frames, so the envelope is a sample of the true peak; the audio is
29% the Rayleigh integral box-averaged onto audio samples.
30%
31% See plate_result_schema for the fields of result.
33t_start = tic;
34P = plate_setup(params);
35plan = plate_fd_plan(params);
36nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
37x = linspace(0, P.plate.Lx, nx);
38y = linspace(0, P.plate.Ly, ny);
39[X, Y] = meshgrid(x, y); % w(i, j) lives at (x(j), y(i))
40fs = P.output.fs;
41notes = {};
42if P.fd.cfl > 1
43 notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
44end
46s = plate_force_footprint(P.excitation, X, Y, dx * dy);
47tt = (0:nsteps) * dt;
48F = plate_force_time(P.excitation, tt);
50% Probe: nearest interior node
51[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
52[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), ny - 1);
53if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
54 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
55end
57% Frames at the steps nearest to nframes equally spaced times in [0, duration]
58frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));
59nf = numel(frame_steps);
60Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
61tf = zeros(1, nf);
63% Listener: every grid node is a patch; the acceleration of a step is added into
64% the audio sample that contains it, delayed by the node's travel time.
65L = plate_listener(P.listener, X, Y, dx * dy, fs);
66bins = L.delay(:) + 1;
67wgt = L.weight(:);
68na = ceil(nsteps * dt * fs) + L.nbins + 1;
69Pbins = zeros(L.nbins, na);
71% Free edges: sparse operators on every node, and the rigid-body projector
72free = strcmp(P.bc, 'free');
73if 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';
78end
80% State
81w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
82w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
83c_old = 1 - P.damping.sigma0 * dt;
84c_new = 1 / (1 + P.damping.sigma0 * dt);
85k2dt2 = P.kappa^2 * dt^2;
86s1dt2 = 2 * P.damping.sigma1 * dt;
87fdt2 = dt^2 / P.mu;
88blowup = 1e3 * P.plate.h; % deflection beyond which the run is declared unstable
89diverged_step = 0;
90kf = 1;
92for nstep = 1:nsteps
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
102 wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
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
109 % Probe displacement and curvatures from the 3 x 3 neighbourhood
110 w_p(nstep + 1) = wn(ip, jp);
111 wxx_p(nstep + 1) = (wn(ip, jp + 1) - 2 * wn(ip, jp) + wn(ip, jp - 1)) / dx^2;
112 wyy_p(nstep + 1) = (wn(ip + 1, jp) - 2 * wn(ip, jp) + wn(ip - 1, jp)) / dy^2;
113 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);
115 % Radiated pressure: acceleration at this step, weighted and delayed per node
116 acc = (wn - 2 * w + wm) / dt^2;
117 ia = floor(nstep * dt * fs) + 1;
118 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [L.nbins 1]) * (dt * fs);
120 if kf <= nf && nstep == frame_steps(kf)
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
126 S = plate_stress(P, wxx, wyy, wxy);
127 Wf(:, :, kf) = wn; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
128 kf = kf + 1;
129 end
131 peak = max(abs(wn(:)));
132 if ~(peak < blowup) % also catches NaN and Inf
133 diverged_step = nstep;
134 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
135 break
136 end
137 wm = w; w = wn; Lw_m = Lw;
138end
140% Truncate after a blow-up so the outputs stay finite
141if diverged_step > 0
142 nkeep = diverged_step; % keep the last finite level
143 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);
144 keepf = frame_steps < diverged_step;
145 Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
146end
148% Audio: sum the delay bins, then trim to the simulated duration
149p_audio = zeros(1, na);
150for d = 1:L.nbins
151 p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
152end
153n_audio = min(na, ceil(numel(tt) * dt * fs));
154p_audio = p_audio(1:n_audio);
156Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
157[fw, Aw] = plate_spectrum(tt, w_p);
158[~, Avm] = plate_spectrum(tt, Sp.vm);
160result.method = 'fd';
161result.grid = struct('x', x, 'y', y);
162result.t = tt;
163result.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);
164result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
165if isempty(tf)
166 result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
167else
168 result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
169end
170result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
171result.audio = struct('fs', fs, 'p', p_audio, 't0', L.t0);
172result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', plan.dt_max, 'cfl', P.fd.cfl, ...
173 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
174 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
175end