# Ways for faster torque interpolation

2 views (last 30 days)
Siddharth Jain on 9 Mar 2023
Commented: Jan on 9 Mar 2023
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)
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
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))
xlabel('Time')
title('Driver gear angular velocity vs time')
grid on
% nexttile
% plot(t,y(:,1))
% xlabel('Time')
% title('Driver gear angular velocity vs time')
% grid on
% Driven gear graphs
nexttile
plot(t,y(:,4))
xlabel('Time')
title('Driven gear angular velocity vs time')
grid on
% nexttile
% plot(t,y(:,3))
% xlabel('Time')
% title('Driven gear angular velocity vs time')
% grid on
% nexttile
% plot(t,T_a(1:601))
% ylabel('Torque (Nm)')
% xlabel('Time (sec)')
% title('Torque vs time')
% grid on
My torque file - torque_for_Sid_60s.mat

Jan on 9 Mar 2023
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).
Jan on 9 Mar 2023
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

### Categories

Find more on Programming in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!