2004981Plate dynamics lab 2: a hammer on a hung plate by Rayleigh-Ritz and by Mindlin-Reissner finite elementsJeremy Magland 1function [ke, me] = plate_fem_element(dx, dy, params, P)
2% PLATE_FEM_ELEMENT Stiffness and mass of one Q4 Mindlin-Reissner element.
3%
4% [ke, me] = plate_fem_element(dx, dy, params, P) integrates, with 2 x 2
5% Gauss points, the bending energy of the rotations and the transverse-
6% shear energy of a rectangular element dx x dy with three unknowns per
7% node (w, theta_x, theta_y):
8%
9% U = 1/2 integral ( kappa' Db kappa + gamma' Ds gamma ) dA,
10% kappa = (d theta_x / dx, d theta_y / dy, d theta_x / dy + d theta_y / dx),
11% gamma = (w_x + theta_x, w_y + theta_y),
12% Db = D [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2], Ds = k G h I, k = 5/6,
13%
14% and the kinetic energy with the translational mass rho h and the rotary
15% inertia rho h^3 / 12. Full integration of the shear term is what the
16% simulator this lab follows does; on a thin plate it locks, which is why
17% that simulator needs its 80 x 80 mesh to converge (see the README).
18% Degrees of freedom are ordered (w, theta_x, theta_y) per node, nodes
19% anticlockwise from the lower-left corner.
21mat = params.material;
22h = params.plate.thickness;
23nu = mat.nu;
24bend = P.D * [1, nu, 0; nu, 1, 0; 0, 0, 0.5 * (1 - nu)];
25shear_modulus = mat.E / (2 * (1 + nu));
26kappa = 5 / 6;
27shear = kappa * shear_modulus * h * eye(2);
28rho = mat.rho;
29ke = zeros(12, 12);
30me = zeros(12, 12);
31gp = 1 / sqrt(3);
32detj = dx * dy / 4;
33for xi = [-gp, gp]
34 for eta = [-gp, gp]
35 [n, dxi, deta] = plate_fem_shape(xi, eta);
36 dnx = dxi * 2 / dx;
37 dny = deta * 2 / dy;
38 bb = zeros(3, 12);
39 bs = zeros(2, 12);
40 for a = 1:4
41 w = 3 * a - 2; rx = 3 * a - 1; ry = 3 * a;
42 bb(1, rx) = dnx(a);
43 bb(2, ry) = dny(a);
44 bb(3, rx) = dny(a);
45 bb(3, ry) = dnx(a);
46 bs(1, w) = dnx(a); bs(1, rx) = n(a);
47 bs(2, w) = dny(a); bs(2, ry) = n(a);
48 end
49 ke = ke + (bb' * bend * bb + bs' * shear * bs) * detj;
50 for a = 1:4
51 for b = 1:4
52 nn = n(a) * n(b) * detj;
53 me(3 * a - 2, 3 * b - 2) = me(3 * a - 2, 3 * b - 2) + rho * h * nn;
54 rotary = rho * h^3 / 12 * nn;
55 me(3 * a - 1, 3 * b - 1) = me(3 * a - 1, 3 * b - 1) + rotary;
56 me(3 * a, 3 * b) = me(3 * a, 3 * b) + rotary;
57 end
58 end
59 end
60end
61end