Ways for faster torque interpolation

1 次查看（过去 30 天）
Siddharth Jain
Siddharth Jain 2023-3-9
评论： Jan 2023-3-9
I am modelling a single stage helical gear transmission, in which I am applying an input torque to the driver gear. However, the torque interpolation in my code takes the longest time and hence a 60 sec simulation takes more than 12 hours to run. Is there a way to interpolate the torque faster or through a more efficient method?
My main code-
function [yp] = reduced_t(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;
% Common parameters
% e_a = zeros(size(t)); %circumferential displacement of the driver gear (m)
% e_p = zeros(size(t)); %circumferential displacement of the driver gear (m)
R_a = 51.19e-3; %Radius
R_p = 139.763e-3; %Radius
e_a = (pi*2*R_a*tan(beta))/(2*cos(P));
e_p = (pi*2*R_p*tan(beta))/(2*cos(P));
% for i = 2:length(t)
% % theta_a_vec(i) = theta_a_vec(i-1) - speed*2*pi*(t(i) - t(i-1)); % iterative assignment
% e_a(i) = e_a(i-1) + theta_a_vec(i)*R_a;
% e_p(i) = e_p(i-1) - theta_a_vec(i)*R_p;
% end
e = 0;
ks = 0.9e8 + 20000*sin(2*pi*Freq*t); %Shear stiffness
Cs = 0.1 + 0.001*sin(2*pi*Freq*t); %Shear damping
m_a = 0.5;
I_a = 0.5*m_a*(R_a^2); %Moment of inertia of driver gear (Kg.m3)
% Driven gear
m_p = 0.5; %mass of driven gear (kg)
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
% 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_p; %circumferential displacement at the meshing point on the driven gear (m)
yp = zeros(4,1); %vector of 4 equations
% 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*yp(3); %circumferential dynamic excitation force at the meshing point (N)
%Time interpolation
time = 0:0.00001:0.06;
Torque = interp1(time,T_a,t);
% Torque = spline(time,T_a,t);
Opp_Torque = 68.853; % Average torque value
% %Time interpolation using parfor loop
% time = 0:0.00001:0.6;
% Torque = zeros(size(t));
% parfor i=1:length(t)
% Torque(i) = interp1(time,T_a,t(i));
% end
% Opp_Torque = 68.853; % Average torque value
%Driver gear equations
yp(1) = y(2);
yp(2) = (Torque- Cs*cos(beta)*R_a*(R_a*y(2) + R_p*y(4)) - ks*cos(beta)*R_a*(R_a*y(1)+R_p*y(3)) -ks*cos(beta)*R_a*(e_a-e_p-e))/I_a;
%Driven gear equations
yp(3) = y(4);
yp(4) = (Opp_Torque*i - Cs*cos(beta)*R_p*(R_a*y(2) + R_p*y(4)) - ks*cos(beta)*R_p*(R_a*y(1)+ R_p*y(3)) -ks*cos(beta)*R_p*(e_a-e_p-e))/I_p;
end
My command window-
tic
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:6001); %Torque (Nm)
% speed = 1000/60;
Freq = 1000*20/60;
tspan=0:0.00001:0.06;
y0 = [0;104.719;0;104.719/3];
[t,y] = ode45(@(t,y) reduced_t(t,y,T_a),tspan,y0);
% TE = theta_p*R_p - theta_a_vec*R_a; %transmission error
toc
%Driver gear graphs
nexttile
plot(t,y(:,2))
ylabel('angular velocity (rad/s)')
xlabel('Time')
title('Driver gear angular velocity vs time')
axis padded
grid on
% nexttile
% plot(t,y(:,1))
% 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(:,4))
ylabel('angular velocity (rad/s )')
xlabel('Time')
title('Driven gear angular velocity vs time')
axis padded
grid on
% nexttile
% plot(t,y(:,3))
% ylabel('angular velocity (rad/s)')
% xlabel('Time')
% title('Driven gear angular velocity vs time')
% axis padded
% grid on
% nexttile
% plot(t,T_a(1:601))
% ylabel('Torque (Nm)')
% xlabel('Time (sec)')
% title('Torque vs time')
% axis padded
% grid on
My torque file - torque_for_Sid_60s.mat

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

回答（1 个）

Jan
Jan 2023-3-9
interp1 is slower than griddedInterpolant. You can created the interpolant outside the integration to save more time.
A sever problem remains: The linear interpolation is not smooth (deifferentable), but Matlab's ODE integrators are defined for smooth functions only. So you code can confuse the stepsize control. On one hand this can reduce the stepsizes dratsically and increase the computation times. On the other hand, the huge number of function evaluations accumulates the rounding errors. In consequence, the result can be dominated by rounding errors, such that the shown code is a poor random number generator.
Use a smooth function (e.g. a cubic interpolation).
  4 个评论
显示 2更早的评论
Siddharth Jain
Siddharth Jain 2023-3-9
torque_interp = griddedInterpolant(time, torque_range);
Torque = torque_interp(t);
should I use as above?
Torque = griddedInterpolant(time,T_a,t)
gives me error that the value should be scalar
Jan
Jan 2023-3-9
Does the first do, what you want? You cancompare it with the output of interp1. If so:
time = 0:0.00001:0.06;
torque_interp = griddedInterpolant(time, T_a, 'pchip'); % or 'spline'
[t,y] = ode45(@(t,y) reduced_t(t,y, torque_interp),tspan,y0);
function [yp] = reduced_t(t,y,torque_interp)
...
Torque = torque_interp(t);
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 2-D and 3-D Plots 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by