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
10 个评论
Walter Roberson
2023-2-21
Are you getting ode45 messages about having problems integrating?
Walter Roberson
2023-2-21
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
2023-2-21
编辑：Siddharth Jain
2023-2-21
Walter Roberson
2023-2-21
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
2023-2-21
Star Strider
2023-2-21
编辑：Star Strider
2023-2-21
Can you attach ‘torque_for_Sid_60s.mat’?
It will likely be necessary to rule out that it be the problem.
Siddharth Jain
2023-2-21
Star Strider
2023-2-21
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
2023-2-21
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.
回答（3 个）
Walter Roberson
2023-2-21
1 个投票
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
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.
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 个评论
Siddharth Jain
2023-2-22
编辑：Siddharth Jain
2023-2-22
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
2023-2-22
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
2023-2-22
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)
d_Fy_expr_dt = vpa(diff(Fy_expr,t), 5) % 'cos' Argument in 'F_y'
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)
d_Fz_expr_dt = vpa(diff(Fz_expr,t), 5) % 'sin' Argument in 'F_z'
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
2023-2-22
@Siddharth Jain, those are very impressive and thoughtful answers from @Star Strider and @Torsten and @Walter Roberson.
Siddharth Jain
2023-3-2
Star Strider
2023-3-2
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
2023-3-2
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.
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
类别
在 帮助中心 和 File Exchange 中查找有关 Assembly 的更多信息
产品
2023-2-21
2023-3-2
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
