ME3227_FinalProject.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ME 3227 Final Project
% This code allows the user to automatically calculate and visualize the
% required diameters if the design requirements change. The requirements that can be changed are the gear weight, the number of teeth, the pressure angle, and the diametric pitch.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Given Parameters: 
% Shaft Material - AISI 4130 Q&T CD
E = 29*10^6; %psi (Young's Modulus)
sut = 118*10^3; %psi (ultimate strength)
sy = 102*10^3; %psi (yield strength)
p = 0.3; %(poisson's ratio)
ga_w = 4; %lbf - make as input value
ga_t = 50; %teeth number - input
ga_a = 20; %deg - input
ga_d = 4; %diametric pitch - input
gc_w = 2; %lbf - make as input value
gc_t = 25; %teeth number - input
gc_a = 20; %deg - input
gc_d = 4; %diametric pitch - input
ss = 3000; %rpm (shaft speed)
T = 3500; %lbf-in (torque)
% Reliability for fatigue life = 0.99
% nd, Safety factor for fatigue life = 2
% nd, safety factor for all other criteria = 1
% r/d = 0.02 - woodruff key
gear_def = 0.01; % deflection at gears - from table 7-2
bearing_slope = 0.001; %rad
gear_slope = 0.005; %rad 
d_a = ga_t/ga_d; % diameter of gear a
d_c = gc_t/gc_d; % diameter of gear c
Fct = T/(d_c/2); % tangential force on gear c
Fat = T/(d_a/2); % tangential force on gear a
Fcn = Fct * tand(gc_a); % Normal force on gear c
Fan = Fat * tand(ga_a); % Normal force on gear a
dtheta_dx = 0.00058; % rad/in

% Find Reaction Forces - XY plane
Rby = ((Fcn*3)+(Fan*9))/6;
Rdy = Fan + Fcn - Rby; 

% Find Reaction Forces - XZ plane
Rbz = ((Fat*9)-(Fct*3))/6;
Rdz = Fat - Fct - Rbz; 

% Calculate C1 and C2 - XY Plane
c1y = (117*Fan-36*Rby+4.5*Fcn)/6;
c2y = 4.5*Fan - 3*c1y; 

% Calculate C1 and C2 - XZ Plane
c1z = (117*Fat-36*Rbz-4.5*Fct)/6;
c2z = 4.5*Fat - 3*c1z;

% Singularity Functions - x = 0 - 9 inches
x = linspace(0,9,10);
for i=0:9
   if i<(3)
       q_xy(i+1) = 0;
       M_xy(i+1) = 0;
       EIdy_x(i+1) = -Fan/2*(i)^2+c1y;
       EIy_x(i+1) = -Fan/6*(i)^3+c1y*i+c2y;
       M_xz(i+1) = 0;
       EIdz_x(i+1) = -Fat/2*(i)^2+c1z;
       EIz_x(i+1) = -Fat/6*(i)^3+c1z*i+c2z;
   elseif i>=3 && i<6
       q_xy(i+1) = -Fan*(i)^-1 + Rby*(i-3)^(-1);
       M_xy(i+1) = -Fan*i + Rby*(i-3);
       EIdy_x(i+1) = -Fan/2*(i)^2+Rby/2*(i-3)^2+c1y;
       EIy_x(i+1) = -Fan/6*(i)^3+Rby/6*(i-3)^3+c1y*i+c2y;
       M_xz(i+1) = -Fat*i + Rbz*(i-3);
       EIdz_x(i+1) = -Fat/2*(i)^2+Rbz/2*(i-3)^2+c1z;
       EIz_x(i+1) = -Fat/6*(i)^3+Rbz/6*(i-3)^3+c1z*i+c2z;
   elseif i>=6 && i<9 
       q_xy(i+1)  = -Fan*i^-1+Rby*(i-3)^-1-Fcn*(i-6)^-1;
       M_xy(i+1)  = -Fan*i+Rby*(i-3)-Fcn*(i-6);
       EIdy_x(i+1) = -Fan/2*(i)^2+Rby/2*(i-3)^2-Fcn/2*(i-6)^2+c1y;
       EIy_x(i+1) = -Fan/6*(i)^3+Rby/6*(i-3)^3-Fcn/6*(i-6)^3+c1y*i+c2y;
       M_xz(i+1)  = -Fat*i+Rbz*(i-3)+Fct*(i-6);
       EIdz_x(i+1) = -Fat/2*(i)^2+Rbz/2*(i-3)^2+Fct/2*(i-6)^2+c1z;
       EIz_x(i+1) = -Fat/6*(i)^3+Rbz/6*(i-3)^3+Fct/6*(i-6)^3+c1z*i+c2z;
   else
       q_xy(i+1) = -Fan*i^-1 +Rby*(i-3)^-1-Fcn*(i-6)^-1+Rdy*(i-9)^-1;
       M_xy(i+1) = -Fan*i +Rby*(i-3)-Fcn*(i-6)+Rdy*(i-9);
       EIdy_x(i+1) = -Fan/2*(i)^2+Rby/2*(i-3)^2-Fcn/2*(i-6)^2+Rdy/2*(i-9)^2+c1y;
       EIy_x(i+1) = -Fan/6*(i)^3+Rby/6*(i-3)^3-Fcn/6*(i-6)^3+Rdy/2*(i-9)^3+c1y*i+c2y;
       M_xz(i+1) = -Fat*i +Rbz*(i-3)+Fct*(i-6)+Rdz*(i-9);
       EIdz_x(i+1) = -Fat/2*(i)^2+Rbz/2*(i-3)^2+Fct/2*(i-6)^2+Rdz/2*(i-9)^2+c1z;
       EIz_x(i+1) = -Fat/6*(i)^3+Rbz/6*(i-3)^3+Fct/6*(i-6)^3+Rdz/2*(i-9)^3+c1z*i+c2z;
   end
   i = i+1; 
end

%Create plots
 figure(1) %position vs. deflection in xy plane
 plot(x,EIy_x)
 xlabel('Axial Postion')
 ylabel('EI times Deflection')
 figure(2) %position vs. deflection in xz plane
 plot(x,EIz_x)
 xlabel('Axial Postion')
 ylabel('EI times Deflection')
 figure(3) %position vs. slope in xy plane
 plot(x,EIdy_x)
 xlabel('Axial Postion')
 ylabel('EI times Slope')
 figure(4) %position vs. slope in xz plane
 plot(x,EIdz_x)
 xlabel('Axial Postion')
 ylabel('EI times Slope')

 % Find total slope and deflection
delta_prime = sqrt((EIdy_x).^2+(EIdz_x).^2);
delta = sqrt((EIz_x).^2+(EIy_x).^2);

% plot vs. axial position
figure(5)
plot(x,delta_prime)
xlabel('Axial Postion')
ylabel('EI times Slope')
figure(6)
plot(x,delta)
xlabel('Axial Postion')
ylabel('EI times Deflection')
 
% Calculate Resultant Slopes and Deflections
delta_a = delta(1);
delta_c = delta(7);
theta_a = delta_prime(1);
theta_b = delta_prime(4);
theta_c = delta_prime(7);
theta_d = delta_prime(10);

% Diameter of the shaft due to slope at the bearings - Points B/D
D_b = ((theta_b*64)/(E*pi*bearing_slope))^(1/4);
D_d = ((theta_d*64)/(E*pi*bearing_slope))^(1/4);
d_shaft_1 = max(D_b,D_d);

% Diameter of the shaft due to slope at the gears - Points A/C
D_a = ((theta_a*64)/(E*pi*gear_slope))^(1/4);
D_c = ((theta_c*64)/(E*pi*gear_slope))^(1/4);
d_shaft_2 = max(D_a,D_c);

% Diameter of the shaft due to deflection at the gears - Points A/C
D_a1 = ((delta_a*64)/(E*pi*gear_def))^(1/4);
D_c1 = ((delta_c*64)/(E*pi*gear_def))^(1/4);
d_shaft_3 = max(D_a1,D_c1);

%Diameter of the shaft due to torsional wind-up
G = E/(2*(1+p)); 
d_shaft_4 = ((T*32)/(G*pi*dtheta_dx))^(1/4);

%Fatigue/Failure Criteria - Modified Goodman 
% Modified Endurance Limit 
ka = (2.7)*(sut*10^-3)^-0.265; %Cold-Drawn
kb = 0.879*1^-0.107; %diameter guess of 1 in
kc = 1; %due to combined loading
kd = 1; %ambient temperature
ke = 0.814; % 99% reliability
kf = 1; %given 
se_prime = 0.5*sut;
se=ka*kb*kc*kd*ke*kf*se_prime;

%Total Bending Moment  
M_tot = sqrt((M_xy).^2+(M_xz).^2);
M_max = max(M_tot); %Occurs at Point C 
% woodruff key: r/d = 0.02
% diameter guess = 1 in., so r = 0.02 in
q = 0.7; %From chart 6-20
qs = 0.77; % From chart 6-21
kt = 2.14; %given
kts = 3.0; %given
kf_notch = q*(kt-1)+1;
kfs =q*(kts-1)+1;
syms d
sig_a = (kf*M_max*32)/(pi*d^3);
sig_m = 0; 
tau_a = 0;
tau_m = (16*kfs*T)/(pi*d^3);
sig_a_prime = sqrt((sig_a)^2+3*(tau_a)^2);
sig_m_prime = sqrt((sig_m)^2+3*(tau_m)^2);
n = 2; % design factor

%Modified Goodman 
% (sig_a_prime/se)+(sig_m_prime/sut)=(1/n)
% ^ use this equation to derive diameter for d_shaft_5
d_shaft_5 = ((((32*M_max*kf_notch)/se)+(kfs*T*16*sqrt(3))/(sut))*(n/pi))^(1/3);


% Critical speed of shaft 
% w_operating(shaft speed) *nd = w1; 
g = 386.24; % in/sec^2
nd = 1.5; 
ss1 = ss*0.1047; % convert to rad/in
w1 = ss1 * nd; 
y1 = -EIy_x(1);
y2 = EIy_x(7);
d_shaft_6 = ((((w1)^2/g)*4096*(ga_w*y1^2+gc_w*y2^2))/((pi*E)*64*(ga_w*y1+gc_w*y2)))^(1/12);

% Determine required diameter for the shaft
d_shaft_total = [d_shaft_1,d_shaft_2,d_shaft_3,d_shaft_4,d_shaft_5,d_shaft_6];
d_shaft_required = max(d_shaft_total);