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

1 次查看（过去 30 天）
Siddharth Jain
Siddharth Jain 2022-11-30
评论： Siddharth Jain 2022-11-30
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)

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

回答（1 个）

VBBV
VBBV 2022-11-30
Ran in:
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 个评论
显示 4更早的评论
VBBV
VBBV 2022-11-30
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 2022-11-30
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;

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Gears 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by