1function [Kg, ps] = plate_prestress(P, dom, mnt, ops)
2% PLATE_PRESTRESS In-plane stress of the plate's own weight, and its effect on bending.
3%
4% [Kg, ps] = plate_prestress(P, dom, mnt, ops) solves the plane-stress
5% problem of the plate loaded by the in-plane part of its weight,
6% rho h g (gx, gy) per unit area, with the in-plane constraints of the mount
7% (the holes a thread passes through, the jaws of a vise, the pinched patch),
8% and returns the sparse N x N operator Kg of the membrane term
9%
10% U_N = 1/2 integral( N_xx w_x^2 + N_yy w_y^2 + 2 N_xy w_x w_y ) dA
11%
12% that the in-plane stress resultants N add to the bending energy. This
13% is how gravity acts on the vibration of a hanging plate: a tilted element
14% under tension is pulled back, and for a rigid rotation of the plate about
15% its hole the term integrates to the compound-pendulum stiffness m g d
16% exactly (d the distance from the hole to the centre of mass), so the
17% swing of a hung plate and the small stiffening of its lowest flexural
18% modes come out of the same operator. Parts of the plate above the hole
19% are in compression and the term softens them slightly. Without an
20% in-plane gravity component (a horizontal plate) Kg is empty.
21%
22% Discretization. The in-plane displacement (u, v) lives on the nodes.
23% The direct strains are first differences along the grid edges,
24% e_xx = (u_R - u_L) / dx on the x-edges and e_yy on the y-edges, the shear
25% strain gamma = u_y + v_x and the coupling e_xx e_yy are taken at the cell
26% centres from the averaged corner differences (ops from plate_bending).
27% The strain energy
28%
29% U = h/2 sum_edges C11 w_e e_xx^2 + h/2 sum_edges C11 w_e e_yy^2
30% + h sum_cells A_c ( C12 e_xx e_yy + G/2 gamma^2 ),
31% C11 = E / (1 - nu^2), C12 = nu C11, G = E / (2 (1 + nu)),
32%
33% with the edge weights w_e their share of cell area, has no spurious zero
34% modes: only the two translations and the rotation are strain-free. The
35% stress resultants are the derivatives of U with respect to the strains,
36% n_xx = dU/de_xx on each x-edge, n_yy on each y-edge, n_xy = dU/dgamma on
37% each cell (each N times its area), and Kg is assembled from the same
38% difference operators applied to w:
39%
40% Kg = Bx' diag(n_xx) Bx + By' diag(n_yy) By + Bxc' diag(n_xy) Byc + Byc' diag(n_xy) Bxc.
41%
42% Because both U and U_N use the same strain operators, the identity
43% r' Kg r = m g d for the rigid tilt r holds to rounding, not just to
44% truncation error.
45%
46% Constraints (mnt.fix) are imposed by a change of variables: a node with
47% one fixed direction keeps its displacement along the other, a node with
48% two keeps none. A plate on one thread can still rotate about the hole;
49% the load has no moment about it (the plate hangs with its centre of mass
50% below the hole, plate_mount), so fixing the tangential displacement of
51% the node farthest from the hole removes the rotation at zero reaction.
52%
53% ps returns the solution for inspection: u, v (ny x nx), the reactions at
54% the fixed nodes and the tension of each thread, and the node-averaged
55% stress resultants Nxx, Nyy, Nxy [N/m].
57nx = dom.nx; ny = dom.ny;
58N = nx * ny;
59ps = struct('solved', false);
60if P.g == 0 || (mnt.gx == 0 && mnt.gy == 0)
61 Kg = sparse(N, N);
62 return
63end
64idx = ops.idx; n = ops.n;
65Bx = ops.Bx; By = ops.By; Bxc = ops.Bxc; Byc = ops.Byc; Cx = ops.Cx; Cy = ops.Cy;
66wx = ops.wx; wy = ops.wy; area = ops.area; ncell = ops.ncell;
67nex = numel(wx); ney = numel(wy);
68E = P.material.E; nu = P.material.nu; h = P.plate.h;
69C11 = E / (1 - nu^2); C12 = nu * C11; G = E / (2 * (1 + nu));
71% Stiffness on [u; v]
72Kuu = h * C11 * (Bx' * spdiags(wx, 0, nex, nex) * Bx);
73Kvv = h * C11 * (By' * spdiags(wy, 0, ney, ney) * By);
74Kuv = h * C12 * area * (Bxc' * Byc);
75Gam = [Byc, Bxc]; % gamma = u_y + v_x at the cells
76Kps = [Kuu, Kuv; Kuv', Kvv] + h * G * area * (Gam' * Gam);
78% Load: the in-plane weight of each node
79Mi = P.mu * ops.mu;
80f = [Mi * P.g * mnt.gx; Mi * P.g * mnt.gy];
82% Constraints, in the plate's node numbering
83fix = mnt.fix;
84loc = zeros(N, 1); loc(idx) = (1:n)';
85fixn = loc(fix(:, 1));
86if any(fixn == 0)
87 error('plate_prestress: a constraint sits on a node off the plate');
88end
89fix_mount = unique(fixn); % the mount's own constraints, for the reactions
90if mnt.free_rotation
91 % the tangential displacement of the node farthest from the pivot is fixed
92 pk = fixn(1);
93 X = dom.X(idx); Y = dom.Y(idx);
94 d2 = (X - X(pk)).^2 + (Y - Y(pk)).^2;
95 d2(ops.unk == 0) = -1;
96 [~, far] = max(d2);
97 t = [-(Y(far) - Y(pk)), X(far) - X(pk)];
98 t = t / norm(t);
99 fixn = [fixn; far];
100 fix = [fix; idx(far) t];
101end
102% Rotation matrix T from the reduced unknowns to [u; v]
103one = zeros(n, 2); cnt = zeros(n, 1);
104for q = 1:numel(fixn)
105 k = fixn(q);
106 cnt(k) = cnt(k) + 1;
107 one(k, :) = fix(q, 2:3) / norm(fix(q, 2:3));
108end
109if any(cnt > 2)
110 error('plate_prestress: more than two constraints on one node');
111end
112both = cnt == 2;
113single = cnt == 1;
114none = cnt == 0;
115kn = find(none); ks = find(single);
116nq = 2 * numel(kn) + numel(ks);
117ti = [kn; n + kn; ks; n + ks];
118tj = [(1:numel(kn))'; numel(kn) + (1:numel(kn))'; 2 * numel(kn) + (1:numel(ks))'; 2 * numel(kn) + (1:numel(ks))'];
119tv = [ones(2 * numel(kn), 1); -one(ks, 2); one(ks, 1)]; % the direction perpendicular to the fixed one
120T = sparse(ti, tj, tv, 2 * n, nq);
121Kred = T' * Kps * T;
122fred = T' * f;
123q = Kred \ fred;
124U = T * q;
125u = U(1:n); v = U(n + 1:end);
126if any(~isfinite(U))
127 error('plate_prestress: the in-plane solve failed; the mount does not hold the plate in its plane');
128end
129react = Kps * U - f; % force the constraints apply, [x; y] per node
131% Stress resultants times area at the edges and cells
132exx = Bx * u; eyy = By * v;
133exxc = Bxc * u; eyyc = Byc * v; gam = Gam * U;
134nxx = h * (C11 * (wx .* exx) + C12 * area * (Cx' * eyyc));
135nyy = h * (C11 * (wy .* eyy) + C12 * area * (Cy' * exxc));
136nxy = h * G * area * gam;
137Kgu = Bx' * spdiags(nxx, 0, nex, nex) * Bx + By' * spdiags(nyy, 0, ney, ney) * By ...
138 + Bxc' * spdiags(nxy, 0, ncell, ncell) * Byc + Byc' * spdiags(nxy, 0, ncell, ncell) * Bxc;
139Pz = spdiags(ops.unk, 0, n, n);
140Kgu = Pz * Kgu * Pz;
141[ki, kj, kv] = find(Kgu);
142Kg = sparse(idx(ki), idx(kj), kv, N, N);
144% Node-averaged resultants for display: each edge's value shared by its two
145% nodes, each cell's by its four corners, divided by the areas they cover
146ps.solved = true;
147uf = zeros(ny, nx); vf = zeros(ny, nx);
148uf(idx) = u; vf(idx) = v;
149ps.u = uf; ps.v = vf;
150Ex = abs(Bx) * (dom.dx / 2); % 1/2 at each node of an x-edge
151Ey = abs(By) * (dom.dy / 2);
152sx = Ex' * nxx; ax = Ex' * wx;
153sy = Ey' * nyy; ay = Ey' * wy;
154ic = (1:ncell)';
155Qc = sparse([ic; ic; ic; ic], ops.corners(:), 0.25 * ones(4 * ncell, 1), ncell, n);
156sxy = Qc' * nxy; ac = Qc' * (area * ones(ncell, 1));
157Nxx = zeros(ny, nx); Nyy = zeros(ny, nx); Nxy = zeros(ny, nx);
158Nxx(idx) = sx ./ max(ax, eps); Nyy(idx) = sy ./ max(ay, eps); Nxy(idx) = sxy ./ max(ac, eps);
159Nxx(~dom.inside) = NaN; Nyy(~dom.inside) = NaN; Nxy(~dom.inside) = NaN;
160ps.Nxx = Nxx; ps.Nyy = Nyy; ps.Nxy = Nxy;
161ps.nxx = nxx; ps.nyy = nyy; ps.nxy = nxy; % the resultants times area, per edge and cell: the factors of Kg
162% Reactions and thread tensions: the force along -ed (up) at each fixed node
163rx = react(1:n); ry = react(n + 1:end);
164ps.react_x = rx; ps.react_y = ry;
165ps.up = -(rx * mnt.ed(1) + ry * mnt.ed(2)); % upward force each constrained node carries
166ps.fixed_nodes = idx(fix_mount);
167ps.mass_check = sum(ps.up(fix_mount)) / (mnt.mass * P.g); % 1 when the constraints carry the weight
168ps.energy = 0.5 * (U' * Kps * U);
169end