% fig_08_31b
% circumferential stress component at edge of cylindrical
% hole with internal pressure in orthotropic elastic body
% compare to isotropic case
% Jaeger and Cook (1979)
% equation (8.146)

clear all, clf reset; % clear memory and figures
p = 1; % magnitude of internal pressure (MPa)
TH = 0:pi/180:pi; THD = TH*180/pi; % Angle theta
C2T = cos(2*TH); C4T = cos(4*TH);
STTI = p*ones(size(TH)); % isotropic case

E1 = 40000; E2 = 20000; % Young's moduli (MPa)
G = 10000; % Shear modulus (MPa)
nu12 = 0.2; nu21 = nu12*E2/E1; % Poisson's ratios
c1 = E2/E1; c2 = (E2/G) - 2*nu21;
cond = c2^2 - 4*c1 % NOTE: cond must be zero or positive
a1 = 0.5*c2 + 0.5*sqrt(cond); a2 = c1/a1;
g1 = (a1^0.5 -1)/(a1^0.5 +1); g2 = (a2^0.5 -1)/(a2^0.5 +1);
NUM = 1+(g1+g2)^2-3*g1^2*g2^2-2*(g1+g2)*(1-g1*g2)*C2T-2*g1*g2*C4T;
DEN = (1+g1^2-2*g1*C2T).*(1+g2^2-2*g2*C2T);
STTA = p*NUM./DEN; % anisotropic case

plot(THD,STTI,'k-',THD,STTA,'k-.'); 
xlabel('theta (degrees)'); ylabel('stress (MPa)');
legend('STTI','STTA'); axis ([0 180 0.75 1.25]);