%Mason Averill
%ME-550 Fall 2020, 10/27
%Design Problem 1


%Input Parameters

%Equivalent Loads, Torques, and Moments Applied
P_axial=40;%in terms of p
T=500;%in terms of p
M_y=-1100;%in terms of p
M_z=100;%in terms of p


%Material Properties

Yield_strength=36300;%Based on A-36 steel
v=0.26;%Based on A-36 steel

d=0.1;%lower bound for d
delta_d=0.5;%increment in d per iteration

d_upper=20;%upper bound for d

size=round((d_upper-d)/delta_d+1);


p_store_tensile_max_normal_stress=zeros(1,size);

p_store_tensile_max_shear_stress=zeros(1,size);
p_store_compression_max_shear_stress=zeros(1,size);


p_store_tensile_max_normal_strain=zeros(1,size);
p_store_compression_max_normal_strain=zeros(1,size);

p_store_tensile_octahedral=zeros(1,size);
p_store_compression_octahedral=zeros(1,size);

p_store_tensile_DET=zeros(1,size);
p_store_compression_DET=zeros(1,size);


d_store=zeros(1,size);

i=1;

while(d<d_upper+delta_d)

width=2*d;
height=d;
thickness=0.1*d;


%Cross-Section Properties
I_zz=2*(1/12*thickness*(height+thickness)^3)+2*(1/12*(width-thickness)*(thickness)^3+(width-thickness)*(thickness)*(height/2)^2);

I_yy=2*(1/12*(thickness)*(width+thickness)^3)+2*(1/12*(height-thickness)*(thickness)^3+(height-thickness)*(thickness)*(width/2)^2);

I_yz=0;

Area_mean=width*height;

Area=(width+thickness)*(height+thickness)-(width-thickness)*(height-thickness);

%Location of Neutral Axis

argument=((M_y*I_zz-M_z*I_yz)/(M_z*I_yy-M_y*I_yz));

theta=double(atan(-argument));

theta_degrees=theta*180/pi();

%Locate Points of Max Compression and Tension


z_tension_max=-((width/2)+(thickness/2));
y_tension_max=((height/2)+(thickness/2));

z_compression_max=((width/2)+(thickness/2));
y_compression_max=-((height/2)+(thickness/2));


%Compute Stresses

%Axial Stress

Stress_axial=P_axial/Area;


%Shear Stress (same at both points of interest)

Stress_shear=T/(2*Area_mean*(thickness));

%Bending Stress

    %Max Tensile Stress

    Stress_bending_tension=(M_z*y_tension_max)/I_zz+(M_y*z_tension_max)/I_yy;
    
    %Max Compressive Stress

    Stress_bending_compression=(M_z*y_compression_max)/I_zz+(M_y*z_compression_max)/I_yy;


%Sum Stresses at Both Locations of Interest

    %Stress at Max Tensile Stress Location
    
    Stress_x_tension=Stress_axial+Stress_bending_tension;
    
    %Stress at Max Compression Stress Location

    Stress_x_compression=Stress_axial+Stress_bending_compression;
    
    
  %Now lets find principal stresses  
    
  principal_stress_max_tensile_1=(Stress_x_tension/2)+sqrt((Stress_x_tension/2)^2+(Stress_shear)^2); %Max 
  
  principal_stress_max_tensile_2=(Stress_x_tension/2)-sqrt((Stress_x_tension/2)^2+(Stress_shear)^2); %Min (and less than 0)
  
  
  principal_stress_max_compression_1=(Stress_x_compression/2)+sqrt((Stress_x_compression/2)^2+(Stress_shear)^2); %Max 
  
  principal_stress_max_compression_2=(Stress_x_compression/2)-sqrt((Stress_x_compression/2)^2+(Stress_shear)^2); %Min (and less than 0)
  
  %Solve for P based on Maximum Normal Stress failure criterion (for brittle
%materials, i.e. tension controls)

        %At Max Tensile Location
        p_solved_tensile_max_normal_stress=Yield_strength/principal_stress_max_tensile_1;


