1. System Modeling: The Jeffcott Rotor

The Jeffcott rotor is a classic representation for analyzing rotor-bearing systems, treating the rotor as a concentrated mass on a flexible, massless shaft. The following MATLAB implementation defines the state-space equations considering viscous damping, stiffness, and unbalance forces.

% System parameters
rotorMass = 12.5;       % kg
shaftStiffness = 8e5;   % N/m
dampingCoeff = 1200;    % Ns/m
frictionCoef = 0.08;
clearance = 1e-4;       % m
angularVel = 1200;      % rad/s

% State vector: [x, dx, y, dy]
rotorDynamics = @(t, s) [s(2); 
   (-dampingCoeff*s(2) - shaftStiffness*s(1) + unbalanceForce(t) + rubForces(s(1), s(3), 1))/rotorMass;
   s(4); 
   (-dampingCoeff*s(4) - shaftStiffness*s(3) + unbalanceForce(t) + rubForces(s(1), s(3), 2))/rotorMass];
   

2. Nonlinear Rub-Impact Formulation

Rub-impact occurs when the rotor displacement exceeds the radial clearance, triggering contact forces. A common approach involves applying Coulomb friction models proportional to the contact depth.

function F = rubForces(x, y, axis)
   radialDisp = sqrt(x^2 + y^2);
   if radialDisp > clearance
       % Coulomb friction calculation
       contactForce = -frictionCoef * shaftStiffness * (radialDisp - clearance);
       if axis == 1
           F = contactForce * (x / radialDisp);
       else
           F = contactForce * (y / radialDisp);
       end
   else
       F = 0;
   end
end
   

3. Numerical Integration and Signal Processing

Using the ode45 solver allows for efficient calculation of the system's time-history. Following simulation, frequency domain analysis is performed to isolate harmonic components.

% Simulation configuration
tRange = [0, 20];
initialConditions = [1e-5, 0, 1e-5, 0];
solverOptions = odeset('RelTol', 1e-7, 'MaxStep', 0.01);

[time, state] = ode45(rotorDynamics, tRange, initialConditions, solverOptions);

% FFT Analysis
sampleFreq = 1 / (time(2) - time(1));
spectrum = abs(fft(state(:,1))) / length(time);
frequencies = (0:length(time)-1) * (sampleFreq / length(time));
   

4. Nonlinear Dynamics Identification

To characterize transition regimes, such as transition to chaos, bifurcation diagram are generated by sweeping the rotor speed parameter and plotting steady-state Poincarè points.

speedSweep = 800:50:2000;
poincareMap = zeros(length(speedSweep), 50);

for i = 1:length(speedSweep)
   % Perform sub-simulation at specific RPM
   [~, S] = ode45(@(t,s) systemModel(t, s, speedSweep(i)), [0 10], init);
   poincareMap(i, :) = S(end-49:end, 1);
end
   

5. Best Practices for Numerical Stability

• Stiffness Matrix: Ensure that your stiffness and damping matricse are positive definite to avoid solver divergence.
• Integration Error: If encountering non-physical spikes in results, reduce RelTol to 1e-9 and verify that the rub-impact velocity threshold is correctly implemented.
• Spectral Leakage: When performing FFT on short-duration signals, apply a Hanning or Hamming window to minimize energy leakage at frequency boundaries.