Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
1
  • Loading branch information
Leahy committed Dec 8, 2017
1 parent 88b3482 commit 5bdd3c7
Show file tree
Hide file tree
Showing 3 changed files with 298 additions and 35 deletions.
161 changes: 142 additions & 19 deletions ME3227_FinalProject.asv
@@ -1,8 +1,8 @@
% 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)
sut = 118*10^3; %psi (ultimate stress)
sy = 102*10^3; %psi (yield stress)
p = 0.3; %(poisson's ratio)
ga_w = 4; %lbf - make as input value
ga_t = 50; %teeth number - input
Expand All @@ -25,38 +25,161 @@ 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);
Fan = Fat * tand(ga_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
c1 = (117*Fan-36*Rby+4.5*Fcn)/6;
c2 = 4.5*Fan - 3*c1;
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
for i = 0:length(i)
if (0<=i) && (i<=3)
EIdy_x = -Fan/2*(i)^2+c1;
EIy_x = -Fan/6*(i)^3+c1*i+c2;

elseif (3<i) && (i<=6)
EIdy_x = -Fan/2*(i)^2+Rby/2*(i-3)^2+c1;
EIy_x = -Fan/6*(i)^3+Rby/6*(i-3)^2+c1*i+c2;
elseif (6<i) && (i<=9)
EIdy_x = -Fan/2*(i)^2+Rby/2*(i-3)^2-Fcn/2*(i-6)^2+c1;
EIy_x = -Fan/6*(i)^3+Rby/6*(i-3)^2-Fcn/6*(i-6)^3+c1*i+c2;
x = linspace(0,9,10);
for i=0:9
if i<(3)
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
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
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
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
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*10^3)*sut^-0.265; %Cold-Drawn
kb = 0.879^-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 = 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

% Diameter of the shaft due to slope at the bearings - Points B,D
%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 = ((n^2*(((kf^2*M_max^2*32^2)/(pi^2*se^2))+((3*16^2*kfs^2*T^2)/(pi^2*sut^2)))))^(1/9);


% Critical speed of shaft
% w_operating(shaft speed) *nd = w1;
g = 386.24; % in/sec^2
nd = 1.5;
ss1 = ss*0.1047;
w1 = ss1 * nd;

158 changes: 142 additions & 16 deletions ME3227_FinalProject.m
Expand Up @@ -25,37 +25,163 @@ 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);
Fan = Fat * tand(ga_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
c1 = (117*Fan-36*Rby+4.5*Fcn)/6;
c2 = 4.5*Fan - 3*c1;
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
for x = linspace(0,9,10)
if (0<=x) && (x<=3)
EIdy_x = -Fan/2*(x)^2+c1;
EIy_x = -Fan/6*(x)^3+c1*x+c2;
elseif (3<x) && (x<=6)
EIdy_x = -Fan/2*(x)^2+Rby/2*(x-3)^2+c1;
EIy_x = -Fan/6*(x)^3+Rby/6*(x-3)^2+c1*x+c2;
elseif (6<x) && (x<=9)
EIdy_x = -Fan/2*(x)^2+Rby/2*(x-3)^2-Fcn/2*(x-6)^2+c1;
EIy_x = -Fan/6*(x)^3+Rby/6*(x-3)^2-Fcn/6*(x-6)^3+c1*x+c2;
x = linspace(0,9,10);
for i=0:9
if i<(3)
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
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
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
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
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*10^3)*sut^-0.265; %Cold-Drawn
kb = 0.879^-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 = 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

% Diameter of the shaft due to slope at the bearings - Points B,D
%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 = ((n^2*(((kf^2*M_max^2*32^2)/(pi^2*se^2))+((3*16^2*kfs^2*T^2)/(pi^2*sut^2)))))^(1/9);


% Critical speed of shaft
% w_operating(shaft speed) *nd = w1;
g = 386.24; % in/sec^2
nd = 1.5;
ss1 = ss*0.1047;
w1 = ss1 * nd;
y1 = -EIy_x(1);
y2 = EIy_x(7);

14 changes: 14 additions & 0 deletions parameters.m
@@ -0,0 +1,14 @@
function = parameters(ga_w,ga_t,ga_a,ga_d,gc_w,gc_t,gc_a,gc_d)
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
%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

end

0 comments on commit 5bdd3c7

Please sign in to comment.