基于Floquet理论的2DOF周期系统涡动颤振图绘制问题求助

我采用Floquet理论实现带周期系数的2自由度（2DOF）系统稳定性边界，试图在Matlab中通过特征值识别颤振、发散等不稳定区域并绘制涡动颤振图（whirl flutter diagram），但所得图表与参考结果不符。我的计算步骤为：定义系统周期$T=1/(2\pi f)$，遍历周期内不同时间点$t$，构造特征多项式求解特征值分析根的特性，目前对Floquet理论的应用方式存疑，附上Matlab代码如下：

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.01:5:140;
f_yaw= 0.01:5:140;

% 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; t, K_psi, K_theta];

               % Proximity check for 1/Ω divergence
               if abs(real(r(k)) - freq_1_over_Omega) < 1e-5
                  Rdivergence_map = [Rdivergence_map; t, 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;

核心问题与修正方向 #

• Floquet理论应用错误：当前直接对每个时间点的时变矩阵求解特征值，并非Floquet理论的正确用法。Floquet分析需计算单值矩阵（Monodromy Matrix）——系统状态矩阵在一个完整周期内的转移矩阵，再通过单值矩阵的特征值（Floquet乘子）判断稳定性：若任意乘子的模大于1，系统不稳定。
• 周期定义错误：系统的周期由转子旋转频率Ω决定，应为$T=2\pi/\Omega$，而非俯仰/偏航频率的最大周期。你的系统是因转子旋转产生周期系数，周期应与转子转速绑定。
• 状态空间建模缺失：Floquet分析要求将二阶系统转化为一阶状态空间形式（如令$x=[\theta, \psi, \dot{\theta}, \dot{\psi}]^T$），再求解状态转移矩阵的周期演化。
• 不稳定判据逻辑错误：
  • 颤振判据：应基于Floquet乘子的模，而非某时刻特征值的实部
  • 发散判据：不能仅通过det(K_matrix)<0判断，需结合Floquet分析或平均刚度矩阵的正定性分析
• 变量定义混淆：K_psi和K_theta的赋值逻辑颠倒，且应对应系统的结构刚度参数，而非由俯仰/偏航频率推导得出。

关键修正步骤示例 #

1. 构建一阶状态空间模型
   将二阶系统$M(t)\ddot{q} + D(t)\dot{q} + K(t)q = 0$转化为一阶状态方程：
   $$\dot{x} = A(t)x$$
   其中$x=[q_1, q_2, \dot{q_1}, \dot{q_2}]^T$，$A(t)$为4×4时变矩阵：
   $$A(t) = \begin{bmatrix}
   0 & 0 & 1 & 0 \
   0 & 0 & 0 & 1 \
   -M^{-1}(t)K(t){1,:} & -M^{-1}(t)D(t){1,:} \
   -M^{-1}(t)K(t){2,:} & -M^{-1}(t)D(t){2,:}
   \end{bmatrix}$$
2. 计算单值矩阵
   使用Matlab的ode45求解状态转移矩阵的周期演化：从$t=0$到$t=T$，初始转移矩阵为单位矩阵，求解$\dot{\Phi}(t) = A(t)\Phi(t)$，得到$\Phi(T)$即为单值矩阵。
3. 求解Floquet乘子并判断稳定性
   计算单值矩阵的特征值，若任何特征值的模大于1，则标记当前参数组合为不稳定。
4. 遍历参数绘制颤振图
   遍历不同的结构刚度（$K_\psi, K_\theta$），对每组参数执行上述分析，标记不稳定区域并绘图。

内容的提问来源于stack exchange，提问作者Ookin Etemadi