code says
V = 1000; % the rotor forward velocity [knots]
Where is this propellor? That is about 1.5 Mach at sea level. Fast!
The Mathieu Equation—Stability for 2DOF whirlflutter system
15 views (last 30 days)
Show older comments
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;
4 Comments
Umar
on 29 Jul 2024
Hi @Nikko,
Please see my response to your comments below.
You asked, “Is my approach for determining the Stiffnesses and angles correct”
Based on the analysis of your code, it seems that your approach for determining the stiffnesses and angles is correct. You correctly calculate the stiffness matrices at each time point and use them to determine the stability regions (flutter and divergence) of the system. However, it's important to note that without knowing the specific requirements and constraints of your problem, it's difficult to determine if your approach is entirely correct. It's always a good idea to validate your results against known analytical solutions or experimental data if available.
Now, addressing your next question about “how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly”
Looking at the provided code, it seems that your goal is to analyze the stability of a rotor system with multiple blades. The parameters and equations used in the code are related to the dynamics of the rotor system. So, to apply the Floquet technique and the Mathieu equation, you need to make the following modifications to the code by defining the Floquet parameters since they are related to the periodicity of the system. In this case, the periodicity is determined by the time periods for pitch and yaw, which are defined as T_pitch and T_yaw respectively. These parameters can be adjusted to change the periodicity of the system. Also, calculate the Floquet exponents since they are the eigenvalues of the Floquet matrix, which is derived from the Mathieu equation. In your code, the Floquet exponents are calculated using the roots function. To change the process, you can modify the equations used to calculate the matrices M_matrix, D_matrix, and K_matrix. These matrices represent the inertia, damping, and stiffness of the system respectively. By changing these equations, you can modify the behavior of the system and analyze its stability using the Floquet exponents. Finally, analyze the stability the stability of the system which is analyzed by checking the real and imaginary parts of the Floquet exponents. If the real part is positive and the imaginary part is non-positive, it indicates instability. You can modify the conditions for instability based on your specific requirements. Remember to test and validate the changes to ensure they produce the desired results.
Good luck!
Answers (1)
William Rose
on 30 Jul 2024
Edited: William Rose
on 30 Jul 2024
[Edit: Replaced scan of Figure 9 with a better scan that does not cut off the Y axis label.]
You seek advice for code to reproduce the figure in your original post. That figure is Figure 9 on p.173 of NASA Technical Note D-7677, 1974, here. Also available here https://ntrs.nasa.gov/api/citations/19740019387/downloads/19740019387.pdf. (Page numbers are the page numbers printed in the document, which may not be the same as the page numbers of the scan.)
When I run your code, I get the figure below. Figure 9 from the Technical Note is next to it.
This figure has no axis labels and no legend. Please provide both. The range on the horizontal and vertical axes is 0 to 6 x 10^4. The axis ranges on Figure 9 are 0 to 15. Therefore I recommend that you resolve the discrepancy between the axes of your figure and the figure you are trying to replicate.
The Tech Note says (p.93) "a 1/rev divergence occurs where 
or 
is near 1/rev..." You identify regions that correspond to 1/Rev Divergence (labelled "1" in Figure 9) as regions where abs(real(r(k))-freq_1_over_Omega) < 1e-5. Why did you choose 1e-5 as your criterion for "near 1/rev"? Why not use 0.1 or some other number?
or is near 1/rev..." You identify regions that correspond to 1/Rev Divergence (labelled "1" in Figure 9) as regions where abs(real(r(k))-freq_1_over_Omega) < 1e-5. Why did you choose 1e-5 as your criterion for "near 1/rev"? Why not use 0.1 or some other number?Since time is measured with the dimensionless variable ψ (see p.xi), does that mean the Ω=rotor rotational speed should always equal 
radians per unit of dimensionless time? The value for Omega in your code is Omega=V/R which yields a value of 1133.
radians per unit of dimensionless time? The value for Omega in your code is Omega=V/R which yields a value of 1133. Your code comments say that the units for V is knots and units for Omega is RPM. What are the units for rotor blade radius R in your code? How does [knots]/[R]=[RPM]? I suggest you re-evaluate your equation for Omega.
The equations related to Figure 9 are described on pp.90-93 of Technical Note. Regarding Figure 9, the Technical Note says on p.93:
Equation 147 on p.93 is a set of two coupled second order linear differential equations for 
=rotor shaft yaw angle at pivot, and 
=rotor shaft pitch angle at pivot (defined on p.ix), that vary periodically. The derivatives in equation 147 are taken with respect to the rotor blade angle (see p.xi and p.xii). This, and the fact that the rotor under consideration in Figure 9 has two blades, explains why the coefficients have period=pi. The Technical Note says on p.v "quantities are made dimensionless with ρ, Ω, and R (air density, rotor rotational speed, and rotor radius)".
=rotor shaft yaw angle at pivot, and =rotor shaft pitch angle at pivot (defined on p.ix), that vary periodically. The derivatives in equation 147 are taken with respect to the rotor blade angle (see p.xi and p.xii). This, and the fact that the rotor under consideration in Figure 9 has two blades, explains why the coefficients have period=pi. The Technical Note says on p.v "quantities are made dimensionless with ρ, Ω, and R (air density, rotor rotational speed, and rotor radius)".Equation 147 is similar to Mathieu's equation, but more complicated. I do not know if your code correctly applies Floquet theory to find the regions of stability for the system described by equation 147. I do realize that Floquet theory is useful for characterizing the stability of a system described by Mathieu's equation.
In your comment above, you posted p.100 of the same NASA Technical Note. It says a typical value for V is 325 knots. Your V is >3x larger. Since you want your results to match the figure, and they don't match now, you might as well use the value for V that was used to make the figure.
5 Comments
Umar
on 31 Jul 2024
Hi @Nikko,
I made modifications to your code and experimenting with it. However, I am too tired now, please take a look at the code below to give you some clues to find solution to your problem, when you execute code, it will display “Index in position 2 exceeds array bounds”. However, pay close attention to plot. You have to change parameters in code to find flaw in your code and fix the problem. At the moment, I am too tired. Hope to hear good news.
% 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));
% Frequency ranges
f_pitch = 0.01:5:80;
f_yaw = 0.01:5:80;
% 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 = [];
% Define the value of H_u
H_u = 1; % Replace with the actual value
% Define the values of T_pitch and T_yaw
T_pitch = [1, 2, 3]; % Replace with the actual values
T_yaw = [4, 5, 6]; % Replace with the actual values
% Define the values of I_psioverI_b, I_thetaoverI_b,
C_thetaoverI_b, C_psioverI_b, M_b, M_u, gamma, h, and V
I_psioverI_b = 0.5; % Replace with the actual value I_thetaoverI_b = 0.7; % Replace with the actual value C_thetaoverI_b = 0.3; % Replace with the actual value C_psioverI_b = 0.2; % Replace with the actual value M_b = 0.9; % Replace with the actual value M_u = 0.8; % Replace with the actual value gamma = 0.6; % Replace with the actual value h = 0.1; % Replace with the actual value V = 0.5; % Replace with the actual value
% Initialize the matrices to store the results
unstable = [];
divergence_map = [];
Rdivergence_map = [];
% 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 = 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 = 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)
if size(unstable, 1) < k
unstable = [unstable; zeros(k - size(unstable, 1), 3)];
end
unstable(k, :) = [phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
freq_1_over_Omega = 0.1; % Replace with the actual value
if abs(real(r(k)) - freq_1_over_Omega) < 0.1
if size(Rdivergence_map, 1) < k
Rdivergence_map = [Rdivergence_map; zeros(k -
size(Rdivergence_map, 1), 3)];
end
Rdivergence_map(k, :) = [phi, K_psi, K_theta];
end
end
end
end
% Divergence
if (real(det(K_matrix)) < 0)
if size(divergence_map, 1) < k
divergence_map = [divergence_map; zeros(k -
size(divergence_map, 1), 3)];
end
divergence_map(k, :) = [phi, K_psi, K_theta];
end
end
end
end
% Plotting
figure;
scatter3(unstable(:,1), unstable(:,2), unstable(:,3), 'r',
'filled');
hold on;
scatter3(Rdivergence_map(:,1), Rdivergence_map(:,2),
Rdivergence_map(:,3), 'g', 'filled');
hold on;
scatter3(divergence_map(:,1), divergence_map(:,2),
divergence_map(:,3), 'b', 'filled');
xlabel('Phi');
ylabel('K_psi');
zlabel('K_theta');
legend('Unstable', '1/Omega Divergence', 'Divergence');
title('Flutter Analysis Results');
Please see attached plot.
Good luck!
William Rose
on 31 Jul 2024
I have looked at various web pages related to Floquet theory including wikipedia and multiple others. This undergraduate thesis by E. Folkers looks more understandable to me than other sources.
In particular, see
Section 3.3 Stability of the Floquet system
Example 3.11 (system matrix A is 2x2)
Section 3.4 Hill's differential equation (system matrix A is 2x2)
Maybe the above will help you. I do not have time to really understand it.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)