Analyze the eigenvalues of a system with periodic coefficients.

I have a four-degree-of-freedom system for which I must consider all parameters dimensionless and ultimately analyze its stability boundaries. I used the ode45 function for numerical solutions and directly obtained the eigenvalues of the final solution matrix. I am very suspicious of my eigenvalues because they are very similar to each other and their values are extremely small. Can someone help me determine if I made a mistake in the dimensionless conversion? And is my script correct or not?

clc
clear all ;

% Define parameters 
N = 2 ; %Number of blades

I_thetaoverI_b =2 ; % Moment of inertia pitch axis over I_b 
I_psioverI_b =2 ; % Moment of inertia yaw axis over I_b 

C_thetaoverI_b = 0.00; % Damping coefficient over I_b 
C_psioverI_b = 0.00; % Damping coefficient over I_b 

h = 0.3; % aircraft center of gravity to rotor hub

hoverR = 0.34;
R = h / hoverR ;

gamma = 4 ; %lock number

V_knots = 325; %the rotor forward velocity [knots]

% Convert velocity from [knots] to [meters per second]
% 1 knot = 0.51444 meters per second
V = V_knots * 0.51444 ;

% Calculate angular velocity in radians per second as μ = V/Ω ∗ R = 1
omega_rad_s =   V / R ;

% Convert angular velocity from radians per second to RPM
% 1 radian per second = (60 / (2 * pi)) RPM
Omega = omega_rad_s * (60 / (2 * pi)) ;

%%%%%%%%% Aerodynamic coefficients calculations
M_u =(0.25*sqrt(2) - 0.25*log(1 + sqrt(2)));
M_b = (0.0625*sqrt(2) - 0.1875*log(1 + sqrt(2)));
H_u = 0.5*log(1 + sqrt(2));

% Frequency ranges
f_pitch= 0.05:0.1:0.45; % Frequency normalized [Cycle/s or Hz]
f_yaw= 0.05:0.1:0.45; % Frequency normalized [Cycle/s or Hz]

divergence_map = [];
Rdivergence_map = [];
unstable = [];

% Modify the loop to iterate over frequency points
for kk1 = 1:length(f_pitch)
    for kk2 = 1:length(f_yaw)

% Angular frequencies for the current frequency points
w_omega_pitch = 2 * pi * f_pitch(kk1);
w_omega_yaw = 2 * pi * f_yaw(kk2);

K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;

tspan = [ 0.0   pi/omega_rad_s^2] ;
 
[ t1 , y1 ] = ode45(@(t, y) two_bladed_system( t , y , omega_rad_s, I_thetaoverI_b, I_psioverI_b, K_theta, K_psi, h, gamma, H_u, M_b, M_u) , tspan , [1; 0; 0; 0]) ;
[ t2 , y2 ] = ode45(@(t, y) two_bladed_system( t , y , omega_rad_s, I_thetaoverI_b, I_psioverI_b, K_theta, K_psi, h, gamma, H_u, M_b, M_u) , tspan , [0; 1; 0; 0]) ;
[ t3 , y3 ] = ode45(@(t, y) two_bladed_system( t , y , omega_rad_s, I_thetaoverI_b, I_psioverI_b, K_theta, K_psi, h, gamma, H_u, M_b, M_u) , tspan , [0; 0; 1; 0]) ;
[ t4 , y4 ] = ode45(@(t, y) two_bladed_system( t , y , omega_rad_s, I_thetaoverI_b, I_psioverI_b, K_theta, K_psi, h, gamma, H_u, M_b, M_u) , tspan , [0; 0; 0; 1]) ;

[ nrows1 , ncols1 ] = size(y1) ;
[ nrows2 , ncols2 ] = size(y2) ;
[ nrows3 , ncols3 ] = size(y3) ;
[ nrows4 , ncols4 ] = size(y4) ;

% computation of monodromy matrix
MOO = [y1(nrows1 ,:)', y2(nrows2,:)', y3(nrows3,:)', y4(nrows4,:)'];

% computation of floquet multipliers
eigenvalues= eig(MOO);
EigVals{kk1,kk2} = eigenvalues;

 
            end 
        end


% Display eigenvalues for each frequency combination
for kk1 = 1:length(f_pitch)
    for kk2 = 1:length(f_yaw)
        fprintf('Eigenvalues for f_pitch = %f Hz and f_yaw = %f Hz:\n', f_pitch(kk1), f_yaw(kk2));
        disp(EigVals{kk1, kk2});
    end
end

Eigenvalues for f_pitch = 0.050000 Hz and f_yaw = 0.050000 Hz:

   1.0000 + 0.0000i
   1.0000 - 0.0000i
   1.0000 + 0.0000i
   1.0000 - 0.0000i

Eigenvalues for f_pitch = 0.050000 Hz and f_yaw = 0.150000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

Eigenvalues for f_pitch = 0.050000 Hz and f_yaw = 0.250000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

Eigenvalues for f_pitch = 0.050000 Hz and f_yaw = 0.350000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

Eigenvalues for f_pitch = 0.050000 Hz and f_yaw = 0.450000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

Eigenvalues for f_pitch = 0.150000 Hz and f_yaw = 0.050000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0000i
   1.0000 + 0.0000i

Eigenvalues for f_pitch = 0.150000 Hz and f_yaw = 0.150000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0000i
   1.0000 - 0.0000i

