80c0b6dPlate dynamics lab: plates in realizable mounts under gravity, explicit finite differences on WebGPUJeremy Magland 1function st = plate_static(P, dom, mnt, A0, Mvec, unk)
2% PLATE_STATIC Static equilibrium of the plate in its mount under gravity.
3%
4% st = plate_static(P, dom, mnt, A0, Mvec, unk) solves
5%
6% (A0 + K_contact) w0 = M g_n + contact offsets
7%
8% for the deflection w0 of the plate at rest, with A0 the linear operator
9% (bending, self-weight prestress and the linear springs of the mount,
10% plate_model), M the node masses, g_n the normal component of gravity, and
11% K_contact the springs of the one-sided contacts that are compressed at
12% rest. Which contacts those are is found by an active-set iteration: all
13% contacts start engaged, any found in tension is released, and the solve is
14% repeated until the set no longer changes. A plate that nothing holds in
15% equilibrium (too few pads to rest on, or one standing on its edge that
16% buckles under its own weight) fails the Cholesky factorization, and the
17% error says so.
18%
19% The dynamics are then computed for the departure from w0 (plate_fd_loop):
20% the static forces cancel exactly, the contact forces are linear in the
21% departure until a contact opens, and single precision has only the motion
22% to represent.
23%
24% Fields.
25% w0 N x 1 static deflection (zero at nodes that are not unknowns)
26% actb, actt N x 1 logical: bottom and top contacts engaged at rest
27% gb, gt N x 1 compression of each contact at rest [m]: the bottom
28% contact pushes while gb - w > 0, the top while gt + w > 0
29% fstat_s, fstat_b, fstat_t N x 1 force at rest of the linear springs, the
30% bottom contacts and the top contacts [N]
31% Alin the operator linearized about rest, for the release solve
32% sag largest |w0| [m]
33% notes
35N = dom.nx * dom.ny;
36u = unk(:);
37kb = mnt.kb(:); kt = mnt.kt(:); ob = mnt.ob(:); ot = mnt.ot(:);
38ks = mnt.ks(:);
39fg = Mvec(:) .* u * P.g * mnt.gz;
40notes = {};
41actb = kb > 0; actt = kt > 0;
42w0 = zeros(N, 1);
43iu = find(u);
44if ~any(fg) && ~any(actb) && ~any(actt)
45 Alin = A0;
46else
47 for iter = 1:20
48 kc = kb .* actb + kt .* actt;
49 Alin = A0 + spdiags(kc, 0, N, N);
50 rhs = fg + kb .* actb .* ob - kt .* actt .* ot;
51 w0 = zeros(N, 1);
52 w0(iu) = plate_static_solve(Alin(iu, iu), rhs(iu));
53 nb = (kb > 0) & (ob - w0 > 0);
54 nt = (kt > 0) & (ot + w0 > 0);
55 if isequal(nb, actb) && isequal(nt, actt)
56 break
57 end
58 actb = nb; actt = nt;
59 if iter == 20
60 notes{end + 1} = 'the set of engaged contacts did not settle in 20 iterations; the last one is used';
61 end
62 end
63 if any(kb > 0) && ~all(actb(kb > 0))
64 notes{end + 1} = sprintf('%d support node(s) carry no load at rest: the plate does not touch them', nnz(kb > 0 & ~actb));
65 end
66 if any(kt > 0) && ~all(actt(kt > 0))
67 notes{end + 1} = sprintf('%d upper-felt node(s) are clear of the plate at rest', nnz(kt > 0 & ~actt));
68 end
69end
70st.w0 = w0;
71st.actb = actb; st.actt = actt;
72st.gb = ob - w0; st.gt = ot + w0;
73st.fstat_s = -ks .* w0; % linear springs
74st.fstat_b = kb .* max(st.gb, 0); % bottom contacts, pushing up
75st.fstat_t = -kt .* max(st.gt, 0); % top contacts, pushing down
76st.Alin = Alin;
77st.sag = max(abs(w0));
78st.notes = notes;
79end
81function x = plate_static_solve(A, b)
82[Rc, p] = chol(A);
83if p > 0
84 error('plate_static: the plate has no stable equilibrium in this mount (it tips, swings over or buckles under its own weight); check the pads, the holes or the tilt');
85end
86x = Rc \ (Rc' \ b);
87end