Keep getting a blank graph

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

[Walter Roberson]

Are you getting ode45 messages about having problems integrating?

[Walter Roberson]

plot(t,Fy,'*')

You do not have hold on for that plot. It is constantly going to overwrite the previous plot. 36000 or more times unless the integration ends early.

[Siddharth Jain]

No Ode45 problems.

Instead of a * graph, how can I get a normal line graph. I added hold on

[Walter Roberson]

animatedline()

However you are assuming that each time is evaluated once in increasing order and assuming that you can therefore plot time vs value there and that you will get a nice line graph as it progresses through time. That is not at all how ode45 works. ode45 always evaluates at least 6 different locations per iteration.

Imagine you are hiking on a mountain in the dark, and that you have a pinpoint strobe flashlight. You point it in various directions and strobe and each time you see how high the terrain is at that one point. How can you walk safely? Well you can create a careful search pattern that includes deliberately pointing in directions you are not currently intending to walk, doing so because the results in that direction give you information about how steep the path is locally and how it is turning. You might perhaps want to create a map of everything you saw, but more likely would be that you only wanted to chart the path you ended up taking.

[Siddharth Jain]

I used animatedline() but still get a blank graph. My orignal code has 6 Dof, but the modified code I pasted here is considering only 2DOF as I am trying to find the problems I am havihng with it

[Star Strider]

Can you attach ‘torque_for_Sid_60s.mat’?

It will likely be necessary to rule out that it be the problem.

[Siddharth Jain]

https://lunet-my.sharepoint.com/:u:/g/personal/wssj4_lunet_lboro_ac_uk/EVA2dOl4fShApZJz1KMN2DYBNlSIaYytmcf_IGZE9menBg?e=JAkr9e

[Star Strider]

You should be able to upload it here using the ‘paperclip’ icon just to the right of the Σ icon in the top Toolstrip. I am reluctant to use outside websites for these purposes.

[Siddharth Jain]

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

[Torsten]

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

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.

[Siddharth Jain]

so do I use spline instead of interp1?

[Walter Roberson]

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]

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.

.

[Siddharth Jain]

Hi, Thank you for your help.

Can you also help me with this-

How can I add the time derivative part to Fy and Fz, currently I did the derivative on paper and then added it to the code (speed*2*pi). But I feel it is the wrong way of doing it.

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

The research paper I am follwing says this-

[Star Strider]

Part of it I understand.

I believe the first terms in those expressions should be:

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

I assume that the derivatives are calculated elsewhere, so I have no idea how to correct the second term in those expressions to include them. (although ‘beta’ should multiply them).

I would use the Symbolic Math Toolbox to be certain that the derivative expressions are correct.

I leave that to you, because you understand your code. (Mechaincal engineering is not an area of my expertise.)

[Siddharth Jain]

no I just did the derivatives part on paper and and added them to the code, is there way to do the derivative in the equation itself. so for instance, do the time dertivative of (y_ac-y_pc-e) inside the Fy equation

[Star Strider]

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.

.

[William Rose]

@Siddharth Jain, those are very impressive and thoughtful answers from @Star Strider and @Torsten and @Walter Roberson.

[Siddharth Jain]

Hi @Star Strider

I am trying to convert the torque file for other time simulations as well from .mat to .txt. Could you please explain how you did this for 0.06 sec?

[Star Strider]

I do not understand.

The file already is a .txt file. If you want to do it for other simulations, first load the .mat file into your workspace and then use writetable to write it to a .txt file.

The ‘time’ vector is derived from the ‘torque0.06s.txt’ file as the first column, and goes from 0 to 0.06 seconds.. The code uses that as the ‘tspan’ argument to the ode45 call.

[Siddharth Jain]

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

[Star Strider]

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]

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