Untitled

function [yp] = reduced2(t,y,T_a)  
beta = 13*(pi/180);                      %Helix angle (rad)
speed = 1000/60;                         % Speed in Revs/sec

Freq = 1000*20/60;                       %Speed divided by gear ratio

% Common parameters
theta_a_vec = -speed*2*pi*t;     %torsional angle of driver gear (rad)                    
e = 0;                          %circumferential relative displacement of the teeth (m)
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
e_a = 48e-6;                    %circumferential displacement of the driver gear (m)
e_p = 48e-6;                    %circumferential displacement of the driver gear (m)
ks = 0.9e3 + 20*sin(2*pi*Freq*t);                   %Shear stiffness 
Cs = 0.1 + 0.001*sin(2*pi*Freq*t);                   %Shear damping


% Driver gear  
m_a = 0.5;                      %mass of driver gear (kg)
c_ay=500;                       %Damping of driver gear in y direction (Ns/m)
c_az=500;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 8e7;                      %Stiffness of driver gear in y direction (N/m)
k_az= 5e7;                      %Stiffness of driver gear in z direction (N/m)
R_a = 51.19e-3;                 %Radius
I_a = 0.5*m_a*(R_a^2);          %Moment of inertia of driver gear (Kg.m3)
y_ac= e_a + theta_a_vec*R_a;    %circumferential displacement at the meshing point on the driver gear (m)
y_pc = e_p - theta_a_vec*R_a;   %circumferential displacement at the meshing point on the driven gear (m)
z_a = e_a*tan(beta);  
z_p = e_p*tan(beta);  
z_ac = (z_a-y_ac)*tan(beta);    %axial displacements of the meshing point on the driving gear  (m)
z_pc = (z_p+y_pc)*tan(beta);    %axial displacements of the meshing point on the driven gear   (m)

yp = zeros(12,1);                           %vector of 12 equations
% Excitation forces 
Fy=ks*cos(beta)*(y_ac-y_pc-e) + Cs*cos(beta)*2*R_a*speed*2*pi;                %axial dynamic excitation force at the meshing point (N)
Fz=ks*sin(beta)*(z_ac-z_pc-z) - Cs*sin(beta)*2*tan(beta)*R_a*speed*2*pi;      %circumferential dynamic excitation force at the meshing point (N)


%Time interpolation - Needed for solution of equations (It basically uses
%your torque and prescribed time matrices to generate a time matrix to be
%used in the ODE solver)
time = 0:0.00001:0.06;
Torque = interp1(time,T_a,t);
Opp_Torque = 68.853;                        % Average torque value

%Driver gear equations
yp(1) = y(2); 
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (Torque - Fy*R_a + Opp_Torque)/I_a;        % angular acceleration theta_a double dot (rad/s2)

% Driven gear  
m_p = 0.5;                                  %mass of driven gear (kg)
c_py=500;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=500;                                   %Damping of driven gear in z direction (Ns/m)
k_py= 9.5e7;                                %Stiffness of driven gear in y direction (N/m)
k_pz =9e7;                                  %Stiffness of driven gear in z direction (N/m)
R_p = 139.763e-3;                           %Radius
I_p = 0.5*m_p*(R_p^2);                      %Moment of inertia of driven gear (Kg.m3)
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio


%Driven gear equations
yp(7) = y(8);  
yp(8) = (-Fy - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (-Fz - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (-Torque*i - Fy*R_p + Opp_Torque*i)/I_p;      % angular accelration theta_p double dot (rad/s2)
end