Eigenvalues for f_pitch = 0.150000 Hz and f_yaw = 0.250000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.150000 Hz and f_yaw = 0.350000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.150000 Hz and f_yaw = 0.450000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.250000 Hz and f_yaw = 0.050000 Hz:

   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.250000 Hz and f_yaw = 0.150000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.250000 Hz and f_yaw = 0.250000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.250000 Hz and f_yaw = 0.350000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.250000 Hz and f_yaw = 0.450000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.350000 Hz and f_yaw = 0.050000 Hz:

   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.350000 Hz and f_yaw = 0.150000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.350000 Hz and f_yaw = 0.250000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.350000 Hz and f_yaw = 0.350000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0001i
   1.0000 - 0.0001i

Eigenvalues for f_pitch = 0.350000 Hz and f_yaw = 0.450000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.450000 Hz and f_yaw = 0.050000 Hz:

   1.0000 + 0.0000i
   1.0000 + 0.0000i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.450000 Hz and f_yaw = 0.150000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.450000 Hz and f_yaw = 0.250000 Hz:

   1.0000 + 0.0001i
   1.0000 - 0.0001i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.450000 Hz and f_yaw = 0.350000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

Eigenvalues for f_pitch = 0.450000 Hz and f_yaw = 0.450000 Hz:

   1.0000 + 0.0002i
   1.0000 - 0.0002i
   1.0000 + 0.0002i
   1.0000 - 0.0002i

function [dy_dt] = two_bladed_system(t, y, omega_rad_s, I_thetaoverI_b, I_psioverI_b, K_theta, K_psi, h, gamma, H_u, M_b, M_u)

ct = cos(2*omega_rad_s^2*t);
st = sin(2*omega_rad_s^2*t);

    % Define inertia matrix [M]
    M11 = I_thetaoverI_b + 1 + cos(2*omega_rad_s^2*t);
    M12 = -sin(2*omega_rad_s^2*t);
    M21 = -sin(2*omega_rad_s^2*t);
    M22 = I_psioverI_b + 1 - cos(2*omega_rad_s^2*t);

    M = [M11 M12; M21 M22];

    % Define damping matrix [D]
    D11 =  h^2*gamma*H_u*(1 - cos(2*omega_rad_s^2*t)) - gamma*M_b*(1 + cos(2*omega_rad_s^2*t)) - (2 + 2*h*gamma*M_u)*sin(2*omega_rad_s^2*t);
    D12 =  h^2*gamma*H_u*sin(2*omega_rad_s^2*t) + gamma*M_b*sin(2*omega_rad_s^2*t) - 2*(1 + cos(2*omega_rad_s^2*t)) - 2*h*gamma*M_u*cos(2*omega_rad_s^2*t);
    D21 =  h^2*gamma*H_u*sin(2*omega_rad_s^2*t) + gamma*M_b*sin(2*omega_rad_s^2*t) + 2*(1 - cos(2*omega_rad_s^2*t)) - 2*h*gamma*M_u*cos(2*omega_rad_s^2*t);
    D22 =  h^2*gamma*H_u*(1 + cos(2*omega_rad_s^2*t)) - gamma*M_b*(1 - cos(2*omega_rad_s^2*t)) + (2 + 2*h*gamma*M_u)*sin(2*omega_rad_s^2*t);

    D = [D11 D12; D21 D22];

    % Define stiffness matrix [K]
    K11 = K_theta - h*gamma*H_u*(1 - cos(2*omega_rad_s^2*t)) + gamma*M_u*sin(2*omega_rad_s^2*t);
    K12 = -h*gamma*H_u*sin(2*omega_rad_s^2*t) + gamma*M_u*(1 + cos(2*omega_rad_s^2*t));
    K21 = -h*gamma*H_u*sin(2*omega_rad_s^2*t) - gamma*M_u*(1 - cos(2*omega_rad_s^2*t));
    K22 = K_psi - h*gamma*H_u*(1 + cos(2*omega_rad_s^2*t)) - gamma*M_u*sin(2*omega_rad_s^2*t);

    K = [K11 K12; K21 K22];

    % Calculate the inverse matrix-vector product and assign correctly
    acceleration = -inv(M) * (D * [y(2); y(4)] + K * [y(1); y(3)]);

    % System equations
    dy_dt = [y(2) ; acceleration(1); y(4); acceleration(2)] ;

end

Answer 1

[Prathamesh]

Hi @Nikoo,

I understand that the customer have a four-degree-of-freedom system. But the eigenvalues are very similar to each other and their values are extremely small.

Correct the Angular Frequency Term in “two_bladed_system” Function:

• The argument for the “cos” and sin functions, currently uses “2*omega_rad_s^2*t”. Based on standard dimensionless analysis this should be dimensionally consistent with an angle which is “2*omega_rad_s*t”

Change

ct = cos(2*omega_rad_s^2*t);    to       ct = cos(2*omega_rad_s*t); 
st = sin(2*omega_rad_s^2*t);    to       st = sin(2*omega_rad_s*t);

Correct the Time Span for “ode45” Integration:

• For Floquet analysis, “ode45” must integrate over exactly one period of the system's periodic coefficients. Your current tspan uses “pi/omega_rad_s^2”.

Change

tspan = [0.0 pi/omega_rad_s^2];     to           tspan = [ 0.0 pi/omega_rad_s ] ;

Please refer below documentation link for “ode45”

https://www.mathworks.com/help/matlab/ref/ode45.html