- 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:
Help needed to find error
1 次查看(过去 30 天)
显示 更早的评论
*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 个评论
回答(1 个)
Torsten
2023-3-18
移动:Torsten
2023-3-18
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:
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 个评论
Torsten
2023-3-28
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.
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!