How to set the theorem of floquet for Linear variational equation with periodic equation?
I’m trying to analyze the stability regions of the equation of motion of a system with periodic coefficients (whirlflutter). I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.
If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?
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 [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]
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
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));
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= 5:3:140;
f_yaw= 5:3:140;
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
phi_steps = linspace(0, pi, 100); % Time steps within one period
for phi = phi_steps
% Angular frequencies for the current time point
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 = C_thetaoverI_b + 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 = C_psioverI_b + 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);
% Find the roots of the polynomial equation
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) < 0.1
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; phi, 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’m trying to analyze the stability regions of the equation of motion of a system with periodic coefficients (whirlflutter). I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.
If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?
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 [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]
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
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));
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= 5:3:140;
f_yaw= 5:3:140;
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
phi_steps = linspace(0, pi, 100); % Time steps within one period
for phi = phi_steps
% Angular frequencies for the current time point
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 = C_thetaoverI_b + 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 = C_psioverI_b + 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);
% Find the roots of the polynomial equation
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) < 0.1
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; phi, 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’m trying to analyze the stability regions of the equation of motion of a system with periodic coefficients (whirlflutter). I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.
If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?
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 [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]
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
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));
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= 5:3:140;
f_yaw= 5:3:140;
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
phi_steps = linspace(0, pi, 100); % Time steps within one period
for phi = phi_steps
% Angular frequencies for the current time point
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 = C_thetaoverI_b + 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 = C_psioverI_b + 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);
% Find the roots of the polynomial equation
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) < 0.1
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; phi, 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; stability theory, floquet theorem, characteristic multiplier, periodic coeffcients MATLAB Answers — New Questions
