Getting an empty figure for ( amplitude vs rotor speed)

% 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

Answer 1

[Torsten]

You get amp = NaN as result, and NaN values are not plotted.

[Amira]

@Torsten how shall solve NaN values and thanks in advance!!

[Torsten]

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.