1function S = plate_stress(P, wxx, wyy, wxy)
2% PLATE_STRESS Surface stresses of a Kirchhoff-Love plate from its curvatures.
3%
4% S = plate_stress(P, wxx, wyy, wxy) takes the second derivatives of the
5% displacement (arrays of any common size) and returns, at the surface z = h/2,
6% S.sxx, S.syy, S.txy in-plane stress components [Pa]
7% S.vm von Mises equivalent stress [Pa]
8% S.sp principal stress of largest magnitude, signed
9% (tension positive) [Pa]
10%
11% In thin-plate theory the in-plane strains vary linearly through the
12% thickness, eps_xx = -z w_xx, eps_yy = -z w_yy, gamma_xy = -2 z w_xy, and
13% plane stress gives
14% sigma_xx = -E z / (1 - nu^2) (w_xx + nu w_yy)
15% sigma_yy = -E z / (1 - nu^2) (w_yy + nu w_xx)
16% tau_xy = -E z / (1 + nu) w_xy
17% Stresses are largest at the surfaces z = +-h/2 and vanish at the mid-plane;
18% the two surfaces carry equal and opposite stress. P.Kb and P.Ks are the
19% factors E h / (2 (1 - nu^2)) and E h / (2 (1 + nu)) (see plate_setup).
20%
21% Because stress is a second derivative of the primary unknown, any error in
22% the computed displacement is amplified here.
24nu = P.material.nu;
25S.sxx = -P.Kb * (wxx + nu * wyy);
26S.syy = -P.Kb * (wyy + nu * wxx);
27S.txy = -P.Ks * wxy;
29S.vm = sqrt(S.sxx.^2 - S.sxx .* S.syy + S.syy.^2 + 3 * S.txy.^2);
31c = (S.sxx + S.syy) / 2;
32r = sqrt(((S.sxx - S.syy) / 2).^2 + S.txy.^2);
33s1 = c + r;
34s2 = c - r;
35S.sp = s1;
36pick = abs(s2) > abs(s1);
37S.sp(pick) = s2(pick);
38end