2004981Plate dynamics lab 2: a hammer on a hung plate by Rayleigh-Ritz and by Mindlin-Reissner finite elementsJeremy Magland 1function C = plate_rr_contact(params, P, freqs, zeta, phi_hit, t)
2% PLATE_RR_CONTACT The prescribed Hertz pulse and the modal states it leaves.
3%
4% C = plate_rr_contact(params, P, freqs, zeta, phi_hit, t) applies the
5% force history
6%
7% F(s) = I sin(pi s / T)^(3/2) / (T Gamma(5/4) / (sqrt(pi) Gamma(7/4))), 0 <= s <= T,
8%
9% of duration T and impulse I from plate_hertz (the rounded pulse of a
10% Hertz contact on a rigid target) to every mode, whose response is the
11% Duhamel integral with the damped impulse response
12%
13% h(u) = e^(-alpha u) sin(wd u) / wd, alpha = zeta omega, wd = omega sqrt(1 - zeta^2).
14%
15% The state at separation (C.q0, C.v0) is a 64-point Gauss-Legendre
16% quadrature of the integral over [0, T]; the states at the audio samples
17% inside the contact (C.qs, C.vs at C.ts) and on a 600-point plotting grid
18% (C.qg, C.vg at C.tg) are 32-point quadratures over [0, s]. The plate
19% does not move the hammer back in this model: the force is what it would
20% be on a rigid target, which is why the peak force and the impulse differ
21% from the coupled finite-element contact (see plate_fem_contact).
23H = P.hertz;
24T = H.duration;
25I = H.impulse;
26nrm = T * H.shape_integral;
27freqs = freqs(:)';
28omega = 2 * pi * freqs;
29zeta = zeta(:)';
30alpha = zeta .* omega;
31wd = omega .* sqrt(max(1 - zeta.^2, 1e-12));
32phi = phi_hit(:)';
33nm = numel(freqs);
35% terminal state
36[z, wz] = plate_gauss_legendre(64);
37s = 0.5 * T * (z + 1);
38ws = 0.5 * T * wz;
39force = I * sin(pi * s / T).^1.5 / nrm;
40u = T - s;
41E = exp(-alpha' * u);
42SW = sin(wd' * u);
43CW = cos(wd' * u);
44h = E .* SW ./ wd';
45hd = E .* (CW - (alpha' ./ wd') .* SW);
46q0 = phi .* (h * (ws .* force)')';
47v0 = phi .* (hd * (ws .* force)')';
49% states inside the contact: at the audio samples and on the plotting grid
50ts = t(t < T);
51[qs, vs, fs_] = rr_states(ts, T, I, nrm, alpha, wd, phi);
52tg = linspace(0, T, 600);
53[qg, vg, fg] = rr_states(tg, T, I, nrm, alpha, wd, phi);
55C.model = 'prescribed Hertz pulse (rigid-target contact)';
56C.tc = T;
57C.q0 = q0;
58C.v0 = v0;
59C.ts = ts;
60C.qs = qs;
61C.vs = vs;
62C.fs = fs_;
63C.tg = tg;
64C.qg = qg;
65C.vg = vg;
66C.fg = fg;
67C.plate_v = phi * vg;
68C.plate_w = phi * qg;
69C.hammer_v = zeros(1, 0);
70C.hammer_x = zeros(1, 0);
71C.impulse = I;
72C.peak_force = H.peak_force;
73C.hammer_exit_v = NaN;
74C.plate_exit_v = phi * v0';
75C.relative_separation_v = NaN;
76C.effective_restitution = params.hammer.restitution;
77C.notes = {sprintf('prescribed Hertz pulse: %.4g us long, %.4g N peak, %.4g mN s', T * 1e6, H.peak_force, I * 1e3)};
78end
80function [q, v, f] = rr_states(times, T, I, nrm, alpha, wd, phi)
81% modal states at the given times within the pulse, by quadrature over [0, time]
82nm = numel(alpha);
83nt = numel(times);
84q = zeros(nm, nt);
85v = zeros(nm, nt);
86f = zeros(1, nt);
87[z, wz] = plate_gauss_legendre(32);
88for j = 1:nt
89 tj = times(j);
90 if tj <= 0
91 continue
92 end
93 s = 0.5 * tj * (z + 1);
94 ws = 0.5 * tj * wz;
95 force = I * sin(pi * s / T).^1.5 / nrm;
96 u = tj - s;
97 E = exp(-alpha' * u);
98 SW = sin(wd' * u);
99 CW = cos(wd' * u);
100 h = E .* SW ./ wd';
101 hd = E .* (CW - (alpha' ./ wd') .* SW);
102 q(:, j) = phi' .* (h * (ws .* force)');
103 v(:, j) = phi' .* (hd * (ws .* force)');
104 if tj < T
105 f(j) = I * sin(pi * tj / T)^1.5 / nrm;
106 end
107end
108end