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;

[William Rose]

@Nikoo,

code says

V = 1000; % the rotor forward velocity [knots]

Where is this propellor? That is about 1.5 Mach at sea level. Fast!

[Nikoo]

the parameters has been taken from this table but I was changing the inputs to analyze the behavior of diagram.

[Nikoo]

@William Rose The inputs dont seem very suspicious to me. my questions are more about the implementation process and approach for these boundries. However, the above document reviews all the propeller constants.

[Umar]

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!

Answer 1

[William Rose]

@Nikoo,

[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?

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.

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)".

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.

[Nikoo]

Thank you, @William Rose and @Umar, for your hints.

I tried to implement them.

First, I increased the frequency range to show from around zero to 15.

Second, I changed the calculation method for omega (rotor rotational speed), considering that the speed unit was in knots and the radius in meters.

Third, for finding the unstable/divergence regions, I wrote down the polynomial equation and then attempted to examine the roots.

Regarding the 1/Ω divergence, I considered the value to be close to zero and 0.1.

I have three different conditions, but the regions are the same. Do you see any major issues in my updated approach? However, the boundaries are not being correctly identified.

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= 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 = [];

% 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)
               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;

[Umar]

Hi @Nikko,

Based on the detailed steps you provided for your updated approach in analyzing the rotor system dynamics, it seems like you have made significant modifications to the frequency range, calculation method for rotor rotational speed, and examination of unstable/divergence regions. Let's break down your approach and address any potential major issues.

Frequency Range: Increasing the frequency range is a good step to ensure that you cover a wider range of dynamics. However, it's important to make sure that the range is appropriate for your specific problem. You mentioned that you increased the range to show from zero to 15, but it's unclear if this range is suitable for your system. Make sure to choose a range that captures the relevant frequencies for your problem.

Calculation of Omega: You mentioned that you changed the calculation method for omega (rotor rotational speed) to consider the speed unit in knots and the radius in meters. It's important to ensure that the conversion is done correctly. From the code provided, it seems that you have correctly converted the velocity from knots to meters per second and then calculated the angular velocity in radians per second. However, it's always a good idea to double-check the conversion formulas and units to avoid any errors.

Examination of Roots: You mentioned that you attempted to examine the roots of the polynomial equation to find the unstable and divergence regions. From the code provided, it seems that you have correctly calculated the roots using the roots function. However, it's important to note that the roots of a polynomial equation can be complex numbers. In your code, you are checking if the real part of the root is greater than zero and the imaginary part is less than or equal to zero to identify unstable regions. This approach seems reasonable, but it's important to consider all possible cases and ensure that you are correctly identifying the unstable regions.

Boundary Identification: You mentioned that the boundaries are not being correctly identified. It's difficult to pinpoint the exact issue without more information or specific error messages. However, it's possible that the issue lies in the logic used to identify the boundaries. Double-check the conditions and proximity checks you are using to determine the boundaries and ensure that they are correctly implemented.

To further investigate the issues and provide more specific guidance, it would be helpful to have more information about the specific problem you are trying to solve and any error messages or unexpected results you are encountering.

To provide example code that correctly identifies boundaries, you can refine your approach by double-checking the accuracy of matrix calculations, root analysis, and proximity checks. Here is an example snippet for plotting:

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;

By refining these aspects and ensuring consistency in your calculations and analysis, you can enhance the accuracy and effectiveness of your updated approach in identifying unstable regions in the rotor system dynamics. Now let us wait what @William Rose has to say about your recent modifications made to the code.

[William Rose]

@Nikoo,

Can you add to your posting above the plot that is generated by your code? Just run the code in the window. Please include axis labels for the plot, with units, if the units are not dimensionless. Figre 9 in NASA Tech Note D-7677 has Kx* and Ky* on the axes. You have K_psi and K_theta. Does your K_psi equal the Tech Note's Kx*? Does your K_theta equal the Tech Note's Ky*? Read the Tech Note carefully to be sure.

I have never used Floquet theory to analyze a linear system with coefficients that vary periodically. I will have to do some reading and thinking to give an opinion on whether the approach you use above is reasonable.

Here is my description of the algorithm in your code, as it appears to me, from reading your code.

___________________________

An outer pair of loops iterates over f_pitch and f_yaw, which are defined as

f_pitch= 0.01:5:80;
f_yaw= 0.01:5:80;

Then there is an innermost loop, a third "for" loop. It loops over t: 100 time steps, which make one complete one cycle of whichever is slower, f_pitch or f_yaw. Inside each pass of the loop over t, w_omega_pitch, w_omega_yaw, K_psi, and K_theta are computed. They do not depend on t, so they could be computed outside the loop for efficiency. The blade angle phi at time t is also computed. Then three 2x2 matrices are computed: inertia matrix M, damping matrix D, and stiffness matrix K. Each matrix depends on the current blade angle, phi. The elements of matrices M, D, and K are used to construct a 4th order polynomial. The roots of this polynomial are found. Each root is analyzed. If root r(k) has re(r(k))>0 and imag(r(k))<=0, the root is labeled unstable, and a row is added to array "unstable", with the K_psi and K_theta values for the current pass.

