Keep getting a blank graph

Question

Siddharth Jain 2023-2-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1916150-keep-getting-a-blank-graph

评论： Star Strider 2023-3-2

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Siddharth Jain 2023-2-21

It says max size 5mb, the file is 22mb. Can you please access it through the link?

Torsten 2023-2-21

编辑：Torsten 2023-2-21

@Siddharth Jain

As already noted several times before, Fy in function "reduced_t" is a scalar value (same as t). Thus your plot of Fy in function "reduced_t" will be a single point on a big empty white sheet.

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2023-2-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1916150-keep-getting-a-blank-graph#answer_1176275

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 个评论
显示无隐藏无

Siddharth Jain 2023-2-21

so do I use spline instead of interp1?

Walter Roberson 2023-2-21

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.

请先登录，再进行评论。

Answer 2

Star Strider 2023-2-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1916150-keep-getting-a-blank-graph#answer_1176540

在 MATLAB Online 中打开

torque0.06s.txt

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

A = 6001×2 table

Time torque(0:0.06s) _______ _______________ 0 67.657 1e-05 67.67 2e-05 67.683 3e-05 67.697 4e-05 67.712 5e-05 67.727 6e-05 67.743 7e-05 67.76 8e-05 67.777 9e-05 67.794 0.0001 67.812 0.00011 67.831 0.00012 67.849 0.00013 67.868 0.00014 67.888 0.00015 67.907

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

y1_limits = 1×2

-97.3687 0.0000

% TE = theta_p*R_p - theta_a_vec*R_a; %transmission error

toc

Elapsed time is 0.978614 seconds.

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 个评论
显示 7更早的评论隐藏 7更早的评论

Star Strider 2023-2-22

在 MATLAB Online 中打开

This is how I would caalculate them —

Considering that:

speed = 1000/60;

theta_a_vec = -speed*2*pi*t;

e_a = theta_a_vec*R_a;

e_p = -theta_a_vec*R_p;

and:

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)

and:

beta = 13*(pi/180); %Helix angle (rad)

z = e*tan(beta); %axial tooth displacement caused by internal excitation (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)

syms e R_a R_p t

beta = 13*(pi/180); %Helix angle (rad)

speed = 1000/60;

theta_a_vec = -speed*2*pi*t;

e_a = theta_a_vec*R_a;

e_p = -theta_a_vec*R_p;

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)

Fy_expr = beta*(y_ac - y_pc - e)

Fy_expr =

d_Fy_expr_dt = vpa(diff(Fy_expr,t), 5) % 'cos' Argument in 'F_y'

d_Fy_expr_dt =

z = e*tan(beta); %axial tooth displacement caused by internal excitation (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)

Fz_expr = beta*(z_ac - z_pc - z);

Fz_expr = vpa(Fz_expr, 5)

Fz_expr =

d_Fz_expr_dt = vpa(diff(Fz_expr,t), 5) % 'sin' Argument in 'F_z'

d_Fz_expr_dt =

I tend to let the Symbolic Math Toolbox do all the ‘heavy lifting’ largely because I am prone to algebra errors and it is not. It is always possible to get these as functions of the various arguments and then use that, if the arguments change.

I am not criticizing your derivations. I am only providing expressions for what I believe the arguments to the second terms of those experssions should be.

.

Siddharth Jain 2023-3-2

I meant for eg running the simulationn for 6 sec. Then, I need to make a new txt file

Star Strider 2023-3-2

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.

请先登录，再进行评论。

Answer 3

Torsten 2023-2-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1916150-keep-getting-a-blank-graph#answer_1176525

在 MATLAB Online 中打开

I don't understand why you ask the same question again.

It has already been anwered here:

https://de.mathworks.com/matlabcentral/answers/1915190-how-to-include-iterative-values-for-calculations?s_tid=prof_contriblnk

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Keep getting a blank graph

10 个评论
显示 8更早的评论隐藏 8更早的评论

回答（3 个）

2 个评论
显示无隐藏无

9 个评论
显示 7更早的评论隐藏 7更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Keep getting a blank graph

10 个评论 显示 8更早的评论隐藏 8更早的评论

回答（3 个）

2 个评论 显示 无隐藏 无

9 个评论 显示 7更早的评论隐藏 7更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

10 个评论
显示 8更早的评论隐藏 8更早的评论

2 个评论
显示无隐藏无

9 个评论
显示 7更早的评论隐藏 7更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论