80c0b6dPlate dynamics lab: plates in realizable mounts under gravity, explicit finite differences on WebGPUJeremy Magland 1function mdl = plate_model(params)
2% PLATE_MODEL The discrete plate in its mount: operators, statics and time step.
3%
4% mdl = plate_model(params) assembles everything about the problem that
5% does not depend on the excitation: the grid (plate_domain), the mount
6% (plate_mount), the bending and Laplacian operators (plate_bending), the
7% self-weight prestress (plate_prestress), the rigid-body basis
8% (plate_rigid_basis), the static equilibrium (plate_static), the reduced
9% rigid-body system and the stable time step. plate_fd_plan reports its
10% numbers before a run and plate_fd_setup builds the run on it.
11%
12% The semi-discrete equation on the unknown nodes is
13%
14% M w_tt + 2 sigma0 M w_t - 2 sigma1 mu Le w_t + A0 w + f_contact(w, w_t) = F(t) s + M g_n
15% A0 = Kb + Kg + Ks
16%
17% with M the node masses, Kb the bending stiffness, Kg the prestress
18% operator, Ks the linear springs of the mount (threads, fingers), Le the
19% Laplacian energy operator, and f_contact the one-sided pads and felts
20% with their dashpots. Writing w = w0 + R c + phi, with w0 the static
21% deflection, R the rigid-body basis and c its coordinates, and phi the
22% flexural remainder (mass-orthogonal to R), the equation splits into
23%
24% G (c_tt + 2 sigma0 c_t) = R' f - Krr c - AR' phi, G = R' M R,
25% phi_tt + 2 sigma0 phi_t = M^-1 (f - A0 phi - AR c) - R (c_tt + 2 sigma0 c_t),
26%
27% where f collects the applied and contact forces about rest and the sigma1
28% loss 2 sigma1 mu Le w_t of the total velocity (so that it enters both
29% parts and stays dissipative), AR = A0 R is the force of a unit rigid
30% displacement (with the thread's torsion added as a force on the spin) and
31% Krr = R' AR the reduced stiffness of the rigid motions (the pendulum
32% stiffness, the thread springs, the torsion). Since Kb R = 0, AR is
33% computed from Kg + Ks alone; both are exact to rounding. The sum
34% R c + phi then obeys exactly the leapfrog of the whole system, and the
35% rigid content of phi, Q' phi, stays zero: every force acts on the whole
36% plate, and the rigid equations are its projection. (Leaving AR c out of
37% the flexural equation, an early version of this code, made the two parts
38% pump each other and the swing of a hung plate grow without bound.)
39%
40% Stable time step. For the bending operator alone the von Neumann bound of
41% the explicit scheme is dt <= 2 / (lambda (sigma1 + sqrt(kappa^2 + sigma1^2)))
42% with lambda = 4/dx^2 + 4/dy^2 the largest eigenvalue of the five-point
43% Laplacian. The extra stiffness of the prestress, the springs and the
44% engaged contacts adds a rate Omega^2 bounded by Gershgorin's theorem on
45% M^-1 (Kg + Ks + K_contact), and the bound becomes
46%
47% dt <= ( -2 sigma1 lambda + sqrt(4 sigma1^2 lambda^2 + 4 (kappa^2 lambda^2 + Omega^2)) )
48% / (kappa^2 lambda^2 + Omega^2),
49%
50% which reduces to the first form when Omega = 0. The dashpots are treated
51% explicitly and need dt < 2 M / c; the plan keeps dt below half of that.
52% The largest eigenvalue of the energy-method bending operator lies just
53% below the interior bound (about 0.97 to 0.99 of it), so the step is safe
54% and slightly conservative.
55%
56% Fields: P, dom, mnt, unk, Mvec, Kb, Le, Lm, Kg, Ks, A0, R, Q, G, Ginv,
57% Krr, AR, ps (prestress solution), st (statics), rigid_hz (frequencies of
58% the rigid-body modes, ascending, NaN for an unstable one) and
59% rigid_mode_names (the coordinate that dominates each), lam_max, Om2, dt_max, dt,
60% nsteps, ncells, nunknowns, notes.
62P = plate_setup(params);
63nx = P.fd.nx; ny = P.fd.ny;
64if nx < 5 || ny < 5
65 error('plate_model: the grid needs at least 5 nodes in each direction');
66end
67dom = plate_domain(P.plate, nx, ny);
68dx = dom.dx; dy = dom.dy;
69N = nx * ny;
70mnt = plate_mount(P, dom);
71unk = dom.inside & ~mnt.held;
72a = dom.m * dx * dy;
73Mvec = P.mu * a(:);
74Mvec(~unk(:)) = 1; % dummy mass at nodes outside the system
75[Kb, Le, ops] = plate_bending(dom, mnt.held, P.material.nu, P.D);
76[Kg, ps] = plate_prestress(P, dom, mnt, ops);
77ksv = mnt.ks(:) .* double(unk(:));
78Ks = spdiags(ksv, 0, N, N);
79A0 = Kb + Kg + Ks;
80Lm = spdiags(1 ./ max(a(:), eps) .* double(unk(:)), 0, N, N) * Le;
81[R, Q, G] = plate_rigid_basis(dom, mnt, Mvec, unk);
82st = plate_static(P, dom, mnt, A0, Mvec, unk);
83nr = mnt.nr;
84notes = [mnt.notes, st.notes];
86% Reduced rigid-body system: bending annihilates R, so only Kg and Ks act on it
87if nr > 0
88 AR = (Kg + Ks) * R;
89 if mnt.ktor > 0
90 % the thread's torsion as a force on the plate: -ktor c3 M r_h / (r_h' M r_h),
91 % which changes only the spin coordinate, so that Krr = R' AR stays exact
92 AR(:, 3) = AR(:, 3) + mnt.ktor * (Mvec .* R(:, 3)) / G(3, 3);
93 end
94 Krr = R' * AR;
95 Krr = (Krr + Krr') / 2;
96 Ginv = inv(G);
97 % modes of the reduced system in the mass-normalized coordinates, so that
98 % the dominant coordinate of each is meaningful across metres and radians
99 Rg = chol(G); % G = Rg' Rg
100 Bm = (Rg' \ Krr) / Rg;
101 [V, Dm] = eig((Bm + Bm') / 2);
102 lam = real(diag(Dm));
103 [lam, order] = sort(lam);
104 V = real(V(:, order));
105 rigid_hz = sqrt(max(lam, 0)) / (2 * pi);
106 rigid_hz(lam < -1e-9 * max(abs(lam))) = NaN;
107 rigid_mode_names = cell(1, nr);
108 for k = 1:nr
109 [~, kd] = max(abs(V(:, k)));
110 rigid_mode_names{k} = mnt.rigid_names{kd};
111 end
112 if any(isnan(rigid_hz))
113 notes{end + 1} = 'a rigid-body motion has negative stiffness: the plate is in an unstable position and the run will diverge';
114 end
115else
116 AR = zeros(N, 0); Krr = zeros(0, 0); Ginv = zeros(0, 0); rigid_hz = zeros(0, 1); rigid_mode_names = {};
117end
119% Stable time step
120lam_max = 4 / dx^2 + 4 / dy^2;
121kc = mnt.kb(:) .* double(st.actb) + mnt.kt(:) .* double(st.actt);
122Kx = abs(Kg) + spdiags(abs(ksv) + kc, 0, N, N);
123row = sum(Kx, 2);
124Om2 = max(row(unk(:)) ./ Mvec(unk(:)));
125if isempty(Om2)
126 Om2 = 0;
127end
128s1 = P.damping.sigma1;
129om2b = P.kappa^2 * lam_max^2 + Om2;
130dt_max = (-2 * s1 * lam_max + sqrt(4 * s1^2 * lam_max^2 + 4 * om2b)) / om2b;
131cv = (mnt.cb(:) + mnt.ct(:) + mnt.cs(:)) ./ Mvec;
132cmax = max(cv(unk(:)));
133if cmax > 0 && dt_max > 0.5 / cmax
134 dt_max = 0.5 / cmax;
135 notes{end + 1} = 'the support damping, not the bending, limits the time step';
136end
137dt = P.fd.cfl * dt_max;
138nsteps = ceil(P.fd.duration / dt);
140mdl.P = P; mdl.dom = dom; mdl.mnt = mnt; mdl.unk = unk;
141mdl.nx = nx; mdl.ny = ny; mdl.dx = dx; mdl.dy = dy; mdl.N = N;
142mdl.a = a; mdl.Mvec = Mvec;
143mdl.Kb = Kb; mdl.Le = Le; mdl.Lm = Lm; mdl.Kg = Kg; mdl.Ks = Ks; mdl.A0 = A0;
144mdl.R = R; mdl.Q = Q; mdl.G = G; mdl.Ginv = Ginv; mdl.Krr = Krr; mdl.AR = AR; mdl.nr = nr;
145mdl.ps = ps; mdl.st = st; mdl.ops = ops;
146mdl.rigid_hz = rigid_hz; mdl.rigid_mode_names = rigid_mode_names;
147mdl.lam_max = lam_max; mdl.Om2 = Om2;
148mdl.dt_max = dt_max; mdl.dt = dt; mdl.nsteps = nsteps;
149mdl.ncells = N; mdl.nunknowns = nnz(unk);
150mdl.notes = notes;
151end