Keep getting a blank graph
1 view (last 30 days)
Show older comments
My code is trying to model the torsional degrees of freedom of a helical gear pair. I am unable to plot Force Fy vs time as it keeps coming up as a blank graph. My main code is-
function [yp] = reduced_t(t,y,T_a)
beta = 13*(pi/180); %Helix 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
theta_a_vec = -speed*2*pi*t;
e_a = theta_a_vec*R_a;
e_p = -theta_a_vec*R_p;
% theta_a_vec = zeros(size(t)); %torsional angle of driver gear (rad)
% theta_a_vec(1) = 0;
% e_a(1) = 0;
% e_p(1) = 0;
% 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
% theta_p = 22.9; %torsional angle of driven gear (rad)
theta_p = pi/6 + speed*2*pi*t;
%e = 20e-6; %circumferential relative displacement of the teeth (m)
e = 0;
z = e*tan(beta); %axial tooth displacement caused by internal excitation (m)
ks = 0.9e3 + 20*sin(2*pi*Freq*t); %Shear stiffness
Cs = 0.1 + 0.01*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)
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)
z_a = e_a*tan(beta); % axial displacement of driver gear (m)
z_p = e_p*tan(beta); % axial displacement of driven gear (m)
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)
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)
plot(t,Fy,'*')
axis padded
grid on
%Time interpolation
time = 0:0.00001:0.06;
% Torque = interp1(time,T_a,t);
Torque = interp1(time,T_a,t);
% sampling_frequency = 100e3;
% dT = 1/sampling_frequency;
% t_max = 0.6;
% time = 0:dT:t_max;
% Torque = T_a(1:numel(time));
Opp_Torque = 68.853; % Average torque value
%Driver gear equations
yp(1) = y(2);
yp(2) = (Fy*R_a)/I_a;
% 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= 9.5e7; %Stiffness of driven gear in y direction (N/m)
k_pz =9e7; %Stiffness of driven gear in z direction (N/m)
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
%Driven gear equations
yp(3) = y(4);
yp(4) = (Fy*R_p)/I_p; % angular accelration (rad/s2)
end
My command window code is-
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,0,0,0];
[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(:,1))
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('Driver gear angular acceleration vs time')
axis padded
grid on
%Driven gear graphs
nexttile
plot(t,y(:,3))
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('Driven gear angular acceleration vs time')
axis padded
grid on
nexttile
plot(t,T_a(1:6001))
ylabel('Torque (Nm)')
xlabel('Time (sec)')
title('Torque vs time')
axis padded
grid on
Answers (3)
Walter Roberson
on 21 Feb 2023
You use linear interpolation inside your ode function. Linear interpolation has discontinuous derivatives. The mathematics of Runge Kutta methods requires that the first two derivatives are continuous.
2 Comments
Walter Roberson
on 21 Feb 2023
That would be mathematically robust, but not necessarily physically reasonable. It has the physical implication that at time T-delta that you are already pretty close to where you will be at time T, a smooth quadratic transition at most. This is fine for some models but not fine for impulses. If the torque at time T is not knowable at time T-delta then spline is not the right model.
If you have what is effectively impulses, then instead you need to process only one interval per ode45 call, so that the torque is constant slope over the call.
Star Strider
on 21 Feb 2023
First, you only use the first 0.06 seconds of the torque vector, so create a file with just those values. (I did that offline, as ‘torque0.06s.txt’ and uploaded it.)
Second, the strange plot result was due to your calling plot within your ‘reduced_t’ function. It plotted a single ‘*’ in each iteration because that is what you told it to do. (I commented it out.)
The other plots appear to be appropriate.
I changed ‘reduced_t’ and the call to it to incorporate the imported table data —
tic
A = readtable('torque0.06s.txt', 'VariableNamingRule','preserve')
time = A{:,1};
T_a = A{:,2}; %Torque (Nm)
speed = 1000/60;
Freq = 1000*20/60;
% tspan=0:0.00001:0.06;
tspan = time;
y0 = [0,0,0,0];
[t,y] = ode45(@(t,y) reduced_t(t,y,time,T_a),tspan,y0);
y1_limits = [min(y(:,1)) max(y(:,1))]
% TE = theta_p*R_p - theta_a_vec*R_a; %transmission error
toc
figure
tiledlayout(3,1)
%Driver gear graphs
nexttile
plot(t,y(:,1))
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('Driver gear angular acceleration vs time')
axis padded
grid on
%Driven gear graphs
nexttile
plot(t,y(:,3))
ylabel('angular accelration (rad/s2)')
xlabel('Time')
title('Driven gear angular acceleration vs time')
axis padded
grid on
nexttile
plot(t,T_a)
ylabel('Torque (Nm)')
xlabel('Time (sec)')
title('Torque vs time')
axis padded
grid on
function [yp] = reduced_t(t,y,time,T_a)
beta = 13*(pi/180); %Helix 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
theta_a_vec = -speed*2*pi*t;
e_a = theta_a_vec*R_a;
e_p = -theta_a_vec*R_p;
% theta_a_vec = zeros(size(t)); %torsional angle of driver gear (rad)
% theta_a_vec(1) = 0;
% e_a(1) = 0;
% e_p(1) = 0;
% 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
% theta_p = 22.9; %torsional angle of driven gear (rad)
theta_p = pi/6 + speed*2*pi*t;
%e = 20e-6; %circumferential relative displacement of the teeth (m)
e = 0;
z = e*tan(beta); %axial tooth displacement caused by internal excitation (m)
ks = 0.9e3 + 20*sin(2*pi*Freq*t); %Shear stiffness
Cs = 0.1 + 0.01*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)
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)
z_a = e_a*tan(beta); % axial displacement of driver gear (m)
z_p = e_p*tan(beta); % axial displacement of driven gear (m)
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)
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)
% plot(t,Fy,'*') % Please Do Not Use 'plot' Calls Within Your ODE Function!
% axis padded
% grid on
%Time interpolation
% time = 0:0.00001:0.06;
% Torque = interp1(time,T_a,t);
Torque = interp1(time,T_a,t);
% sampling_frequency = 100e3;
% dT = 1/sampling_frequency;
% t_max = 0.6;
% time = 0:dT:t_max;
% Torque = T_a(1:numel(time));
Opp_Torque = 68.853; % Average torque value
%Driver gear equations
yp(1) = y(2);
yp(2) = (Fy*R_a)/I_a;
% 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= 9.5e7; %Stiffness of driven gear in y direction (N/m)
k_pz =9e7; %Stiffness of driven gear in z direction (N/m)
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
%Driven gear equations
yp(3) = y(4);
yp(4) = (Fy*R_p)/I_p; % angular accelration (rad/s2)
end
This appears to work correctly. If there are problems with it, I will do my best to help you solve them.
.
9 Comments
Star Strider
on 2 Mar 2023
The time seems to be dependent on the data in the files, since those are both used in the ordinary diffrential equaiton function. You will have to create them for the required time, since I have no idea how those data are calculated.
Torsten
on 21 Feb 2023
I don't understand why you ask the same question again.
It has already been anwered here:
t=0:0.00001:0.06;
beta = 13*(pi/180); %Helix 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
theta_a_vec = -speed*2*pi*t;
e_a = theta_a_vec*R_a;
e_p = -theta_a_vec*R_p;
% theta_a_vec = zeros(size(t)); %torsional angle of driver gear (rad)
% theta_a_vec(1) = 0;
% e_a(1) = 0;
% e_p(1) = 0;
% 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
% theta_p = 22.9; %torsional angle of driven gear (rad)
theta_p = pi/6 + speed*2*pi*t;
%e = 20e-6; %circumferential relative displacement of the teeth (m)
e = 0;
z = e*tan(beta); %axial tooth displacement caused by internal excitation (m)
ks = 0.9e3 + 20*sin(2*pi*Freq*t); %Shear stiffness
Cs = 0.1 + 0.01*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)
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)
z_a = e_a*tan(beta); % axial displacement of driver gear (m)
z_p = e_p*tan(beta); % axial displacement of driven gear (m)
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)
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)
plot(t,Fy)
grid on
0 Comments
See Also
Categories
Find more on Gears in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!