Code covered by the BSD License

# Rotor Dynamics toolbox (RotFE)

### Izhak Bucher (view profile)

Toolbox models rotating, elastic shafts with disks

Rot=reduce(Rot)
```function Rot=reduce(Rot)
%Rot=reduce(rotor)
%
% perfrom model order reduction for a rotor
%
% - Using modal reduction
% - adds the rigid body modes to the total count
%
% I. Bucher, 8-7-98, Rev.  1.0
% I. Bucher, 28-7-99, Rev.  1.5
%
% INPUT:
%   Rot  (output of rotfe)
% with the following fields,
%   Reduct, Reduct.Nmodes,
% OUTPUT:
%   reduced order matrices M,G,K, D and projected right hand side arguments
%        Fu_cos  Fu_sin
%   Rot.Reduct.flag=1  marks the fact that a reduction took place

if isfield(Rot,'Reduct') % reduce model ?

%% Check for rigid body modes
dd=sort(eig((Rot.K+Rot.K.')/2));
Rot.Reduct.Nmodes=Rot.Reduct.Nmodes+ sum( abs(dd)<1e-2 ); %  add the rigid body modes
% to the total count
[vv dd]=eigl(Rot.K,Rot.M,Rot.Reduct.Nmodes,[1 .001 .001],1);

% mass orthogonalisation
vv=orthg(full(vv));
%mr=diag( Q.'*Rot.M*Q );
%vv=Q*diag( 1./sqrt(mr) );

Rot.M=vv'*Rot.M*vv;
Rot.K=vv'*Rot.K*vv;
Rot.G=vv'*Rot.G*vv;
Rot.Kst=vv'*Rot.Kst*vv;
Rot.KH=vv'*Rot.KH*vv;
Rot.D=vv'*Rot.D*vv;
if ~isempty(Rot.Fu_cos),   Rot.Fu_cos=vv.'*Rot.Fu_cos; end
if ~isempty(Rot.Fu_sin), Rot.Fu_sin=vv.'*Rot.Fu_sin; end

Rot.T=vv;
Rot.dd=sqrt(dd);
Rot.Reduct.flag=1;  % mark the fact that the M,G,K are reduced
Rot.Reduct.Type='modal';  % mark the fact that the M,G,K are reduced

end

function Q=orthg(A)
% find an orthonormal basis for A
%   % QR decomposition
[Q,R,E]=qr(A);
%   % Determine r = effective rank
tol = eps*norm(A,'fro');
r = sum(abs(diag(R)) > tol);
r = r(1); % fix for case where R is vector.
%   % Use first r columns of Q.
if r > 0
Q = Q(:,1:r);