Roots that are found to be unstable are also tested for being possibly in the 1/Rev divergence zone. The test used is abs(real(r(k))-freq_1_over_Omega)<0.1. If the test is true, then a row is added to array "Rdivergence", with the current values of K_psi and K_theta.

Then stiffness matrix K is analyzed at this time step, t. If the real part of the determinant of K is <0, then a row is added to array "divergence_map", with the current values of K_psi and K_theta.

After all roots for the current time step (and associated blade angle) are analyzed, the t loop increments by 1, and the analysis is repeated, with the blade angle, phi, further around.

After all 100 time steps have been analyzed for the current f_pitch and f_yaw, the next f_pitch, f_yaw combination is analyzed.

____________________

Here are some thoughts on the algorithm described above. If my description of the algorithm above is wrong, then my thoughts that follow are probably also wrong.

• I like how you allocate empty arrays for the different regimes, and then you add a row to the appropriate array each time the eigenvalues turn out a certain way. Then, after all points in the regime of interest are analyzed, you use scatter to plot the different arrays. That is a nice way of doing it.
• I think you should have an array, initially empty, called "stable". If a point is not unstable and not divergent, then a row should be added to the "stable " array, with the cirrent K_psi and K_theta. The stable points should be plotted with a distinct symbol and color, on the final plot.
• What are the units for frequencies f_pitch and f_yaw? As I mentioned above, the Tech Note uses dimensionless time for the analysis. Do the ranges for f_pitch and f_yaw correspond to Kx*=0 to 15, and Ky*=0 to 15? If not, why not? After all, you are trying to reproduce Figure 9.
• Since there are 4 roots to a 4th order polynomial, it is possible that, for a single time t, and associated blade angle phi, some roots are unstable, and some roots are stable. But that does not make sense, because, at a given value of f_yaw and f_pitch and t, the system is either unstable or not. The system cannot be stable and unstable at the same time. But your analysis allows for such a non-sensical result. I think the solution to this is to do the stability analysis by looking at all the roots. If any root is unstable, then the system is unstable for that combination of parameters.
• I have not evaluated your 4th order polynomial to see if it correctly represents the characteristic equation for the system.
• Even if you fix the issue related to analyzing each root separately, there is still another similar issue: You analyze the system at 100 time values, for each comination of f_pitch and f_yaw. K_phi and K_theta do not depend on time. They only depend on f_pitch and f_yaw, at least according to your formulas. It is possible with your algorithm, that there will be 50 unstable points at a particular f_pitch and f_yaw, and 50 unstable points. Or some other combination of 100 total stable and unstable points, for the same f_pitch and f_yaw. But that does not make sense, because, for a given f_pitch and f_yaw, the system should be either stable or unstable, not both. Likewise, the algorithm could find flutter at 50 values of t, and not flutter at the other 50 values of t. That too does not make sense, for the same reason. I am not sure how to fix this. I think the problem indicates that your implementation of Floquet theory to this system is incorrect. I do not know what the correct approach is. As I said above, I have never learned or applied Floquet theory.
• I have not verified the correctness of your formulas for K_psi and K_theta, as functions of f_pitch and f_yaw.
• I have not verified the correctness of your criterion for identifying divergence from analysis of the K matrix.
• I haved not verified the correctness of your criterion for identifying instability from anaysis of the roots of the characteristic polynomial.
• You identify unstable points as being in the 1/Rev divergence region by analyzing the real parts of the roots of the 4th order polynomial. Let's assume r(k)=a+ib, where i=sqrt(-1). (I assume that the polynomial corresponds to the characteristic equation, but I have not verified that.) I recommend that you consider the imaginary part, b, rather than the the real part, a. For a standard linear differential equaiton, the solution has an oscillatory component iff one or more roots are complex. If time is in seconds, then the frequency of oscillation (in rad/s) equals b, the imaginary component of the root. Equation 147 in the Tech Note is a differential equation with repsect to dimensionless time (equal to blade rotation angle). Therefore the frequecies of the solutions will also be dimensionless, or relative to 2 pi radians of blade rotation. Therefore I recommend that you analyze the imaginary part, b, for unstable points. If b equals approximately 1 or 2 for any of the roots, then there is a component of the solution that is close in frequency to the whole propellor (b=1) or every-blade (b=2) rotation frequency, and this means the system is in the 1/Rev divergence region.
• You identify points as being in the Divergence region if the real part of the determinant of the K (stiffness) matrix is <0. This suprises me, because, for many linear systems, an eigenvalue with a positive real part indicates instability (a solution that grows exponentially with time). I admit that I do not understand how the determinant of the stiffness matrix relates to instability, so maybe your method is fine. But you might want to check this.

I am too tired to think about this any more. I hope my analysis is helpful to you. Please remember that I could be wrong. You probably know much more about all this than I do.

[Umar]

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]

@Nikoo,

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.

http://www.12000.org/my_notes/liapunov_floquet_transformation/bMATH_2018_FolkersE.pdf

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.