1function S = plate_stress(P, wxx, wyy, wxy)
2% PLATE_STRESS Surface stresses of a plate from its curvatures.
3%
4% S = plate_stress(P, wxx, wyy, wxy) takes the second derivatives of the
5% transverse displacement (arrays of any common size) and returns, at the
6% surface z = h/2,
7% S.sxx, S.syy, S.txy in-plane stress components [Pa]
8% S.vm von Mises equivalent stress [Pa]
9% S.sp principal stress of largest magnitude, signed
10% (tension positive) [Pa]
11%
12% The in-plane strains vary linearly through the thickness,
13% eps_xx = -z w_xx, eps_yy = -z w_yy, gamma_xy = -2 z w_xy, and plane
14% stress gives
15% sigma_xx = -E z / (1 - nu^2) (w_xx + nu w_yy)
16% sigma_yy = -E z / (1 - nu^2) (w_yy + nu w_xx)
17% tau_xy = -E z / (1 + nu) w_xy
18% Stresses are largest at the surfaces z = +-h/2 and vanish at the
19% mid-plane; the two surfaces carry equal and opposite stress. P.Kb and
20% P.Ks are the factors E h / (2 (1 - nu^2)) and E h / (2 (1 + nu)).
21%
22% For the Mindlin-Reissner method the curvatures are the gradients of the
23% rotations, -d(theta_x)/dx and so on; the same formula then gives the
24% bending stress, without the transverse-shear stress, which is small at
25% these thicknesses.
27nu = P.material.nu;
28S.sxx = -P.Kb * (wxx + nu * wyy);
29S.syy = -P.Kb * (wyy + nu * wxx);
30S.txy = -P.Ks * wxy;
32S.vm = sqrt(S.sxx.^2 - S.sxx .* S.syy + S.syy.^2 + 3 * S.txy.^2);
34c = (S.sxx + S.syy) / 2;
35r = sqrt(((S.sxx - S.syy) / 2).^2 + S.txy.^2);
36s1 = c + r;
37s2 = c - r;
38S.sp = s1;
39pick = abs(s2) > abs(s1);
40S.sp(pick) = s2(pick);
41end