1function C = plate_fem_contact(params, P, freqs, zeta, phi_hit, t)
2% PLATE_FEM_CONTACT The hammer and the plate as one system during the collision.
3%
4% C = plate_fem_contact(params, P, freqs, zeta, phi_hit, t) integrates the
5% coupled equations of plate_fem_rhs from the hammer touching the plate
6% at its speed (x_h = 0, v_h = speed, the plate at rest) until the two
7% separate, delta = x_h - phi' q falling through zero. All the modes take
8% part, the three rigid-body ones included, so that the very light plate
9% can move away from the hammer while they are in contact: this is what
10% makes the force lower and the contact longer than the rigid-target
11% pulse of plate_rr_contact, and the hammer's momentum transfer smaller.
12%
13% The integrator is classical fourth-order Runge-Kutta with a fixed step,
14% the Hertz duration estimate divided by 2000; at 0.05 us the highest
15% mode kept (below 0.47 fs) advances by less than a hundredth of a radian
16% per step. Separation is located by bisection within the step in which
17% delta changes sign. The stored steps and their derivatives give a cubic
18% Hermite interpolant, from which the states at the audio samples inside
19% the contact (C.qs, C.vs, C.fs at C.ts), and the histories on a 600-point
20% grid (C.tg: the force C.fg, the hammer C.hammer_x and C.hammer_v, the
21% plate under it C.plate_w and C.plate_v, and the modal states C.qg, C.vg
22% for the frames) are read. The impulse is the trapezoidal integral of
23% the force; the effective restitution is the relative separation speed
24% over the approach speed.
26hm = params.hammer;
27H = P.hertz;
28omega = 2 * pi * freqs(:);
29zeta = zeta(:);
30phi = phi_hit(:)';
31nm = numel(omega);
32twozw = 2 * zeta .* omega;
33k = H.contact_k;
34hc = H.hunt_crossley;
35mass = hm.mass;
36est = H.duration;
37nsub = 2000;
38h = est / nsub;
39nmax = 8 * nsub;
41y = zeros(2 * nm + 2, 1);
42y(end - 1) = 1e-12;
43y(end) = hm.speed;
44Ys = zeros(2 * nm + 2, nmax + 1);
45Fs = zeros(2 * nm + 2, nmax + 1);
46Ys(:, 1) = y;
47Fs(:, 1) = plate_fem_rhs(y, phi, nm, omega, twozw, k, hc, mass);
48found = false;
49tc = 0;
50yc = y;
51nsteps = 0;
52for step = 1:nmax
53 [ynew, fnew] = rk4_step(y, h, phi, nm, omega, twozw, k, hc, mass);
54 sep = ynew(end - 1) - phi * ynew(1:nm);
55 if step * h >= 1e-10 && sep < 0
56 % bisection on the step length from the last stored state
57 lo = 0;
58 hi = h;
59 for b = 1:60
60 mid = 0.5 * (lo + hi);
61 ym = rk4_step(y, mid, phi, nm, omega, twozw, k, hc, mass);
62 if ym(end - 1) - phi * ym(1:nm) < 0
63 hi = mid;
64 else
65 lo = mid;
66 end
67 if hi - lo < 1e-18
68 break
69 end
70 end
71 tc = (step - 1) * h + hi;
72 yc = rk4_step(y, hi, phi, nm, omega, twozw, k, hc, mass);
73 found = true;
74 nsteps = step - 1;
75 break
76 end
77 y = ynew;
78 Ys(:, step + 1) = y;
79 Fs(:, step + 1) = fnew;
80 nsteps = step;
81end
82if ~found
83 error('plate_fem_contact: the hammer did not separate within %.3g us (8 Hertz durations)', 8 * est * 1e6);
84end
85Ys = Ys(:, 1:nsteps + 1);
86Fs = Fs(:, 1:nsteps + 1);
87fc = plate_fem_rhs(yc, phi, nm, omega, twozw, k, hc, mass);
89% states at the audio samples inside the contact and on the plotting grid
90ts = t(t < tc);
91Yq = hermite(ts, h, tc, Ys, Fs, yc, fc);
92tg = linspace(0, tc, 600);
93Yg = hermite(tg, h, tc, Ys, Fs, yc, fc);
94fs_ = zeros(1, numel(ts));
95for j = 1:numel(ts)
96 fs_(j) = plate_fem_force(Yq(:, j), phi, nm, k, hc);
97end
98fg = zeros(1, 600);
99for j = 1:600
100 fg(j) = plate_fem_force(Yg(:, j), phi, nm, k, hc);
101end
103C.model = 'coupled Hunt-Crossley/Hertz contact with the moving plate';
104C.tc = tc;
105C.q0 = yc(1:nm)';
106C.v0 = yc(nm + 1:2 * nm)';
107C.ts = ts;
108C.qs = Yq(1:nm, :);
109C.vs = Yq(nm + 1:2 * nm, :);
110C.fs = fs_;
111C.tg = tg;
112C.qg = Yg(1:nm, :);
113C.vg = Yg(nm + 1:2 * nm, :);
114C.fg = fg;
115C.hammer_x = Yg(end - 1, :);
116C.hammer_v = Yg(end, :);
117C.plate_w = phi * C.qg;
118C.plate_v = phi * C.vg;
119C.impulse = trapz(tg, fg);
120C.peak_force = max(fg);
121C.hammer_exit_v = yc(end);
122C.plate_exit_v = phi * yc(nm + 1:2 * nm);
123C.relative_separation_v = C.hammer_exit_v - C.plate_exit_v;
124C.effective_restitution = -C.relative_separation_v / max(hm.speed, 1e-9);
125C.steps = nsteps;
126C.notes = {sprintf('coupled contact: %.4g us long, %.4g N peak, %.4g mN s; the hammer leaves at %.3g m/s into the plate, the plate under it at %.3g m/s (%d Runge-Kutta steps of %.3g ns)', ...
127 tc * 1e6, C.peak_force, C.impulse * 1e3, C.hammer_exit_v, C.plate_exit_v, nsteps, h * 1e9)};
128end
130function [ynew, fnew] = rk4_step(y, h, phi, nm, omega, twozw, k, hc, mass)
131k1 = plate_fem_rhs(y, phi, nm, omega, twozw, k, hc, mass);
132k2 = plate_fem_rhs(y + 0.5 * h * k1, phi, nm, omega, twozw, k, hc, mass);
133k3 = plate_fem_rhs(y + 0.5 * h * k2, phi, nm, omega, twozw, k, hc, mass);
134k4 = plate_fem_rhs(y + h * k3, phi, nm, omega, twozw, k, hc, mass);
135ynew = y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
136fnew = plate_fem_rhs(ynew, phi, nm, omega, twozw, k, hc, mass);
137end
139function Y = hermite(times, h, tc, Ys, Fs, yc, fc)
140% cubic Hermite interpolation of the stored trajectory at the given times in [0, tc]
141n = size(Ys, 2) - 1; % whole steps stored
142Y = zeros(size(Ys, 1), numel(times));
143for j = 1:numel(times)
144 tj = min(max(times(j), 0), tc);
145 i = floor(tj / h);
146 if i < n
147 t0 = i * h;
148 hh = h;
149 y0 = Ys(:, i + 1); f0 = Fs(:, i + 1);
150 y1 = Ys(:, i + 2); f1 = Fs(:, i + 2);
151 else
152 t0 = n * h;
153 hh = tc - t0;
154 y0 = Ys(:, n + 1); f0 = Fs(:, n + 1);
155 y1 = yc; f1 = fc;
156 end
157 if hh <= 0
158 Y(:, j) = y0;
159 continue
160 end
161 s = (tj - t0) / hh;
162 h00 = 2 * s^3 - 3 * s^2 + 1;
163 h10 = s^3 - 2 * s^2 + s;
164 h01 = -2 * s^3 + 3 * s^2;
165 h11 = s^3 - s^2;
166 Y(:, j) = h00 * y0 + h10 * hh * f0 + h01 * y1 + h11 * hh * f1;
167end
168end