Help needed to find error

Question

Siddharth Jain 2023-3-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1931000-help-needed-to-find-error

评论： Torsten 2023-3-28

*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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-3-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1931000-help-needed-to-find-error#answer_1195775

移动：Torsten 2023-3-18

在 MATLAB Online 中打开

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Torsten 2023-3-18

编辑：Torsten 2023-3-18

在 MATLAB Online 中打开

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 2023-3-28

编辑：Siddharth Jain 2023-3-28

在 MATLAB Online 中打开

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