Help with ODE tolerance options

1 次查看（过去 30 天）
Siddharth Jain
Siddharth Jain 2023-3-15
编辑： Siddharth Jain 2023-3-16
I am modelling a single stage helical gear transmission, and want to adjust the ODE tolerance options to allow for faster simulation times. Is there a way to figure out the tolerance values accurate for my model. I have tried various values for REltol and Abstol, but cannot get accurate results with that.
My main code-
function [yp] = reduced2(t,y,T_a)
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
time = 0:0.00001*2:10;
% Torque = interp1(time,T_a,t)/1000;
% torque_range = T_a;
% torque_interp = griddedInterpolant(time, torque_range);
% Torque = torque_interp(t) / 1000;
Torque = griddedInterpolant(time, T_a);
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 command window code-
tic
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:2:1000001); %Torque (Nm) chnages with time step
speed = 1000/60;
tspan=0:0.00001*2:10; %can try 4 or 5
y0 = [0;0;0;0;0;-104.719;0;0;0;0;0;104.719/3];
opts = odeset('RelTol',1e-1,'AbsTol',1e-1);
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);
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
% nexttile
% plot(t,T_a(1:60001))
% ylabel('Torque (Nm)')
% xlabel('Time (sec)')
% title('Torque vs time')
% axis padded
% grid on
My torque file- torque_for_Sid_60s - Copy.mat
  4 个评论
显示 2更早的评论
Jan
Jan 2023-3-15
If run time matters, move the creation of the griddedInterpolant out of the loop, as I've suggested in your other question.
A fast dirty solution is not useful, because the applied method of integrating a non-smooth function causes numerical troubles. Starting with the fine-tuning by choosing optimal tolerances is not the point: The approach is not accurate.
Mike Croucher
Mike Croucher 2023-3-15
I didn't notice you had already suggested that @Jan - sorry. That's pretty much all my answer below is about.

请先登录，再进行评论。

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

回答（1 个）

Mike Croucher
Mike Croucher 2023-3-15
编辑：Mike Croucher 2023-3-15
Thanks so much for including all of your code and data. It allowed me to run everything on my machine and use the profiler to figure out what was going on. In your function reduced2, you do this
time = 0:0.00001*2:10;
Torque = griddedInterpolant(time, T_a);
These are both constants for every call to reduced2 that ode45 makes. ode45 called reduced2 over 10,000 times and every time it recomputes these constant expressions. So, let's lift them out and just compute them once.
Instead of doing
tic
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);
oldtime = toc
I do
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*2:10;
Torque = griddedInterpolant(time, T_a);
[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0,opts);
newtime = toc
I've created a new function reduced3 that doesn't take T_a. The only reason you needed T_a was to create Torque. I've done that outside the function and sent in in as an argument. I now get the following times.
oldtime =
31.9675
newtime =
0.1238
Speedup is 258.130906x
I've also checked that all results are equal. Here is my reduced3
function [yp] = reduced3(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
and here is how we use it with your original version included for comparison
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:2:1000001); %Torque (Nm) chnages with time step
speed = 1000/60;
tspan=0:0.00001*2:10; %can try 4 or 5
y0 = [0;0;0;0;0;-104.719;0;0;0;0;0;104.719/3];
opts = odeset('RelTol',1e-1,'AbsTol',1e-1);
tic
[t,y] = ode45(@(t,y) reduced2(t,y,T_a),tspan,y0,opts);
oldtime = toc
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*2:10;
Torque = griddedInterpolant(time, T_a);
[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0,opts);
newtime = toc
fprintf("Speedup is %fx\n",oldtime/newtime)
fprintf("Are all results equal?\n")
all(tnew==t)
all(ynew==y,'all')
% 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
  10 个评论
显示 8更早的评论
Mike Croucher
Mike Croucher 2023-3-15
Sorry.
Change
[tnew,ynew] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0);
to
[t,y] = ode45(@(t,y) reduced3(t,y,Torque),tspan,y0);
Siddharth Jain
Siddharth Jain 2023-3-16
编辑：Siddharth Jain 2023-3-16
Hi @Mike Croucher, thank you for all your help. Can you also help me with plotting KS vs time in my code
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Gears 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by