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

2 views (last 30 days)
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)

Answers (1)

VBBV
VBBV on 30 Nov 2022
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
  6 Comments
VBBV
VBBV on 30 Nov 2022
you have to update the function input arguments as below
function yp = reduced(t,y,T_a)
Check my updated code
Siddharth Jain
Siddharth Jain on 30 Nov 2022
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;

Sign in to comment.

Categories

Find more on Gears in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!