indicating regions of flutter (F), and divergence (D)

I tried to implemented the stability boundries of a 2DOF system while my coefficents are periodic.
Initially, I wrote my homogeneous matrix and directly examined the eigenvalues, but the final plot didn’t match the reference.
clc
clear all ;

% Define parameters

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; % rotor mast height
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 325 ; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

%%%%%%%%%%%the flap moment aerodynamic coefficients for small V
%M_b = -1*(1 + V^2)/8 ;
%M_u = V/4;
%the propeller aerodynamic coefficients
%H_u = (V^2/2)*log(2/V);

f_pitch= 0.01:5:140;
f_yaw= 0.01:5:140;

phi_steps = linspace(0, pi, 50); % Evaluation points from 0 to pi

divergence_map = [];

Rdivergence_map = [];

unstable = [];

for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
for phi = phi_steps

% Calculate stiffness for the current frequency
w_omega_pitch = 2*pi*f_pitch(i);
w_omega_yaw = 2*pi*f_yaw(j);

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

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];

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

D_matrix = [D11, D12;
D21, D22];

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

K_matrix = [K11, K12;
K21, K22];

A_top = [zeros(2, 2), eye(2)];
A_bottom = [-inv(M_matrix) * K_matrix, -inv(M_matrix) * D_matrix];
A = [A_top; A_bottom];

eigenvalues = eig(A);

% Stability condition
% Flutter

if any(real(eigenvalues) > 0)
unstable = [unstable; K_psi, K_theta];
end

% Divergence condition
if det(K_matrix) < 0
divergence_map = [divergence_map; K_psi, K_theta];
end

% 1/Ω *(Divergence) proximity check
for ev = eigenvalues’
if abs(ev – freq_1_over_Omega) < 1e-2
Rdivergence_map = [Rdivergence_map; K_psi, K_theta];
end
end
end
end
end

% Plot the Flutter and divergence maps
figure;
hold on;

scatter(unstable(:,1), unstable(:,2), ‘filled’);
scatter(divergence_map(:,1), divergence_map(:,2), ‘filled’, ‘r’);
scatter(Rdivergence_map(:, 1), Rdivergence_map(:, 2), ‘filled’, ‘g’);

xlabel(‘K_psi’);
ylabel(‘K_theta’);

title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘ 1/Ω Divergence area’);

hold off;

I’d appreciate it if you could check it. Also, if I need to use the Floquet technique, please let me know how to modify the code.I tried to implemented the stability boundries of a 2DOF system while my coefficents are periodic.
Initially, I wrote my homogeneous matrix and directly examined the eigenvalues, but the final plot didn’t match the reference.
clc
clear all ;

% Define parameters

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; % rotor mast height
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 325 ; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

%%%%%%%%%%%the flap moment aerodynamic coefficients for small V
%M_b = -1*(1 + V^2)/8 ;
%M_u = V/4;
%the propeller aerodynamic coefficients
%H_u = (V^2/2)*log(2/V);

f_pitch= 0.01:5:140;
f_yaw= 0.01:5:140;

phi_steps = linspace(0, pi, 50); % Evaluation points from 0 to pi

divergence_map = [];

Rdivergence_map = [];

unstable = [];

for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
for phi = phi_steps

% Calculate stiffness for the current frequency
w_omega_pitch = 2*pi*f_pitch(i);
w_omega_yaw = 2*pi*f_yaw(j);

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

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];

D_matrix = [D11, D12;
D21, D22];

K_matrix = [K11, K12;
K21, K22];

A_top = [zeros(2, 2), eye(2)];
A_bottom = [-inv(M_matrix) * K_matrix, -inv(M_matrix) * D_matrix];
A = [A_top; A_bottom];

eigenvalues = eig(A);

% Stability condition
% Flutter

if any(real(eigenvalues) > 0)
unstable = [unstable; K_psi, K_theta];
end

% Divergence condition
if det(K_matrix) < 0
divergence_map = [divergence_map; K_psi, K_theta];
end

% 1/Ω *(Divergence) proximity check
for ev = eigenvalues’
if abs(ev – freq_1_over_Omega) < 1e-2
Rdivergence_map = [Rdivergence_map; K_psi, K_theta];
end
end
end
end
end

% Plot the Flutter and divergence maps
figure;
hold on;

xlabel(‘K_psi’);
ylabel(‘K_theta’);

title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘ 1/Ω Divergence area’);

hold off;

I’d appreciate it if you could check it. Also, if I need to use the Floquet technique, please let me know how to modify the code. I tried to implemented the stability boundries of a 2DOF system while my coefficents are periodic.
Initially, I wrote my homogeneous matrix and directly examined the eigenvalues, but the final plot didn’t match the reference.
clc
clear all ;

% Define parameters

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; % rotor mast height
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 325 ; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

%%%%%%%%%%%the flap moment aerodynamic coefficients for small V
%M_b = -1*(1 + V^2)/8 ;
%M_u = V/4;
%the propeller aerodynamic coefficients
%H_u = (V^2/2)*log(2/V);

f_pitch= 0.01:5:140;
f_yaw= 0.01:5:140;

phi_steps = linspace(0, pi, 50); % Evaluation points from 0 to pi

divergence_map = [];

Rdivergence_map = [];

unstable = [];

for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
for phi = phi_steps

% Calculate stiffness for the current frequency
w_omega_pitch = 2*pi*f_pitch(i);
w_omega_yaw = 2*pi*f_yaw(j);

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

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];

D_matrix = [D11, D12;
D21, D22];

K_matrix = [K11, K12;
K21, K22];

A_top = [zeros(2, 2), eye(2)];
A_bottom = [-inv(M_matrix) * K_matrix, -inv(M_matrix) * D_matrix];
A = [A_top; A_bottom];

eigenvalues = eig(A);

% Stability condition
% Flutter

if any(real(eigenvalues) > 0)
unstable = [unstable; K_psi, K_theta];
end

% Divergence condition
if det(K_matrix) < 0
divergence_map = [divergence_map; K_psi, K_theta];
end

% 1/Ω *(Divergence) proximity check
for ev = eigenvalues’
if abs(ev – freq_1_over_Omega) < 1e-2
Rdivergence_map = [Rdivergence_map; K_psi, K_theta];
end
end
end
end
end

% Plot the Flutter and divergence maps
figure;
hold on;

xlabel(‘K_psi’);
ylabel(‘K_theta’);

title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘ 1/Ω Divergence area’);

hold off;

I’d appreciate it if you could check it. Also, if I need to use the Floquet technique, please let me know how to modify the code. aeroelastic, floquet, flutter MATLAB Answers — New Questions

Cart

Cart