Help needed to find error

*Please help me with my issue, it may seem that my question is repeated, but it is not. My code has been updating each week and now I face a different issue. I want help with plotting KS vs time and wether KS is being calculated properly in my code for each time step-

My main code-

function [yp] = reduced4(t,y,Torque)
beta = 13*(pi/180);                      %Helix angle (rad)
P = 20*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;
R_a = 51.19e-3;                             %Radius
R_p = 139.763e-3;                           %Radius
% Common parameters
e_a = (pi*2*R_a*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
e_p = (pi*2*R_p*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
% e = 0;
e = (pi*2*(R_a+R_p)*tan(beta))/(4*cos(P));
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
Cs = 0.1 + 0.0001*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)
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)
% 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)
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
% % 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  
Torque = Torque(t)/1000;
% torque_table = [time', T_a];
% torque = interp1(torque_table(:,1), torque_table(:,2), t, 'linear', 'extrap') / 1000;
Opp_Torque = 68.853/1000;                        % Average torque value
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
ks = 0.9e8 + 20000*sin(2*pi*Freq*t);                   %Shear stiffness 
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
yp(1) = y(2); 
yp(2) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - 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) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - 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 - Cs*cos(beta)*R_a*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_a*(R_a*y(5)+R_p*y(11)) -KS*cos(beta)*R_a*(e_a-e_p-e))/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - 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) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - 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 - Cs*cos(beta)*R_p*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_p*(R_a*y(5)+ R_p*y(11)) -KS*cos(beta)*R_p*(e_a-e_p-e))/I_p;                % angular accelration (rad/s2) 
end

My function call-

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:100001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*5:1;    %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:1;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) reduced4(t,y,Torque),tspan,y0); 
newtime = toc

% Driver gear graphs
nexttile  
plot(t,y(:,2))
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in y direction vs time')
axis padded  
grid on
nexttile  
plot(t,y(:,4))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,6))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driver gear angular velocity vs time')  
axis padded  
grid on
% Driven gear graphs
nexttile  
plot(t,y(:,8))  
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in y direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,10))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,12))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driven gear angular velocity vs time')  
axis padded  
grid on

The torque file is 21mb large and cannot be attached via the paperclip button - torque_for_Sid_60s.mat

It is my sincere request to please help me!

Answer 1

[Torsten]

You didn't include the error message. If it's a failure in integration at a certain time instant, I think we cannot help you. You know best the equations you implemented and what could be the cause of a difficulty.

We can only give you the usual advices:

• Try a different integrator
• Play with the error tolerances for the integrator ('RelTol' and 'AbsTol')
• Check your code for possible coding errors
• Try a simplified model first
• See if the error occurs at the time when the following discontinuity in your equations becomes active:

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

[Siddharth Jain]

No my problem is not with integration error (there was a typo earlier), I want to know how KS can be plotted vs time

[Torsten]

Make a copy of your function, change the output argument of the new function to KS and call this function in a loop with the solution you got from ode45.

[Siddharth Jain]

sorry, I don't get what you mean

[Torsten]

A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:5:100001);                  %Torque (Nm) chnages with time step
speed = 1000/60; 
tspan=0:0.00001*5:1;    %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:1;
Torque = griddedInterpolant(time, T_a);
[t,y] = ode45(@(t,y) reduced4(t,y,Torque),tspan,y0); 
KS = zeros(size(t));
for i=1:numel(t)
  [~,KS(i)] = reduced4(t(i),y(i,:),Torque)
