2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 1function [u, v, acc] = plate_sdof(omega, sigma, p, dt)
2% PLATE_SDOF Exact response of a damped oscillator to piecewise-linear forcing.
3%
4% [u, v, acc] = plate_sdof(omega, sigma, p, dt) integrates
5%
6% u'' + 2 sigma u' + omega^2 u = p(t)
7%
8% from rest, where p is sampled every dt seconds and varies linearly between
9% samples. Returns displacement u, velocity v and acceleration acc at the
10% sample times. The update is the exact solution over each interval
11% (Nigam and Jennings, 1969; Chopra, Dynamics of Structures, sec. 5.2), so
12% there is no time-step stability limit and no phase error: the only
13% approximation is the linear interpolation of p between samples.
14%
15% Requires sigma < omega (underdamped).
17if sigma >= omega
18 error('plate_sdof: the oscillator must be underdamped (sigma = %g >= omega = %g)', sigma, omega);
19end
20wd = sqrt(omega^2 - sigma^2); % damped natural frequency
21e = exp(-sigma * dt);
22s = sin(wd * dt);
23c = cos(wd * dt);
24w2 = omega^2;
25zr = sigma / wd; % zeta / sqrt(1 - zeta^2)
27% Coefficients of u(n+1) = A u(n) + B v(n) + C p(n) + D p(n+1)
28% v(n+1) = A2 u(n) + B2 v(n) + C2 p(n) + D2 p(n+1)
29A = e * (zr * s + c);
30B = e * s / wd;
31C = (2 * sigma / (w2 * dt) + e * (((w2 - 2 * sigma^2) / (w2 * wd * dt) - zr) * s - (1 + 2 * sigma / (w2 * dt)) * c)) / w2;
32D = (1 - 2 * sigma / (w2 * dt) + e * ((2 * sigma^2 - w2) / (w2 * wd * dt) * s + 2 * sigma / (w2 * dt) * c)) / w2;
33A2 = -e * w2 / wd * s;
34B2 = e * (c - zr * s);
35C2 = (-1 / dt + e * ((w2 / wd + sigma / (dt * wd)) * s + c / dt)) / w2;
36D2 = (1 - e * (zr * s + c)) / (w2 * dt);
38n = numel(p);
39u = zeros(size(p));
40v = zeros(size(p));
41for k = 1:n - 1
42 u(k + 1) = A * u(k) + B * v(k) + C * p(k) + D * p(k + 1);
43 v(k + 1) = A2 * u(k) + B2 * v(k) + C2 * p(k) + D2 * p(k + 1);
44end
45acc = p - 2 * sigma * v - w2 * u;
46end