78d40e9Free edges for both finite-difference methodsJeremy Magland 1function [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu)
2% PLATE_MATRICES Sparse Laplacian, stiffness and mass matrices for the plate.
78d40e9Free edges for both finite-difference methodsJeremy Magland 4% [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu) assembles
5% the spatial operators of the finite-difference schemes as sparse matrices
6% acting on the vector w(idx) of unknowns, where idx indexes the ny x nx field
7% with y running fastest (column-major, as w(:) does) and the unknowns form an
8% my x mx array. The semi-discrete equation is
78d40e9Free edges for both finite-difference methodsJeremy Magland 10% M w_tt + 2 sigma0 M w_t - 2 sigma1 L w_t + kappa^2 K w = F(t) M s / mu
78d40e9Free edges for both finite-difference methodsJeremy Magland 12% Simply supported and clamped edges (w = 0 on the edge nodes). The unknowns
13% are the mx = nx-2 by my = ny-2 interior nodes, M is the identity, L the
14% five-point Laplacian and K the thirteen-point biharmonic, with the edges
15% handled exactly as in the ghost-node scheme of plate_fd / plate_pad: D2*D2
16% encodes the odd reflection (w_nn = 0), and for clamped edges the even
17% reflection changes the end coefficient of the fourth difference from 5/dx^4
18% to 7/dx^4. The mixed term 2 w_xxyy never reaches a ghost node. nu is not
19% used. A test in the lab checks K against the explicit stencil.
20%
21% Free edges. Every node is an unknown (mx = nx, my = ny). The free-edge
22% conditions, no bending moment and no Kirchhoff shear on an edge and no
23% twisting moment at a corner, involve third derivatives and do not fit in a
24% ghost layer. K is therefore assembled from the bending energy
25%
26% U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
27%
28% using second differences at the nodes, the cross difference at the cell
29% centres and trapezoidal quadrature weights (the finite-difference energy
30% method; Bushnell, Computers & Structures 1973). The free-edge conditions
31% are then natural boundary conditions of the discrete energy: the scheme is
32% second order inside and first order at the boundary. M holds the trapezoidal
33% node weights (1 inside, 1/2 on an edge, 1/4 at a corner) and L is the
34% Laplacian energy operator with the natural zero-normal-derivative condition.
35% K is symmetric positive semidefinite and its null space is exactly the three
36% rigid-body motions 1, x and y, which plate_rigid_modes removes from the
37% state. All three matrices are divided by dx*dy, so inside the plate M = 1
38% and K and L reduce to the usual stencils, and the largest eigenvalue of
39% M \ K stays below the interior bound (4/dx^2 + 4/dy^2)^2 that plate_fd_plan
40% uses for the time step.
41%
42% Everything is built from Kronecker products so the matrices stay sparse;
43% kron(Ax, By) acts on w(:) as By * W * Ax'.
46 nu = NaN;
47end
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 48switch bc
78d40e9Free edges for both finite-difference methodsJeremy Magland 49 case {'simply-supported', 'clamped'}
50 mx = nx - 2;
51 my = ny - 2;
52 ex = ones(mx, 1);
53 ey = ones(my, 1);
54 D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
55 D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
56 Ix = speye(mx);
57 Iy = speye(my);
58 D4x = D2x * D2x;
59 D4y = D2y * D2y;
60 if strcmp(bc, 'clamped')
61 D4x(1, 1) = 7 / dx^4; D4x(mx, mx) = 7 / dx^4;
62 D4y(1, 1) = 7 / dy^4; D4y(my, my) = 7 / dy^4;
63 end
64 L = kron(Ix, D2y) + kron(D2x, Iy);
65 K = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
66 M = speye(mx * my);
67 mask = false(ny, nx);
68 mask(2:end - 1, 2:end - 1) = true;
69 idx = find(mask);
70 case 'free'
71 if ~isfinite(nu)
72 error('plate_matrices: free edges need Poisson''s ratio nu');
73 end
74 mx = nx;
75 my = ny;
76 % One-dimensional pieces, as triplets so the layout is unambiguous:
77 % second differences at the interior nodes ((n-2) x n), first differences
78 % at the cell midpoints ((n-1) x n), selectors of the interior nodes, and
79 % trapezoidal node weights.
80 kx = (1:nx - 2)'; ky = (1:ny - 2)';
81 d2x = sparse([kx; kx; kx], [kx; kx + 1; kx + 2], [ones(nx - 2, 1); -2 * ones(nx - 2, 1); ones(nx - 2, 1)], nx - 2, nx) / dx^2;
82 d2y = sparse([ky; ky; ky], [ky; ky + 1; ky + 2], [ones(ny - 2, 1); -2 * ones(ny - 2, 1); ones(ny - 2, 1)], ny - 2, ny) / dy^2;
83 jx = (1:nx - 1)'; jy = (1:ny - 1)';
84 d1x = sparse([jx; jx], [jx; jx + 1], [-ones(nx - 1, 1); ones(nx - 1, 1)], nx - 1, nx) / dx;
85 d1y = sparse([jy; jy], [jy; jy + 1], [-ones(ny - 1, 1); ones(ny - 1, 1)], ny - 1, ny) / dy;
86 Sx = sparse(kx, kx + 1, ones(nx - 2, 1), nx - 2, nx);
87 Sy = sparse(ky, ky + 1, ones(ny - 2, 1), ny - 2, ny);
88 tx = ones(nx, 1); tx([1 end]) = 0.5; Tx = spdiags(tx, 0, nx, nx);
89 ty = ones(ny, 1); ty([1 end]) = 0.5; Ty = spdiags(ty, 0, ny, ny);
90 % The four terms of the energy density, each as D' * (weights) * D
91 Kxx = kron(d2x' * d2x, Ty); % w_xx^2 at nodes 2..nx-1, every row
92 Kyy = kron(Tx, d2y' * d2y); % w_yy^2 at nodes 2..ny-1, every column
93 Kc = nu * (kron(d2x' * Sx, Sy' * d2y) + kron(Sx' * d2x, d2y' * Sy)); % 2 nu w_xx w_yy where both exist
94 Kxy = 2 * (1 - nu) * kron(d1x' * d1x, d1y' * d1y); % 2 (1-nu) w_xy^2 at the cell centres
95 K = Kxx + Kyy + Kc + Kxy;
96 L = -(kron(d1x' * d1x, Ty) + kron(Tx, d1y' * d1y));
97 M = kron(Tx, Ty);
98 idx = (1:nx * ny)';
2546961Plate mechanics lab: Kirchhoff-Love plates solved by numbl in the browserJeremy Magland 99 otherwise
100 error('plate_matrices: unknown boundary condition "%s"', bc);
101end
102end