The Mathieu Equation—Stability for 2DOF whirlflutter system

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!