how to replace/import .mat file with parameter values?

I want to replace torque T_a with the data in .mat file (attached) for my code-

function yp = reduced(t,y)  
beta = 13*(pi/180);                      %Helix angle (rad)
speed = 1000/60; 
%theta_a_vec=19; %linspace(0,2*pi,1/speed); 
% for i = 0:1/speed:2*pi 
% theta_a_vec = i; 
% end 
% Common parameters
theta_a_vec = speed*2*pi*t;     %torsional angle of driver gear (rad)
theta_p = 22.9;                 %torsional angle of driven gear (rad)
e = 20e-6;                      %circumferential relative displacement of the teeth (m)
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
e_a = 48e-6;                    %circumferential displacement of the driver gear (m)
e_p = 48e-6;                    %circumferential displacement of the driver gear (m)
ks = 704e3;                     %Shear stiffness 
Cs = 0.1;                       %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)
T_a = 178*sin(theta_a_vec);     %Torque (Nm)
I_a = 7;                        %Moment of inertia of driver gear (Kg.m3)
R_a = 254.523e-3;               %Radius
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)
% 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)
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
yp(1) = y(2); 
yp(2) = (-Fy - 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) = (-Fz - 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) = (-T_a - Fy*R_a)/I_a;                % angular acceleration (rad/s2)
% 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= 8e7;                                  %Stiffness of driven gear in y direction (N/m)
k_pz =5e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 2000;                                 %Moment of inertia of driven gear (Kg.m3)
R_p = 944.2e-3;                             %Radius
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio
%Driven gear equations
yp(7) = y(8);  
yp(8) = (-Fy - 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) = (-Fz - 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) = (-T_a*i - Fy*R_p)/I_p;                       % angular accelration (rad/s2)
end  

My command window code-

speed = 1000/60; 
tspan=0:.00001:1/speed;  
%tspan = [0 12];  
y0 = zeros(12,1);  
[T,Y] = ode45(@reduced,tspan,y0);  
for i = 1:numel(T)  
    dy(i,:) = reduced(T(i),Y(i,:));  
