K法求解颤振方程-MATLAB实现

最新推荐文章于 2025-04-28 07:29:58 发布

气动弹性小菜鸟

最新推荐文章于 2025-04-28 07:29:58 发布

阅读量805

点赞数 9

文章标签： matlab 算法开发语言

本文链接：https://blog.csdn.net/wgh160613/article/details/142870825

版权

简介

颤振是飞行器在气流中运动发生的一种自激振动现象，是结构弹性变形、非定常气动力和惯性力组成的系统不稳定现象。当颤振发生时，飞行器通常在很短的时间就发生破坏，因此预测飞行器的颤振边界对确保飞行安全具有重要意义。工程上常用线性方法计算颤振边界，即采用频域非定常气动力理论、模态降阶法计算结构动力学响应并用线性方法求解颤振方程。这种方法可以快速求解颤振边界，对于非线性效应不强的工况能够准确计算。K法是一种简单易行的求解颤振方程的方法，本文简单介绍其基本原理并给出MATLAB实现。

理论

K法求颤振边界的基本方程如下

$[-\omega^2\boldsymbol{M}+(1+i\boldsymbol{g}_s)\boldsymbol{K}-q_\infty \boldsymbol{Q}_{ik})]q=0$

其中

$\omega$ ：系统频率
$i\boldsymbol{g}_s$ ：为使系统处于临界（简谐运动）状态需要增加的人工阻尼项
$q_\infty$ ：来流动压
$\boldsymbol{M}$ ：质量阵
$\boldsymbol{K}$ ：刚度阵
$\boldsymbol{Q}$ ：气动力影响系数（AIC）矩阵

定义减缩频率

$k=\frac{\omega L}{V}$

则上式可以写为

$[\boldsymbol{M}+\frac{\rho}{2}(\frac{L}{k})^2\boldsymbol{Q}(ik)- \lambda \boldsymbol{K}]\boldsymbol{q}=0$

其中

$\lambda = \frac{1+i\boldsymbol{g}_s}{\omega ^2}$

是上式的特征值。因此，颤振方程的求解问题本质上就是求上式的特征值。为了求解特征值，需要求解给定马赫数条件下一系列减缩频率 $k$ 的AIC矩阵，随后求出每个减缩频率对应的n个特征值，其中n为模态数。颤振频率、来流速度和人工阻尼可以通过下式求出

$\omega_f = \frac{1}{\mathrm{Re}(\lambda)}$

$g_s = \omega_f^2\mathrm{Im}(\lambda) = \frac{\mathrm{Im}(\lambda)}{\mathrm{Re}(\lambda)}$

$V_f=\frac{\omega_f L}{k}$

MATLAB代码

本节给出了K法的MATLAB实现代码，以ZAERO的Flutter Test Case1为测试算例，结构阻尼取0.01

%% Initialization
clear all;
clc;
load QHHSK.mat MatRes
QHHS = MatRes;
clear MatRes

%% Flutter solution with K method
% input model data
REFC = 2.07055;
REFL = REFC / 2;
M = diag([0.24855E-04 0.90881E-05 0.85232E-05 0.79439E-05]);
K = diag([0.11574E+01 0.15822E+02 0.22821E+02 0.12636E+03]);
damp = 0.01;
C = damp * K;

% input operating condition
Ma = 0.45;
rho = 1.073e-07;
Vlist = [5000, 5500, 6000, 6100, 6200, 6500, 6300, 6350, 6400];
klist = [0 0.05, 0.08, 0.10, 0.11, 0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.50];

NoModes = size(M,1);
Nok = length(klist);
RES = zeros(Nok,7,NoModes);
for jj = 1:Nok

    k = klist(jj);
    QHHStmp = QHHS(:,:,jj);
    
    Mbar = (k/REFL)^2*M + 0.5 * rho * QHHStmp;
    [V,D] = eig(K,Mbar);
    [Dlist,iD] = sort(diag(D));
    V = V(iD);

    for nn = 1:NoModes

        RlD = real(Dlist(nn));
        ImD = imag(Dlist(nn));

        if RlD >= 0
            g_s = -ImD/RlD;
            V_f = sqrt(RlD+ImD^2/RlD);
            omega_f = k*V_f/REFL;
        else
            omega_f = 99999999;
            g_s = 99999999;
            V_f = 99999999;
        end

        RES(jj,1,nn) = 1/k;
        RES(jj,2,nn) = V_f;
        RES(jj,3,nn) = g_s;
        RES(jj,4,nn) = 2*pi/omega_f;
        RES(jj,5,nn) = omega_f;
        RES(jj,6,nn) = RlD;
        RES(jj,7,nn) = ImD;

    end


end

%% Write output file
fid = fopen("ResKMethod.txt","w");
fprintf(fid,"     1/k         V_f         g_s                T         omega_f         Rld              Ild\n");
for nn = 1:NoModes
    fprintf(fid,"mode %d:\n",nn);
    for jj = Nok:-1:1
        fprintf(fid,"%8.2f  % 8.6E  % 8.6E  % 8.6E  % 8.6E  % 8.6E  % 8.6E  \n" ...
            ,RES(jj,1,nn),RES(jj,2,nn),RES(jj,3,nn),RES(jj,4,nn),RES(jj,5,nn),RES(jj,6,nn),RES(jj,7,nn));
    end
end
fclose(fid);

%% Plot V-g & V-f diagram
sid = 2;
% V-g diagram
figure(1)
namelist=[];
for i = 1:NoModes
    plot(RES(sid:end,2,i),RES(sid:end,3,i),"LineWidth",1.5,"LineStyle","-","Marker","o");
    hold on
    namelist = [namelist "Mode"+num2str(i)];
end
legend(namelist,"Location","east");
grid on
% ylim([-0.2 1]);
% xlim([00 9000]);
title("V-G Diagram")
xlabel("Velo (m/s)")
ylabel("Artificial Damping ")

% V-f diagram
figure(2)
hold on
namelist=[];
for i = 1:NoModes
    plot(RES(sid:end,2,i),RES(sid:end,5,i),"LineWidth",1.5,"LineStyle","-","Marker","o");
    namelist = [namelist "Mode"+num2str(i)];
end
legend(namelist,"Location","southeast");
grid on
% xlim([00 9000]);
% ylim([0 1000]);
title("V-F Diagram")
xlabel("Velo (m/s)")
ylabel("omega (rad/s)")

结果

通过V-G图人工阻尼由负向正穿越可以判断颤振的发生，系统需要增加正的人工阻尼才能维持简谐振动，即系统本身相对于平衡状态存在负阻尼，系统不稳定。从V-G图可以看出，系统二阶模态在横轴0.6（单位有误）位置发生穿越，该位置即颤振边界。而通过V-F图可以看出二阶模态于一阶模态有耦合的趋势，通过分离这两阶模态的频率有望提升颤振速度。通过对比MATLAB实现与ZAERO计算的数值结果发现，MATLAB实现的计算精度与成熟商软相当。