%Solve for P based on Maximum Shear Stress failure criterion 


        %At Max Tensile Location
        
        p_solved_tensile_max_shear_stress=Yield_strength/(principal_stress_max_tensile_1-principal_stress_max_tensile_2);
        
        %At Max Compression Location


        p_solved_compression_max_shear_stress=Yield_strength/(principal_stress_max_compression_1-principal_stress_max_compression_2);

%Solve for P based on Maximum Normal Strain failure criterion 


        %At Max Tensile Location
        
        p_solved_tensile_max_normal_strain=Yield_strength/(principal_stress_max_tensile_1-v*principal_stress_max_tensile_2);
        
        %At Max Compression Location


        p_solved_compression_max_normal_strain=Yield_strength/(principal_stress_max_compression_1-v*principal_stress_max_compression_2);

        
        %Solve for P based on Maximum Distortion Energy failure criterion 


        %At Max Tensile Location
        
        p_solved_tensile_DET=sqrt(2*(Yield_strength)^2/((principal_stress_max_tensile_1-principal_stress_max_tensile_2)^2+principal_stress_max_tensile_1^2+principal_stress_max_tensile_2^2));
        
        %At Max Compression Location


        p_solved_compression_DET=sqrt(2*(Yield_strength)^2/((principal_stress_max_compression_1-principal_stress_max_compression_2)^2+principal_stress_max_compression_1^2+principal_stress_max_compression_2^2));

        
        
        %Solve for P based on Octahedral Shear Stress failure criterion 


        %At Max Tensile Location
        
        p_solved_tensile_octahedral=sqrt(2)*Yield_strength/sqrt(((principal_stress_max_tensile_1-principal_stress_max_tensile_2)^2+principal_stress_max_tensile_1^2+principal_stress_max_tensile_2^2));
        
        %At Max Compression Location


        p_solved_compression_octahedral=sqrt(2)*Yield_strength/sqrt(((principal_stress_max_compression_1-principal_stress_max_compression_2)^2+principal_stress_max_compression_1^2+principal_stress_max_compression_2^2));
        
        
        
%Store P's and d

    %Store d
    
    d_store(1,i)=d;
    
    %Store P's

    p_store_tensile_max_normal_stress(1,i)=p_solved_tensile_max_normal_stress;
    
    p_store_tensile_max_shear_stress(1,i)=p_solved_tensile_max_shear_stress;
    p_store_compression_max_shear_stress(1,i)=p_solved_compression_max_shear_stress;
    
    
    p_store_tensile_max_normal_strain(1,i)=p_solved_tensile_max_normal_strain;
    p_store_compression_max_normal_strain(1,i)=p_solved_compression_max_normal_strain;
    
    p_store_tensile_DET(1,i)=p_solved_tensile_DET;
    p_store_compression_DET(1,i)=p_solved_compression_DET;
    
    p_store_tensile_octahedral(1,i)=p_solved_tensile_octahedral;
    p_store_compression_octahedral(1,i)=p_solved_compression_octahedral;

  
    
%increment d and i

d=d+delta_d;
i=i+1;

end 
  
  
figure (1)
plot(d_store,p_store_tensile_max_normal_stress,'r')
xlabel('d (in)')
ylabel('Magnitude of P that will cause yielding (lbf)')
title('Magnitude of P that will cause yielding vs d: Max Tensile Location');
grid on

hold on
plot(d_store,p_store_tensile_max_shear_stress,'b')
plot(d_store,p_store_tensile_max_normal_strain,'g')
plot(d_store,p_store_tensile_DET,'k--o')
plot(d_store,p_store_tensile_octahedral,'m')
legend('Maximum Normal Stress Theory','Maximum Shear Stress Theory','Maximum Normal Strain Theory','Distortion Energy Theory','Octahedral Theory')
xlim([0 d_upper])

hold off

figure (2)
plot(d_store,p_store_compression_max_shear_stress,'b')
xlabel('d (in)')
ylabel('Magnitude of P that will cause yielding (lbf)')
title('Magnitude of P that will cause yielding vs d: Max Compression Location');
grid on

hold on
%plot(d_store,p_store_compression_max_normal_strain,'g')
plot(d_store,p_store_compression_DET,'k--o')
plot(d_store,p_store_compression_octahedral,'m')
legend('Maximum Shear Stress Theory','Distortion Energy Theory','Octahedral Theory')
xlim([0 d_upper])
hold off