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.