Analyze the eigenvalues of a system with periodic coefficients.

9 visualizzazioni (ultimi 30 giorni)
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

Risposte (1)

Prathamesh
Prathamesh il 4 Giu 2025
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”

Categorie

Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by