end  
% Driver gear graphs
nexttile  
plot(T,dy(:,2))
ylabel('(y) up and down accelration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
axis padded  
grid on
nexttile  
plot(T,dy(:,4))  
ylabel('(z) side to side accelration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,6))  
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')  
axis padded  
grid on
% theta_a_vec = speed*2*pi*T; 
% T_a = 178*sin(theta_a_vec); 
% plot(T, T_a) 
% Driven gear graphs
nexttile  
plot(T,dy(:,8))  
ylabel('(y) up and down accelration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,10))  
ylabel('(z) side to side accelration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,12))  
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')  
axis padded  
grid on
nexttile
theta_a_vec = speed*2*pi*T;
T_a = 178*sin(theta_a_vec);
plot(T,T_a)

Answer 1

[VBBV]

A = load(websave('X','https://in.mathworks.com/matlabcentral/answers/uploaded_files/1213353/torque_for_Sid.mat'))

A = struct with fields:
    torque_em: [67.6572 67.6697 67.6830 67.6971 67.7119 67.7273 67.7433 67.7598 67.7769 67.7944 67.8123 67.8306 67.8493 67.8683 67.8876 67.9072 67.9270 67.9472 67.9676 67.9883 68.0093 68.0306 68.0521 68.0740 68.0962 68.1188 68.1417 68.1649 68.1886 … ]

speed = 1000/60; 
tspan=0:.00001:1/speed;  
%tspan = [0 12];  
y0 = zeros(12,1);  
[T,Y] = ode45(@reduced,tspan,y0)  

T = 6001×1

1.0e+00 * 0 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001

Y = 6001×12

0 0 0 0 0 0 0 0 0 0 0 0 -0.0000 -0.0035 0.0000 0.0001 -0.0000 -0.0001 -0.0000 -0.0035 0.0000 0.0001 -0.0000 -0.0000 -0.0000 -0.0141 0.0000 0.0002 -0.0000 -0.0003 -0.0000 -0.0141 0.0000 0.0002 -0.0000 -0.0000 -0.0000 -0.0317 0.0000 0.0003 -0.0000 -0.0006 -0.0000 -0.0317 0.0000 0.0003 -0.0000 -0.0000 -0.0000 -0.0559 0.0000 0.0003 -0.0000 -0.0011 -0.0000 -0.0559 0.0000 0.0003 -0.0000 -0.0000 -0.0000 -0.0862 0.0000 0.0004 -0.0000 -0.0016 -0.0000 -0.0862 0.0000 0.0004 -0.0000 -0.0000 -0.0000 -0.1221 0.0000 0.0005 -0.0000 -0.0024 -0.0000 -0.1221 0.0000 0.0005 -0.0000 -0.0000 -0.0000 -0.1630 0.0000 0.0006 -0.0000 -0.0032 -0.0000 -0.1630 0.0000 0.0006 -0.0000 -0.0000 -0.0000 -0.2081 0.0000 0.0006 -0.0000 -0.0042 -0.0000 -0.2081 0.0000 0.0006 -0.0000 -0.0001 -0.0000 -0.2568 0.0000 0.0007 -0.0000 -0.0054 -0.0000 -0.2568 0.0000 0.0007 -0.0000 -0.0001

for i = 1:numel(T)  
    dy(i,:) = reduced(T(i),Y(i,:));  
end  
% Driver gear graphs
nexttile  
plot(T,dy(:,2))
ylabel('(y) up and down accelration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
axis padded  
grid on
nexttile  
plot(T,dy(:,4))  
ylabel('(z) side to side accelration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,6))  
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')  
axis padded  
grid on
% theta_a_vec = speed*2*pi*T; 
% T_a = 178*sin(theta_a_vec); 
% plot(T, T_a) 
% Driven gear graphs
nexttile  
plot(T,dy(:,8))  
ylabel('(y) up and down accelration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,10))  
ylabel('(z) side to side accelration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,12))  
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')  
axis padded  
grid on
nexttile
theta_a_vec = speed*2*pi*T;
T_a = 178*sin(theta_a_vec); 
plot(T,A.torque_em(1:6001)) % replace T_a with torque_em data in mat file

function yp = reduced(t,y)  
beta = 13*(pi/180);                      %Helix angle (rad)
speed = 1000/60; 
%theta_a_vec=19; %linspace(0,2*pi,1/speed); 
% for i = 0:1/speed:2*pi 
% theta_a_vec = i; 
% end 
% Common parameters
theta_a_vec = speed*2*pi*t;     %torsional angle of driver gear (rad)
theta_p = 22.9;                 %torsional angle of driven gear (rad)
e = 20e-6;                      %circumferential relative displacement of the teeth (m)
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
e_a = 48e-6;                    %circumferential displacement of the driver gear (m)
e_p = 48e-6;                    %circumferential displacement of the driver gear (m)
ks = 704e3;                     %Shear stiffness 
Cs = 0.1;                       %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)
T_a = 178*sin(theta_a_vec);     %Torque (Nm)
I_a = 7;                        %Moment of inertia of driver gear (Kg.m3)
R_a = 254.523e-3;               %Radius
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)
% 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)
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
yp(1) = y(2); 
yp(2) = (-Fy - 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) = (-Fz - 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) = (-T_a - Fy*R_a)/I_a;                % angular acceleration (rad/s2)
% 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= 8e7;                                  %Stiffness of driven gear in y direction (N/m)
k_pz =5e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 2000;                                 %Moment of inertia of driven gear (Kg.m3)
R_p = 944.2e-3;                             %Radius
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio
%Driven gear equations
yp(7) = y(8);  
yp(8) = (-Fy - 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) = (-Fz - 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) = (-T_a*i - Fy*R_p)/I_p;                       % angular accelration (rad/s2)
end  

[VBBV]

The torque_em data variable in the mat file is sampled to a higher value. Can you be specific with replacing the variable data directly or by susbtitution into an equation for variable T in your code?

[Siddharth Jain]

Thnak you!

The torque value can be replaced directly as I do not have an equation for it

So instead of T_a, I want the data in the .mat file

[VBBV]

Ok, then you can follow the below approach

A = load(websave('X','https://in.mathworks.com/matlabcentral/answers/uploaded_files/1213353/torque_for_Sid.mat'))

A = struct with fields:
    torque_em: [67.6572 67.6697 67.6830 67.6971 67.7119 67.7273 67.7433 67.7598 67.7769 67.7944 67.8123 67.8306 67.8493 67.8683 67.8876 67.9072 67.9270 67.9472 67.9676 67.9883 68.0093 68.0306 68.0521 68.0740 68.0962 68.1188 68.1417 68.1649 68.1886 … ]

T_a = A.torque_em;
speed = 1000/60; 
tspan=0:.00001:1/speed;  
%tspan = [0 12];  
y0 = zeros(12,1);  

for k = 1:1000:length(T_a)
    [T,Y] = ode45(@(t,y) reduced(t,y,T_a(k)),tspan,y0);  
    for i = 1:numel(T)
       dy(i,:) = reduced(T(i),Y(i,:),T_a(k)); 
    end     
      
end  
size(dy)

ans = 1×2

6001 12

% Driver gear graphs
nexttile  
plot(T,dy(:,2))
ylabel('(y) up and down accelration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
axis padded  
grid on
nexttile  
plot(T,dy(:,4))  
ylabel('(z) side to side accelration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,6))  
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')  
axis padded  
grid on
% theta_a_vec = speed*2*pi*T; 
% T_a = 178*sin(theta_a_vec); 
% plot(T, T_a) 
% Driven gear graphs
nexttile  
plot(T,dy(:,8))  
ylabel('(y) up and down accelration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,10))  
ylabel('(z) side to side accelration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')  
axis padded  
grid on
nexttile  
plot(T,dy(:,12))  
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')  
axis padded  
grid on

% nexttile
theta_a_vec = speed*2*pi*T;
% % T_a = 178*sin(theta_a_vec); 
% plot(T,) % replace T_a with torque_em data in mat file

function yp = reduced(t,y,T_a)  
beta = 13*(pi/180);                      %Helix angle (rad)
speed = 1000/60; 
%theta_a_vec=19; %linspace(0,2*pi,1/speed); 
% for i = 0:1/speed:2*pi 
% theta_a_vec = i; 
% end 
% Common parameters
theta_a_vec = speed*2*pi*t;     %torsional angle of driver gear (rad)
theta_p = 22.9;                 %torsional angle of driven gear (rad)
e = 20e-6;                      %circumferential relative displacement of the teeth (m)
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
e_a = 48e-6;                    %circumferential displacement of the driver gear (m)
e_p = 48e-6;                    %circumferential displacement of the driver gear (m)
ks = 704e3;                     %Shear stiffness 
Cs = 0.1;                       %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)
% T_a = 178*sin(theta_a_vec);     %Torque (Nm)
I_a = 7;                        %Moment of inertia of driver gear (Kg.m3)
R_a = 254.523e-3;               %Radius
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)
% 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)
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
yp(1) = y(2); 
yp(2) = (-Fy - 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) = (-Fz - 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) = (-T_a - Fy*R_a)/I_a;                % angular acceleration (rad/s2)
% 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= 8e7;                                  %Stiffness of driven gear in y direction (N/m)
k_pz =5e7;                                  %Stiffness of driven gear in z direction (N/m)
I_p = 2000;                                 %Moment of inertia of driven gear (Kg.m3)
R_p = 944.2e-3;                             %Radius
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio
%Driven gear equations
yp(7) = y(8);  
yp(8) = (-Fy - 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) = (-Fz - 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) = (-T_a*i - Fy*R_p)/I_p;                       % angular accelration (rad/s2)
end  

[Siddharth Jain]

I am getting this error

Error using reduced

Too many input arguments.

Error in @(t,y)reduced(t,y,T_a(k))

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

[VBBV]

you have to update the function input arguments as below

function yp = reduced(t,y,T_a)  

Check my updated code

[Siddharth Jain]

Thanks!

One last thing. The torque in yp(6) and yp(12) is not chnaging the graphs behaviour.

yp(6) = (-T_a - Fy*R_a)/I_a;  
yp(12) = (-T_a*i - Fy*R_p)/I_p;