Getting an empty figure for ( amplitude vs rotor speed)

30 次查看（过去 30 天）
Amira
Amira 2025-8-22，16:17
编辑： Torsten 2025-8-22，22:42
Ran in:
% File name: mira.m
% This example has isotropic bearings
% A model with 4 Timoshenko beam elements
%
clear
format short e
close all
set(0,'defaultaxesfontsize',12)
set(0,'defaultaxesfontname','Times New Roman')
set(0,'defaulttextfontsize',12)
set(0,'defaulttextfontname','Times New Roman')
% set the material parameters
E = 2e11;
poisson = 0.3;
G = E/(2*(1+poisson));
rho = 7800;
damping_factor = 0; % no damping in shaft
% Consider the model with 3 equal length elements
% Shaft is 0.4m long
model.node = [1 0.0; 2 0.4/3; 3 2*0.4/3; 4 3*0.4/3];
% Assume shaft type 2 - Timoshenko with gyroscopic effects included
% Solid shaft with 30mm outside diameter
shaft_od = 0.020;
shaft_id = 0.0;
model.shaft = [2 1 2 shaft_od shaft_id rho E G damping_factor; ...
2 2 3 shaft_od shaft_id rho E G damping_factor; ...
2 3 4 shaft_od shaft_id rho E G damping_factor];
% Disk 1 at node 3 has diameter of 280mm and thickness of 70mm
% Note inside diameter of disk is assumed to be the outside diameter of the
% shaft
disk1_od = 0.3;
disk1_id = 0.1;
disk_thick = 0.03;
model.disc = [1 3 rho disk_thick disk1_od shaft_od];
% constant stiffness short isotropic bearing (1NM/m) with damping
% bearings at the ends - nodes 1 and 4
bear_stiff_1 = 5e7;
bear_stiff_2 = 7e7;
bear_stiff_3 = 0;
bear_stiff_4 = 0;
bear_damp_1 = 5e2;
bear_damp_2 = 7e2;
bear_damp_3 = 0;
bear_damp_4 = 0;
model.bearing = [3 1 bear_stiff_1 bear_stiff_2 bear_stiff_3 bear_stiff_4 bear_damp_1 bear_damp_2 bear_damp_3 bear_damp_4; ...
3 4 bear_stiff_1 bear_stiff_2 bear_stiff_3 bear_stiff_4 bear_damp_1 bear_damp_2 bear_damp_3 bear_damp_4];
% draw the rotor
figure(1), clf
picrotor(model)
% plot the Campbell diagram and root locus
Rotor_Spd_rpm = 0:1000:9000.0;
%Rotor_Spd = 2*pi*Rotor_Spd_rpm/60; % convert to rad/s
Rotor_Spd = Rotor_Spd_rpm/60; % convert to Hz
[eigenvalues,~,kappa] = chr_root(model,Rotor_Spd);
figure(2)
NX = 1.5;
damped_NF = 1; % plot damped natural frequencies
plotcamp(Rotor_Spd,eigenvalues,NX,damped_NF,kappa)
figure(3)
plotloci(Rotor_Spd,eigenvalues,NX)
% plot the modes at a given speed
Rotor_Spd_rpm = 3000;
Rotor_Spd = 2*pi*Rotor_Spd_rpm/60; % convert to rad/s
[eigenvalues,eigenvectors,kappa] = chr_root(model,Rotor_Spd);
figure(4)
subplot(221)
plotmode(model,eigenvectors(:,1),eigenvalues(1))
subplot(222)
plotmode(model,eigenvectors(:,3),eigenvalues(3))
subplot(223)
plotmode(model,eigenvectors(:,5),eigenvalues(5))
subplot(224)
plotmode(model,eigenvectors(:,7),eigenvalues(7))
% plot orbits
figure(5)
outputnode = [2 3];
axes('position',[0.2 0.53 0.2 0.2 ])
plotorbit(eigenvectors(:,1),outputnode,'Mode 1',eigenvalues(1))
axes('position',[0.39 0.53 0.2 0.2 ])
plotorbit(eigenvectors(:,3),outputnode,'Mode 2',eigenvalues(3))
axes('position',[0.58 0.53 0.2 0.2 ])
plotorbit(eigenvectors(:,5),outputnode,'Mode 3',eigenvalues(5))
axes('position',[0.2 0.25 0.2 0.2 ])
plotorbit(eigenvectors(:,7),outputnode,'Mode 4',eigenvalues(7))
% Choose DOF where you want the amplitude (e.g., disk node)
dof_measure = 5; % change this to the DOF you want to monitor
mu = 0.01; % unbalance mass [kg]
e = 1e-3; % eccentricity [m]
F0 = mu*e; % unbalance factor
amp = zeros(1,length(Rotor_Spd_rpm));
for i = 1:length(Rotor_Spd_rpm)
Omega = Rotor_Spd_rpm(i)*2*pi/60; % convert rpm -> rad/s
% Solve the dynamic system at each rotor speed
[eigenvalues,eigenvectors,kappa] = chr_root(model,Omega);
% Build mass, stiffness, damping matrices if needed
% Here we use eigenvectors as modal representation
% Simplified: apply unbalance at DOF 'dof_measure'
F = zeros(size(eigenvectors,1),1);
F(dof_measure) = F0*Omega^2;
% Solve for response in modal coordinates
q = eigenvectors * ((eigenvectors' * F) ./ diag(eigenvalues));
% Store amplitude (magnitude at chosen DOF)
amp(i) = abs(q(dof_measure));
end
Rotor_Spd_rpm
Rotor_Spd_rpm =
3000
amp
amp =
NaN
% Plot amplitude vs rotor speed
figure (6)
plot(Rotor_Spd_rpm, amp, 'r','LineWidth',1.5)
xlabel('Rotor speed [rpm]')
ylabel('Amplitude [m]')
title('Unbalance Response: Amplitude vs Rotor Speed')
grid on

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

回答（1 个）

Torsten
Torsten 2025-8-22，20:10
移动：Torsten 2025-8-22，20:10
You get amp = NaN as result, and NaN values are not plotted.
  2 个评论
Torsten
Torsten 2025-8-22，22:40
编辑：Torsten 2025-8-22，22:42
The code looks like an example code to the software for lateral vibrations. If it doesn't work, you should contact the authors for help.
There is a link to Michael Friswell under
http://www.rotordynamics.info/
If you made changes to the code and it worked in its original form, you should ask yourself why the changes resulted in NaN results.

请先登录，再进行评论。

类别

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

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by