MATLAB suite for simulating the lateral dynamics of rotating machines.
Define a rotating machine with a flexible rotor as follows:
P.Rotor = {};
P.Bearing = {};
P.Stator = {};
% Rotor
P.Rotor{1}.Name = 'Main';
P.Rotor{1}.Speed = 1;
P.Rotor{1}.Nodes = [0, 0.25, 0.5, 0.75, 1.0]; % Nodes of the FEM
%% Flexible shaft
P.Rotor{1}.Shaft = {};
P.Rotor{1}.Shaft{1}.Name = 'Shaft';
P.Rotor{1}.Shaft{1}.iNodes = 1:5;
P.Rotor{1}.Shaft{1}.Section.ro = 0.0;
P.Rotor{1}.Shaft{1}.Section.ri = 0.025;
P.Rotor{1}.Shaft{1}.Material.Name = 'steel';
P.Rotor{1}.Shaft{1}.Material.E = 210E9;
P.Rotor{1}.Shaft{1}.Material.rho = 7800;
%% Disc
P.Rotor{1}.Disc = {};
P.Rotor{1}.Disc{1}.Name = 'EndCap';
P.Rotor{1}.Disc{1}.iNode = 3;
P.Rotor{1}.Disc{1}.Type = 'Rigid';
P.Rotor{1}.Disc{1}.Options.bGyro = 1;
P.Rotor{1}.Disc{1}.Material.Name = 'steel';
P.Rotor{1}.Disc{1}.Material.E = Inf; %Rigid disc
P.Rotor{1}.Disc{1}.Material.rho = 7800;
P.Rotor{1}.Disc{1}.Ring.Geometry.R = [0 0.3];
P.Rotor{1}.Disc{1}.Ring.Geometry.t = 15E-3;
% Bearings
P.Bearing{1}.Name = 'LeftBearing';
P.Bearing{1}.kyy = 1E3;
P.Bearing{1}.kxx = 1E3;
P.Bearing{1}.cxx = 1E-3;
P.Bearing{1}.cyy = 1E-3;
P.Bearing{1}.Node{1}.Type = 'rotor';
P.Bearing{1}.Node{1}.iNode = 1;
P.Bearing{1}.Node{1}.iRotor = 1;
P.Bearing{1}.z = 0.0;
P.Bearing{2}.Name = 'RightBearing';
P.Bearing{2}.kyy = 1E3;
P.Bearing{2}.kxx = 1E3;
P.Bearing{2}.cxx = 1E-3;
P.Bearing{2}.cyy = 1E-3;
P.Bearing{2}.Node{1}.Type = 'rotor';
P.Bearing{2}.Node{1}.iNode = 5;
P.Bearing{2}.Node{1}.iRotor = 1;
P.Bearing{2}.z = 1.0;
% Excitation
P.Excite{1}.Name = 'Unbalance';
P.Excite{1}.Type = 'unbalance';
P.Excite{1}.Mode = 'sync';
P.Excite{1}.iRotor = 1;
P.Excite{1}.iDisc = 1;
P.Excite{1}.m = 0.01;
P.Excite{1}.r = 0.1;
P.Excite{1}.Angle = 0.0;
P = setupsystem(P);You can visualise the system with:
plot_system(P);which will generate a plot of the rototing machine like this:
To compute the natural frequencies at a range of speeds:
O = linspace(0,50,100)*2*pi; % Speed in rad/s
[V,d,W] = rotor_eig(P.Model,O);where:
VandWcontain the left and right eigevectorsdcontains the eigenvalues
at each speed.
To convert these eigenvalues and eigenvectors vectors into natural frequencies and mode shapes, run the following:
iNode = 2;
[omega,zeta, modes] = eig2modes(V,d,W);
% sort the modes
[~,ii] = sort(omega(:,1));
omega = omega(ii,:);
zeta = zeta(ii,:);
modes = modes(:,ii,:);
kappa = compute_whirl(P,modes,iNode);where:
omegacontains the natural frequencies in rad/szetacontains the non-dimensionalised damping ratioskappacontains the whirl direction, defined in terms of the motion at the central disc (node 2).
To plot the Campbell diagram with the first 2 modes, run the following
options.xunits = 'Hz';
iPlot = 1:2;
plot_campbell(P,O,omega(iPlot,:),zeta(iPlot,:),kappa(:,iPlot,:),iNode, options); which will generate a plot of the natural frequencies as a function of rotation speed like this:
To plot the shape of the response at each resonance (ie each mode shape), run the following:
plot_whirl(P,O,modes(:,iPlot,:),omega(iPlot,:),kappa(:,iPlot,:));
There is a slider which you can use to change which rotation speed the modes are plotted from.
Compute the response at all frequencies with:
[U,Udot,Uddot] = excitation_frf(P,O);
F = P.Model.Excite.K*U + P.Model.Excite.C*Udot + P.Model.Excite.M*Uddot;
q = rotor_frf_modal(V,d,W,F,O);And then plot the FRF of the disc response (DOF 9 and 10) using:
options.xunits = 'Hz';
options.iPlot = 9:10;
plot_frf(P,q,O,options);