2004981Plate dynamics lab 2: a hammer on a hung plate by Rayleigh-Ritz and by Mindlin-Reissner finite elementsJeremy Magland 1function M = plate_rr_modes(params)
2% PLATE_RR_MODES Modes of the hung plate by the Rayleigh-Ritz method.
3%
4% M = plate_rr_modes(params) minimizes the Kirchhoff-Love bending energy
5%
6% U = D/2 integral ( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1 - nu) w_xy^2 ) dA
7%
8% over the tensor-Legendre basis of order params.rr.basis_order
9% (plate_rr_basis), with the kinetic energy rho h / 2 integral w_t^2 dA, both
10% integrated by Gauss-Legendre quadrature of params.rr.quadrature_order
11% points per side. The hole subtracts its area's mass at its position
12% (its effect on the stiffness is negligible at 0.8 mm) and the thread
13% adds the pendulum spring m g / L there. The generalized eigenproblem
14% K a = omega^2 Mm a is solved densely after a Cholesky reduction; the
15% eigenvectors are mass-normalized.
16%
17% Kept: modes above 20 Hz and below 0.47 fs, at most params.render.modes
18% of them. The pendulum mode of the thread and the two rigid rotations
19% are dropped here (this method has no rigid-body content; the
20% finite-element method keeps them during the collision).
21%
22% Returns M with method 'rr', freqs (1 x nm), vecs (n x nm) the Ritz
23% coefficients, order, support_part (the shape at the hole) and
24% volume_part (the integral of the shape over the plate) per mode, and the
25% shapes and curvatures on a display grid of params.rr.grid points per
26% side: shapes, kxx, kyy, kxy (ny x nx x nm; rows index y), grid_x, grid_y.
28t_start = tic;
29P = plate_setup(params);
30pl = params.plate;
31W = pl.width;
32H = pl.height;
33order = params.rr.basis_order;
34nq = params.rr.quadrature_order;
35mu = params.material.rho * pl.thickness;
36nu = params.material.nu;
37notes = {};
39% quadrature over the plate
40[qx, wx] = plate_gauss_legendre(nq);
41[XX, YY] = ndgrid(qx, qx);
42WW = wx(:) * wx(:)';
43ww = WW(:)' * W * H / 4;
44[b, bxx, byy, bxy] = plate_rr_basis(XX(:) * W / 2, YY(:) * H / 2, order, W, H);
46Mm = mu * ((b .* ww) * b');
47K = P.D * ((bxx .* ww) * bxx' + (byy .* ww) * byy' ...
48 + nu * ((bxx .* ww) * byy' + (byy .* ww) * bxx') ...
49 + 2 * (1 - nu) * ((bxy .* ww) * bxy'));
51% the hole's mass and the thread's spring at the hole
52ph = plate_rr_basis(P.hole_x, P.hole_y, order, W, H);
53Mm = Mm - mu * P.hole_area * (ph * ph');
54K = K + P.thread_k * (ph * ph');
56% generalized symmetric eigenproblem by Cholesky reduction
57Mm = (Mm + Mm') / 2;
58K = (K + K') / 2;
59R = chol(Mm);
60Cm = R' \ K / R;
61Cm = (Cm + Cm') / 2;
62[V, Dm] = eig(Cm);
63[vals, order_] = sort(diag(Dm)');
64V = V(:, order_);
65vecs = R \ V;
67keep = vals > 1e-5;
68vals = vals(keep);
69vecs = vecs(:, keep);
70freqs = sqrt(max(vals, 0)) / (2 * pi);
71audible = (freqs > 20) & (freqs < P.f_max);
72n_audible = nnz(audible);
73count = min(params.render.modes, n_audible);
74sel = find(audible);
75sel = sel(1:count);
76freqs = freqs(sel);
77vecs = vecs(:, sel);
78nm = numel(freqs);
79if count < n_audible
80 notes{end + 1} = sprintf('%d modes lie below 0.47 fs; the first %d are kept (params.render.modes)', n_audible, count);
81end
83% participation of the hole and of the plate's volume (surface quadrature of order max(36, nq))
84support_part = (vecs' * ph)';
85ns = max(36, nq);
86[qs, ws] = plate_gauss_legendre(ns);
87[XS, YS] = ndgrid(qs, qs);
88WS = ws(:) * ws(:)';
89area_w = WS(:) * W * H / 4;
90bs = plate_rr_basis(XS(:) * W / 2, YS(:) * H / 2, order, W, H);
91volume_part = ((vecs' * bs) * area_w)';
93% shapes and curvatures on the display grid
94ng = params.rr.grid;
95grid_x = linspace(-W / 2, W / 2, ng);
96grid_y = linspace(-H / 2, H / 2, ng);
97[YG, XG] = ndgrid(grid_y, grid_x);
98[bg, bgxx, bgyy, bgxy] = plate_rr_basis(XG(:), YG(:), order, W, H);
99shapes = reshape(bg' * vecs, ng, ng, nm);
100kxx = reshape(bgxx' * vecs, ng, ng, nm);
101kyy = reshape(bgyy' * vecs, ng, ng, nm);
102kxy = reshape(bgxy' * vecs, ng, ng, nm);
104n = (order + 1)^2;
105M.method = 'rr';
106M.freqs = freqs;
107M.vecs = vecs;
108M.order = order;
109M.support_part = support_part;
110M.volume_part = volume_part;
111M.shapes = shapes;
112M.kxx = kxx;
113M.kyy = kyy;
114M.kxy = kxy;
115M.grid_x = grid_x;
116M.grid_y = grid_y;
117M.ndof = n;
118M.n_rigid = 0;
119M.n_found = numel(vals);
120M.plate_mass = P.plate_mass;
121M.size_note = sprintf('%d Legendre coefficients (order %d), %d x %d quadrature, dense eigensolve', n, order, nq, nq);
122M.wall_s = toc(t_start);
123M.notes = notes;
124end