Ways for faster torque interpolation

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

Answer 1

[Jan]

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).

[Siddharth Jain]

thank you for your answer, I am getting errors when i am using griddedInterpolant. How can i create the interpolant outside?

I had also tried cubic interpolation with spline, but it also takes lot of time.

[Jan]

@Siddharth Jain: Please post your code and the error messages. It is easier to fix an error thn to guess, what the error is.

You create the griddedInterpolant exactly like it is shown inthe examples in the documentation. Then you provide the interpolant object as parameter, not T_a.

[Siddharth Jain]

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]

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