The Mathieu Equation—Stability for 2DOF whirlflutter system
I am attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. As my coefficients are periodic
Is my approach for determining the Stiffnesses and angles correct, and how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly according to the reference.
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; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;
gamma = 4; % lock number
V = 1000; % 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;
% Frequency ranges
f_pitch= 0.001:0.5:5;
f_yaw= 0.001:0.5:5;
% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
for j = 1:length(T_yaw)
T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
t_steps = linspace(0, T, 100); % Time steps within one period
for t = t_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi / T_pitch(i);
w_omega_yaw = 2 * pi / T_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Calculate matrices at time t using harmonic motion expressions
phi = 2 * pi * t / T; % Phase variation over the period
% 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];
% Compute the system matrices
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);
P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;
P = roots([P0, P1, P2, P3, P4]);
r = 1 * P;
%Flutter
for k = 1:length(r)
if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 1e-5
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; t, K_psi, K_theta];
end
end
end
end
% Plotting section
figure;
hold on;
scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);
xlabel(‘K_psi’);
ylabel(‘K_theta’);
title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);
hold off;I am attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. As my coefficients are periodic
Is my approach for determining the Stiffnesses and angles correct, and how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly according to the reference.
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; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;
gamma = 4; % lock number
V = 1000; % 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;
% Frequency ranges
f_pitch= 0.001:0.5:5;
f_yaw= 0.001:0.5:5;
% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
for j = 1:length(T_yaw)
T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
t_steps = linspace(0, T, 100); % Time steps within one period
for t = t_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi / T_pitch(i);
w_omega_yaw = 2 * pi / T_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Calculate matrices at time t using harmonic motion expressions
phi = 2 * pi * t / T; % Phase variation over the period
% 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];
% Compute the system matrices
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);
P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;
P = roots([P0, P1, P2, P3, P4]);
r = 1 * P;
%Flutter
for k = 1:length(r)
if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 1e-5
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; t, K_psi, K_theta];
end
end
end
end
% Plotting section
figure;
hold on;
scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);
xlabel(‘K_psi’);
ylabel(‘K_theta’);
title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);
hold off; I am attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. As my coefficients are periodic
Is my approach for determining the Stiffnesses and angles correct, and how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly according to the reference.
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; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;
gamma = 4; % lock number
V = 1000; % 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;
% Frequency ranges
f_pitch= 0.001:0.5:5;
f_yaw= 0.001:0.5:5;
% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
for j = 1:length(T_yaw)
T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
t_steps = linspace(0, T, 100); % Time steps within one period
for t = t_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi / T_pitch(i);
w_omega_yaw = 2 * pi / T_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Calculate matrices at time t using harmonic motion expressions
phi = 2 * pi * t / T; % Phase variation over the period
% 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];
% Compute the system matrices
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);
P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;
P = roots([P0, P1, P2, P3, P4]);
r = 1 * P;
%Flutter
for k = 1:length(r)
if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 1e-5
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; t, K_psi, K_theta];
end
end
end
end
% Plotting section
figure;
hold on;
scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);
xlabel(‘K_psi’);
ylabel(‘K_theta’);
title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);
hold off; floquet technic, dynamical system MATLAB Answers — New Questions