end
plot(t,KS)
newtime = toc
% Driver gear graphs
nexttile  
plot(t,y(:,2))
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in y direction vs time')
axis padded  
grid on
nexttile  
plot(t,y(:,4))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driver gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,6))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driver gear angular velocity vs time')  
axis padded  
grid on
% Driven gear graphs
nexttile  
plot(t,y(:,8))  
ylabel('(y) up and down velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in y direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,10))  
ylabel('(z) side to side velocity (m/s)')
xlabel('Time')
title('Driven gear velocity in z direction vs time')  
axis padded  
grid on
nexttile  
plot(t,y(:,12))  
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driven gear angular velocity vs time')  
axis padded  
grid on

function [yp,KS] = reduced4(t,y,Torque)
beta = 13*(pi/180);                      %Helix angle (rad)
P = 20*(pi/180);                         %Pressure angle (rad)
speed = 1000/60; 
Freq = 1000*20/60;
R_a = 51.19e-3;                             %Radius
R_p = 139.763e-3;                           %Radius
% Common parameters
e_a = (pi*2*R_a*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
e_p = (pi*2*R_p*tan(beta))/(2*cos(P));      %circumferential displacement of the driver gear (m)
% e = 0;
e = (pi*2*(R_a+R_p)*tan(beta))/(4*cos(P));
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
Cs = 0.1 + 0.0001*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)
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)
% 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)
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
% % 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  
Torque = Torque(t)/1000;
% torque_table = [time', T_a];
% torque = interp1(torque_table(:,1), torque_table(:,2), t, 'linear', 'extrap') / 1000;
Opp_Torque = 68.853/1000;                        % Average torque value
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
ks = 0.9e8 + 20000*sin(2*pi*Freq*t);                   %Shear stiffness 
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
yp(1) = y(2); 
yp(2) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - 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) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - 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 - Cs*cos(beta)*R_a*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_a*(R_a*y(5)+R_p*y(11)) -KS*cos(beta)*R_a*(e_a-e_p-e))/I_a; %angular acceleration
%Driven gear equations
yp(7) = y(8);  
yp(8) = (-(Cs*cos(beta)*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*(R_a*y(5)+R_p*y(11)) - KS*cos(beta)*(e_a-e_p-e)) - 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) = (-(KS*sin(beta)*(z_a*tan(beta) - e_a*tan(beta) - R_a*y(5)*tan(beta) - z_p*tan(beta) - e_p*tan(beta) - R_p*y(11)*tan(beta) - z) + Cs*sin(beta)*(-tan(beta)*R_a*y(6) - tan(beta)*R_p*y(12))) - 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 - Cs*cos(beta)*R_p*(R_a*y(6) + R_p*y(12)) - KS*cos(beta)*R_p*(R_a*y(5)+ R_p*y(11)) -KS*cos(beta)*R_p*(e_a-e_p-e))/I_p;                % angular accelration (rad/s2) 
end

[Siddharth Jain]

Thanks for that, in the function, do you think that KS is being updated for each time step.

yp = zeros(12,1);                           %vector of 12 equations
ks = 0.9e8 + 20000*sin(2*pi*Freq*t);                   %Shear stiffness 
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

Also can you please help me with one last thing- can TE be plotted in the same way

[Torsten]

Thanks for that, in the function, do you think that KS is being updated for each time step.

In principle, KS depends on t, y(5) and y(11). So KS should be updated in each time step. Of course, if abs(TE) <= b for all times, KS will always be 0.

Also can you please help me with one last thing- can TE be plotted in the same way

KS = zeros(size(t));
TE = zeros(size(t));
for i=1:numel(t)
  [~,KS(i),TE(i)] = reduced4(t(i),y(i,:),Torque)
end
plot(t,[KS,TE])

function [yp,KS,TE] = reduced4(t,y,Torque)
...
end

[Siddharth Jain]

Hi @Torsten, I again need your help. Instead of posting another question, I am asking you here.

How do I plot the FFT of any function say TE in my function call?

[Torsten]

What do you mean by "plot ... in my function call" ?

And here is the call to fft:

https://de.mathworks.com/help/matlab/ref/fft.html

Just pass the vector "TE" to "fft".

[Siddharth Jain]

I meant that I want to check the frequency content of TE, so perfrom its FFT. I am currently doing this to do the FFT, is it correct?

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
nexttile
plot(t,TE,"magenta")
xlabel('Time (s)')
ylabel('Transmission Error TE (meter)')
title('Transmission Error vs time')
axis padded
grid on

TE_fft = fft(TE);
N = length(TE);
amp_spectrum = abs(TE_fft)/N;
amp_spectrum = amp_spectrum(1:N/2+1);
amp_spectrum(2:end-1) = 2*amp_spectrum(2:end-1);
freq = (0:N/2)*(1/(t(2)-t(1)))/N;

nexttile
plot(freq, amp_spectrum, 'magenta')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('FFT of Transmission Error')
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


% %Driver gear graphs
% nexttile  
% plot(t,y(:,2))
% ylabel('(y) up and down velocity (m/s)')
% xlabel('Time')
% title('Driver gear velocity in y direction vs time')
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,4))  
% ylabel('(z) side to side velocity (m/s)')
% xlabel('Time')
% title('Driver gear velocity in z direction vs time')  
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,6))  
% ylabel('angular displacement (rad)')
% xlabel('Time')
% title('Driver gear angular displacment vs time')  
% axis padded  
% grid on


% % Driven gear graphs
% nexttile  
% plot(t,y(:,8),"green")  
% ylabel('(y) up and down velocity (m/s)')
% xlabel('Time')
% title('Driven gear velocity in y direction vs time')  
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,10),"green")  
% ylabel('(z) side to side velocity (m/s)')
% xlabel('Time')
% title('Driven gear velocity in z direction vs time')  
% axis padded  
% grid on

% nexttile  
% plot(t,y(:,12),"green")  
% ylabel('angular displacement (rad)')
% xlabel('Time')
% title('Driven gear angular displacement vs time')  
% axis padded  
% grid on

[Torsten]

I am currently doing this to do the FFT, is it correct?

I have no experience with the usage and interpretation of fft. I think you will have to open a new question.