hello
here some suggestions and mods I did in your code
- remove the initial transient (t< 0.1 s)
- detrend data (so you don't have that large DC output i the fft spectrum)
- display fft magnitude in log scale
Now we see better the fft peaks
(you can use findpeaks to get them)
clearvars
clc
close all
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:600001); %Torque (Nm) chnages with time step
speed = 1000/60;
tspan=0:0.00001*5:6; %can try 4 or 5
y0 = [0;0;0;0;0;104.719;0;0;0;0;0;104.719/3];
tic
%These are constant for every call to reduced2 so compute them here and
%pass them into the modified function reduced3
time = 0:0.00001*5:6;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) spur1(t,y,Torque),tspan,y0);
TE= zeros(size(t));
KS = zeros(size(t));
for i=1:numel(t)
[~,TE(i),KS(i)] = spur1(t(i),y(i,:),Torque);
end
% MN : remove first 0.1 second (transient)
idx = (t<0.1);
t(idx) = [];
TE(idx) = [];
nexttile
plot(t,TE,"magenta")
xlabel('Time (s)')
ylabel('Transmission Error TE (meter)')
title('Transmission Error vs time')
axis padded
grid on
Fs = 1/(t(2)-t(1)); % Sampling frequency
L = length(t); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(detrend(TE), NFFT)/L; % detrend the signal before fft to remove the large DC (and near DC) fft output
f = Fs/2*linspace(0,1,NFFT/2+1);
nexttile;
semilogy(f, 2*abs(Y(1:NFFT/2+1)))
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of TE vs frequency')
axis padded
grid on
% nexttile
% plot(t,KS,"red")
% xlabel('Time (s)')
% ylabel('KS (N)')
% title('KS force vs time')
% axis padded
% grid on
newtime = toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yp,TE,KS] = spur1(t,y,Torque)
yp = zeros(12,1); %vector of 12 equations
beta = 32*(pi/180); %Helix angle (rad) /13
P = 18.5*(pi/180); %Pressure angle (rad)
speed = 1000/60;
Freq = 1000*20/60;
n_a = 22; % no of driver gear teeth
n_p = 59; % no of driven gear teeth
m = 2e-3*cos(beta);
R_a = (m*n_a)/2; %Radius
R_p = (m*n_p)/2; %Radius
i_r = n_p/n_a; % gear ratio
% Driver gear
m_a = 5.3e-2; %mass of driver gear (kg)
c_ay=6.3e2; %Damping of driver gear in y direction (Ns/m)
c_az=1.3e2; %Damping of driver gear in z direction (Ns/m)
k_ay= 1.9e8; %Stiffness of driver gear in y direction (N/m)
k_az= 8.3e6; %Stiffness of driver gear in z direction (N/m)
I_a = 1.7e-5; %Moment of inertia of driver gear (Kg.m3)
% Driven gear
m_p = 2.5e-1; %mass of driven gear (kg)
c_py=3.6e3; %Damping of driven gear in y direction (Ns/m)
c_pz=4.3e2; %Damping of driven gear in z direction (Ns/m)
k_py = 8.3e8; %Stiffness of driven gear in y direction (N/m)
k_pz = 1.2e7; %Stiffness of driven gear in z direction (N/m)
I_p = 3.9e-4; %Moment of inertia of driven gear (Kg.m3)
Torque = Torque(t)/1000;
Opp_Torque = 68.853/1000; % Average torque value
% e = (pi*2*(R_a + R_p)*tan(beta))/(4*cos(P));
e= 20e-6 + 20e-8*sin(2*pi*Freq*t);
ks = 0.9e8 + 0.9e6*sin(2*pi*Freq*t); %Shear stiffness
% Cs = 0.1 + 0.0001*sin(2*pi*Freq*t); %Shear damping
Cs = 2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p)) + 2e-2*0.1*sqrt((ks*I_a*I_p)/(I_a*R_a + I_p*R_p))*sin(2*pi*Freq*t);
TE = -y(11)*R_p + y(5)*R_a ; %transmission error
b = 20e-6; %Backlash
if(TE>b)
h = TE - b;
KS = h*ks;
elseif(TE < -1*b)
h = TE + b;
KS = h*ks;
else
h = 0;
KS = h*ks;
end
Fy = cos(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
Fz = sin(beta)*(Cs*(cos(beta)*(y(2) - y(8) + R_a*y(6) - R_p*y(12)) + sin(beta)*(y(4) - y(10))) + ks*(cos(beta)*(y(1) - y(7) + R_a*y(5) - R_p*y(11)) + sin(beta)*(y(3) - y(9)))); %circumferential excitation force- formula from paper
%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)/I_a; %angular acceleration
%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) = (-Opp_Torque*i_r + Fy*R_p)/I_p; % angular accelration (rad/s2)
end
