1function P = plate_setup(params)
2% PLATE_SETUP Derived quantities for a Kirchhoff-Love plate problem.
3%
4% P = plate_setup(params) returns a copy of params with the derived fields
5% P.D bending stiffness D = E h^3 / (12 (1 - nu^2)) [N m]
6% P.mu mass per unit area mu = rho h [kg/m^2]
7% P.kappa stiffness parameter kappa = sqrt(D / mu) [m^2/s]
8% P.Kb bending-stress factor at the surface z = h/2 [N/m]
9% P.Ks shear-stress factor at the surface z = h/2 [N/m]
10%
11% The plate is params.plate.shape, 'rectangle' or 'ellipse', filling the
12% box [0, Lx] x [0, Ly]: for the ellipse Lx and Ly are the diameters along
13% x and y, and Lx = Ly is a disk. A missing shape means a rectangle.
14%
15% The transverse displacement w(x,y,t) of the plate obeys
16%
17% w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = f(x,y,t) / mu
18%
19% where f is the applied pressure [N/m^2], sigma0 [1/s] is a frequency-independent
20% loss rate and sigma1 [m^2/s] adds losses that grow with spatial frequency
21% (Bilbao, Numerical Sound Synthesis, 2009, ch. 12). The surface stresses follow
22% from the curvatures (see plate_stress).
24P = params;
25pl = params.plate;
26mat = params.material;
28if ~isfield(pl, 'shape')
29 pl.shape = 'rectangle';
30 P.plate.shape = 'rectangle';
31end
32if ~any(strcmp(pl.shape, {'rectangle', 'ellipse'}))
33 error('plate_setup: unknown plate shape "%s" (expected rectangle or ellipse)', pl.shape);
34end
35if pl.Lx <= 0 || pl.Ly <= 0 || pl.h <= 0
36 error('plate_setup: plate dimensions must be positive');
37end
38if pl.h > 0.1 * min(pl.Lx, pl.Ly)
39 warning('plate_setup: thickness is more than a tenth of the shorter side; Kirchhoff-Love thin-plate theory is not reliable here');
40end
41if ~any(strcmp(params.bc, {'simply-supported', 'clamped', 'free', 'pinned-center'}))
42 error('plate_setup: unknown boundary condition "%s" (expected simply-supported, clamped, free or pinned-center)', params.bc);
43end
45P.D = mat.E * pl.h^3 / (12 * (1 - mat.nu^2));
46P.mu = mat.rho * pl.h;
47P.kappa = sqrt(P.D / P.mu);
48P.Kb = mat.E * pl.h / (2 * (1 - mat.nu^2));
49P.Ks = mat.E * pl.h / (2 * (1 + mat.nu));
50end