CN104965991A - 基于传递函数的机翼颤振速度确定方法 - Google Patents

Abstract

本发明涉及基于传递函数的机翼颤振速度确定方法，所述方法利用机翼弯扭振动微分方程和Theodrosen非定常气动力模型，得出长直机翼颤振微分方程，然后对所述机翼颤振微分方程进行Fourier变换，再运用传递函数法求解机翼颤振速度。本发明的有益效果如下：（1）由于本发明利用机翼弯扭振动微分方程来准确描述机翼振动，而没有采用机翼振动的低阶模态来近似描述机翼振动，所以计算的机翼颤振速度更加准确；（2）本发明求解机翼颤振速度的方法与现有技术相比更加简捷。

Description

基于传递函数的机翼颤振速度确定方法

技术领域

本发明涉及一种基于传递函数的机翼颤振速度确定方法，属于飞行器气动弹性技术领域。

背景技术

机翼颤振是指发生在飞机飞行中的动不稳定性，此时的飞行速度称为颤振速度。在飞行达到颤振速度时所发生的自激振动，大多数都会造成灾难性的后果。颤振分析就是求出颤振发生的条件，亦即求出颤振速度，并寻求在飞行器飞行速度范围内避免颤振发生的措施，进而寻找提高颤振临界速度的方法。

进行机翼颤振分析和计算时，一方面需要知道机翼结构的动力学特性，另一方面需要知道机翼结构附近非定常流动的空气动力特性。目前，工程中为了避免分析计算过于复杂，在机翼结构动力学计算时尽量减小机翼振动自由度，通常忽略高阶模态，仅利用机翼振动的低阶模态来近似表示机翼弯曲振动位移h和扭转振动转角α，通常选取前一阶、二阶弯曲模态和前一阶扭转模态来描述机翼的动力学特性，即

$\left\{\begin{array}{c}h=\underset{i=1}{\overset{\mathrm{\∞}}{\Σ}}{h}_i{\mathrm{\ξ}}_i\≈h_1{\mathrm{\ξ}}_1+h_2{\mathrm{\ξ}}_2\\ \mathrm{\α}=\underset{i=1}{\overset{\mathrm{\∞}}{\Σ}}{\mathrm{\α}}_i{\mathrm{\η}}_i\≈{\mathrm{\α}}_1{\mathrm{\η}}_1\end{array}\right.$

然后，再结合Theodrosen非定常气动理论，求解机翼颤振速度。

实际上，机翼是无限多自由度的弹性体，在颤振分析中也没有对自由度的限制。因此，上述计算方法的缺点是：为了避免过于复杂的计算而减小了机翼振动的自由度，对机翼实际振动作出了近似，使计算结果的准确性降低了。

发明内容

本发明所要解决的技术问题是提供一种计算方法简捷、计算结果较准确的基于传递函数的机翼颤振速度确定方法。

本发明解决其技术问题所需采用的技术方案：

一种基于传递函数的机翼颤振速度确定方法，所述方法利用机翼弯扭振动微分方程和Theodrosen非定常气动力模型，得出长直机翼颤振微分方程，然后对所述机翼颤振微分方程进行Fourier变换，再运用传递函数法求解机翼颤振速度；

一、所述方法需要依据的计算公式：

a.利用机翼弯扭振动微分方程建立机翼颤振微分方程：

本发明利用机翼弯扭振动微分方程式(1)来描述机翼弯扭振动：

$\left\{\begin{array}{c}EI\frac{{\∂}^{4}h}{\∂y^4}+m\frac{{\∂}^{2}h}{\∂t^2}+{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}-L_h=0\\ GJ\frac{{\∂}^2\mathrm{\α}}{\∂y^2}+I_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}+{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^{2}h}{\∂t^2}-T_{\mathrm{\α}}=0\end{array}\right.---\left(1\right)$

式中，h为机翼弯曲振动位移，单位为米，α为机翼扭转振动转角，单位为弧度，EI为机翼抗弯刚度，单位牛顿·米²，GJ为机翼抗扭刚度，单位为牛顿·米²，m为机翼单位长度质量，单位为千克，I_α为单位长度机翼绕弹性轴的转动惯量，单位为千克·米²，x_α为机翼弹性轴到机翼截面重心的距离，单位为米，L_k为机翼单位长度的升力，单位为牛顿/米，T_α机翼单位长度的扭矩，单位为牛顿，y为机翼展向坐标，单位为米，t为时间，单位为秒，为机翼弯曲振动位移h对机翼展向坐标y的四阶偏导数，为机翼扭转振动转角α对机翼展向坐标y的二阶偏导数，为机翼弯曲振动位移h对时间t的二阶偏导数，为机翼扭转振动转角α对时间t的二阶偏导数；

b.根据Theodrosen非定常气动力模型，得到机翼单位长度的升力L_k和机翼单位长度的扭矩T_α为下式(2)：

$\left\{\begin{array}{c}{L}_h={\mathrm{\π\ρb}}^2(-\frac{{\∂}^{2}h}{\∂t^2}+V\frac{\∂\mathrm{\α}}{\∂t}-b\stackrel{\‾}{a}\frac{{\∂}^2\mathrm{\α}}{\∂t^2})+2\mathrm{\π}\mathrm{\ρ}VbC\left(k\right)(V\mathrm{\α}-\frac{\∂h}{\∂t}+b(\frac{1}{2}-\stackrel{\‾}{a})\frac{\∂\mathrm{\α}}{\∂t})\\ T_{\mathrm{\α}}={\mathrm{\π\ρb}}^2(-b\stackrel{\‾}{a}\frac{{\∂}^{2}h}{\∂t^2}-Vb(\frac{1}{2}-\stackrel{\‾}{a})\frac{\∂\mathrm{\α}}{\∂t}-b^2(\frac{1}{8}+{\stackrel{\‾}{a}}^2)\frac{{\∂}^2\mathrm{\α}}{\∂t^2})\\ +2{\mathrm{\π\ρVb}}^2(\frac{1}{2}+\stackrel{\‾}{a})C\left(k\right)(V\mathrm{\α}-\frac{\∂h}{\∂t}+b(\frac{1}{2}-\stackrel{\‾}{a})\frac{\∂\mathrm{\α}}{\∂t})\end{array}\right.---\left(2\right)$

式中，V为空速，单位为米/秒，ρ为空气密度，单位为千克/立方米，b为机翼的半弦长，单位为米，C(k)为Theodrosen函数，为减缩频率，ω为圆频率，单位为弧度/秒，为机翼弹性轴到机翼弦长中点的距离占半弦长的百分比；

由于减缩频率k是圆频率ω和空速V的函数，将C(k)写为C(ω，V)，将(2)式代入(1)式得到机翼颤振微分方程式(3)：

$\left\{\begin{array}{c}EI\frac{{\∂}^{4}h}{\∂y^4}+m\frac{{\∂}^{2}h}{\∂t^2}-{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}-{\mathrm{\π\ρb}}^2\left(-\frac{{\∂}^{2}h}{\∂t^2}+V\frac{\∂\mathrm{\α}}{\∂t}-b\stackrel{\‾}{a}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}\right)\\ -2\mathrm{\π}\mathrm{\ρ}VbC\left(\mathrm{\ω},V\right)\left(V\mathrm{\α}-\frac{\∂h}{\∂t}+b\left(\frac{1}{2}-\stackrel{\‾}{a}\right)\frac{\∂\mathrm{\α}}{\∂t}\right)=0\\ GJ\frac{{\∂}^2\mathrm{\α}}{\∂y^2}-I_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}+{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^{2}h}{\∂t^2}-{\mathrm{\π\ρb}}^2\left(-b\stackrel{\‾}{a}\frac{{\∂}^{2}h}{\∂t^2}+Vb\left(\frac{1}{2}-\stackrel{\‾}{a}\right)\frac{\∂\mathrm{\α}}{\∂t}-b^2\left(\frac{1}{8}+{\stackrel{\‾}{a}}^2\right)\frac{{\∂}^2\mathrm{\α}}{\∂t^2}\right)\\ +2{\mathrm{\π\ρVb}}^2\left(\frac{1}{2}+\stackrel{\‾}{a}\right)C\left(\mathrm{\ω},V\right)\left(V\mathrm{\α}-\frac{\∂h}{\∂t}+b\left(\frac{1}{2}-\stackrel{\‾}{a}\right)\frac{\∂\mathrm{\α}}{\∂t}\right)=0\end{array}\right.---\left(3\right)$

c.利用传递函数法求解机翼颤振速度

①对(3)式作Fourier变换，并经过整理可得到下式(4)：

$\{\begin{array}{c}\frac{{\∂}^{4}h}{\∂y^4}=A_1\left(\mathrm{\ω},V\right)h+B_1\left(\mathrm{\ω},V\right)\mathrm{\α}\\ \frac{{\∂}^2\mathrm{\α}}{\∂y^2}=A_2\left(\mathrm{\ω},V\right)h+B_2\left(\mathrm{\ω},V\right)\mathrm{\α}\end{array}---\left(4\right)$

式中，A₁(ω，V)、A₂(ω，V)、B₁(ω，V)、B₂(ω，V)的具体表达式如下：

$\{\begin{array}{c}{A}_1\left(\mathrm{\ω},V\right)=\frac{{\mathrm{\ω}}^2\[\left(m-{\mathrm{\π\ρb}}^2\right)+2\mathrm{\π}\mathrm{\ρ}b\frac{V}{\mathrm{\ω}}C\left(\mathrm{\ω},V\right)i\]}{EI}\\ B_1\left(\mathrm{\ω},V\right)=\frac{{\mathrm{\ω}}^2\[{\mathrm{\π\ρb}}^2\left(b\stackrel{\‾}{a}+{\mathrm{mx}}_{\mathrm{\α}}\right)+\left({\mathrm{\π\ρb}}^2\frac{V}{\mathrm{\ω}}+2{\mathrm{\π\ρb}}^2\frac{V}{\mathrm{\ω}}C\left(\mathrm{\ω},V\right)\left(0.5-\stackrel{\‾}{a}\right)\right)i+2\mathrm{\π\ρb}\frac{V^2}{{\mathrm{\ω}}^2}C\left(\mathrm{\ω},V\right)\]}{EI}\\ A_2\left(\mathrm{\ω},V\right)=\frac{{\mathrm{\ω}}^2\[\left({\mathrm{mx}}_a-{\mathrm{\π\ρb}}^3\stackrel{\‾}{a}\right)+2\mathrm{\π}\mathrm{\ρ}\frac{V}{\mathrm{\ω}}{b}^2\left(0.5+\stackrel{\‾}{a}\right)C\left(\mathrm{\ω},V\right)i\]}{GJ}\\ B_2\left(\mathrm{\ω},V\right)=\frac{{\mathrm{\ω}}^2\[\left(I_{\mathrm{\α}}+{\mathrm{\π\ρb}}^4\left(0.125+{\stackrel{\‾}{a}}^2\right)\right)+\[2\mathrm{\π}\mathrm{\ρ}\frac{V}{\mathrm{\ω}}{b}^3\left(0.5+\stackrel{\‾}{a}\right)C\left(\mathrm{\ω},V\right)\left(0.5-\stackrel{\‾}{a}\right)-{\mathrm{\π\ρb}}^3\frac{V}{\mathrm{\ω}}\left(0.5-\stackrel{\‾}{a}\right)\]i+2{\mathrm{\π\ρb}}^2\frac{V^2}{{\mathrm{\ω}}^2}\left(0.5+\stackrel{\‾}{a}\right)C\left(\mathrm{\ω},V\right)\]}{GJ}\end{array}---\left(5\right)$

其中，虚数

②采用传递函数法求解上述(4)式，确定机翼的颤振速度；

为了便于应用传递函数理论，定义状态变量向量如下式(6)：

$\mathrm{\η}(y,\mathrm{\ω})={\left[\begin{array}{cccccc}h& \frac{\∂h}{\∂y}& \frac{{\∂}^{2}h}{\∂y^2}& \frac{{\∂}^{1}h}{\∂y^1}& \mathrm{\α}& \frac{\∂\mathrm{\α}}{\∂y}\end{array}\right]}^r---\left(6\right)$

式中，T表示向量转置；

从而，将(4)式写成状态方程的形式如下式(7)：

$\frac{\∂\mathrm{\η}(y,\mathrm{\ω})}{\∂y}=F\left(\mathrm{\ω},V\right)\mathrm{\η}\left(y,\mathrm{\ω}\right)+g\left(y,\mathrm{\ω}\right)---\left(7\right)$

式中，

$F\left(\mathrm{\ω},V\right)=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ A_1(\mathrm{\ω},V)& 0& 0& 0& B_1\left(\mathrm{\ω},V\right)& 0\\ 0& 0& 0& 0& 0& 1\\ A_2(\mathrm{\ω},V)& 0& 0& 0& B_2\left(\mathrm{\ω},V\right)& 0\end{array}\right]---\left(8\right)$

g(y，ω)＝0

边界条件为：

M_bη(y＝0，ω)+N_bη(y＝l，ω)＝γ(ω)   (9)

式中，η(0，ω)为η(y，ω)在机翼根端y＝0处的值，η(l，ω)为η(y，ω)在机翼梢端y＝l处的值，l为机翼的半展长，单位为米，M_b为机翼根端边界条件选择矩阵，N_b为机翼梢端边界条件选择矩阵，γ(ω)为由边界条件给定的位移或力组成的向量，其表达式分别为：

$M_b=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],\mathrm{\γ}\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}h(y=0,\mathrm{\ω})\\ \frac{\∂h}{\∂x}\left(y=0,\mathrm{\ω}\right)\\ \mathrm{\α}(y=0,\mathrm{\ω})\\ \frac{{\∂}^{2}h}{\∂x^2}(y=l,\mathrm{\ω})\\ \frac{{\∂}^{3}h}{\∂x^3}(y=l,\mathrm{\ω})\\ \frac{\∂\mathrm{\α}}{\∂x}(y=l,\mathrm{\ω})\end{array}\right\}$

根据传递函数理论，方程(7)式的解为：

$\mathrm{\η}(y,\mathrm{\ω})=H(y,\mathrm{\ω},V)\mathrm{\γ}\left(\mathrm{\ω}\right)+{\∫}_0^{l}G(y,\mathrm{\ξ},\mathrm{\ω},V)g(\mathrm{\ξ},\mathrm{\ω})d\mathrm{\ξ}---\left(10\right)$

式中，G(y，ξ，ω，V)为状态空间方程的域内传递函数，H(y，ω，V)为状态空间方程的边界传递函数，其表达式分别为：

$\left\{\begin{array}{c}G\left(y,\mathrm{\ξ},\mathrm{\ω},V\right)=\left\{\begin{array}{c}H\left(y,\mathrm{\ω},V\right)M_b{e}^{-F\left(\mathrm{\ω},V\right)\mathrm{\ξ}},\mathrm{\ξ}\≤y\\ -H\left(y,\mathrm{\ω},V\right)M_b{e}^{-F\left(\mathrm{\ω},V\right)\left(l-\mathrm{\ξ}\right)},\mathrm{\ξ}>y\end{array}\right.\\ H\left(y,\mathrm{\ω},V\right)=e^{F\left(\mathrm{\ω},V\right)y}{\[M_b+N_b{e}^{F\left(\mathrm{\ω},V\right)l}\]}^{-1}\end{array}\right.---\left(11\right)$

式中，变量ξ∈(0，l)，为机翼展向的坐标，单位为米；

机翼颤振时弯曲振动位移h和扭转振动转角α的振幅为非零常数，即(10)式有非零解，根据传递函数理论可知，(10)式有非零解的充分必要条件为：

det[M_b+N_be^F(ω，V)l]^-1＝0   (12)

令A＝[M_b+N_be^F(ω，V)L]^-1，由于A为复矩阵，其行列式值等于零的必要条件为矩阵行列式值的实部与虚部均为零，即

$\{\begin{array}{c}\mathrm{Re}\[\det A\]=0\\ \mathrm{Im}\[\det A\]=0\end{array}---\left(13\right)$

求解上述(13)式，可以得到满足方程组的空速V和圆频率ω，分别记为V_cz和ω_cz。其中，V_cz即为机翼的颤振速度，ω_cz即为机翼的颤振圆频率。

二、所述方法的具体步骤如下：

步骤(一)：测量长直机翼的下列物理参数：

半弦长b，单位为米；

半展长l，单位为米；

单位长度机翼质量m，单位为千克/米；

机翼弹性轴(2)到机翼弦长中点的距离占半弦长的百分比

机翼弹性轴(2)到机翼截面重心的距离x_α，单位为米；

单位长度机翼绕弹性轴的转动惯量I_α，单位为千克·米²；

机翼的抗弯刚度EI，单位为牛顿·米²；

机翼的抗扭刚度GJ，单位为牛顿·米²；

空气密度ρ，单位为千克/立方米；

步骤(二)：确定飞机机翼的空速V和圆频率ω的大致范围，假设为 $\{\begin{array}{c}V\∈\left(V_0,V_1\right)\\ \mathrm{\ω}\∈\left({\mathrm{\ω}}_0,{\mathrm{\ω}}_1\right)\end{array};$

步骤(三)：在 $\left\{\begin{array}{c}V\∈\left(V_0,V_1\right)\\ \mathrm{\ω}\∈\left({\mathrm{\ω}}_0,{\mathrm{\ω}}_1\right)\end{array}\right.$ 范围内划分合适步长ΔV和Δω，并进行离散，空速V和圆频率ω的取值为： $\{\begin{array}{cc}V=V_0+i\mathrm{\Δ}V& (i=0,1,2,3...)\\ \mathrm{\ω}={\mathrm{\ω}}_0+j\mathrm{\Δ}\mathrm{\ω}& (j=0,1,2,3...)\end{array};$

步骤(四)：取空速V＝V₀，圆频率ω依次取ω₀+jΔω(j＝0，1，2，3…)，将空速V和圆频率ω的取值与步骤(一)中的机翼各物理参数代入(5)式，得到系数A₁(ω，V)、A₂(ω，V)、B₁(ω，V)、_B2(ω，V)的值；

步骤(五)：将系数A₁(ω，V)、A₂(ω，V)、B₁(ω，V)、B₂(ω，V)的值代入(8)式，得到矩阵F(ω，V)的值；

步骤(六)：将矩阵F(ω，V)的值代入(13)式，计算Re[detA]和Im[detA]的值；

步骤(七)：再依次取空速V＝V₀+jΔV(j＝1，2，3…)，重复步骤(四)至步骤(六)，计算Re[detA]和Im[detA]的值；

步骤(八)：确定同时满足(13)式的空速V的值，即为机翼的颤振速度V_cz。

本发明的有益效果如下：

(1)由于本发明利用机翼弯扭振动微分方程来准确描述机翼振动，而没有采用机翼振动的低阶模态来近似描述机翼振动，所以计算的机翼颤振速度更加准确。

(2)本发明求解机翼颤振速度的方法与现有技术相比更加简捷。

附图说明

图1为长直机翼示意图；

图2为长直机翼弦向剖面示意图；

图3为实施例1的Re[detA]的等值线图(V∈(0，50)，ω∈(0，200π))；

图4为实施例1的Im[detA]的等值线图(V∈(0，50)，ω∈(0，200π))；

图5为图3和图4的等值线合成图(V∈(0，50)，ω∈(0，200π))；

图6为实施例1的Re[detA]的等值线图(V∈(30，40)，ω∈(30π，60π))；

图7为实施例1的Im[detA]的等值线图(V∈(30，40)，ω∈(30π，60π))；

图8为图6和图7的等值线合成图(V∈(30，40)，ω∈(30π，60π))。

在图1、2中，1——重心轴，2——弹性轴，3——弹性轴位置，4——重心轴位置。

具体实施方式

实施例1：

为了进一步说明本发明所述方法，本实施例1与参考文献(赵永辉的专著《气动弹性力学与控制》，北京：科学出版社，2006)的计算结果(Vcz＝36.6米/秒)相比较。

本实施例1的具体计算步骤如下：

步骤(一)：采用上述参考文献中所用的长直机翼的物理参数(见下表1)：

表1 长直机翼的物理参数

步骤(二)：确定飞机机翼的空速V和圆频率ω的大致范围为 $\left\{\begin{array}{c}V\∈(0,50)\\ \mathrm{\ω}\∈(0,200\mathrm{\π})\end{array}\right.,$ 划分步长 $\left\{\begin{array}{c}\mathrm{\Δ}V=2\\ \mathrm{\Δ}\mathrm{\ω}=10\mathrm{\π}\end{array}\right.,$ 依据发明内容部分中步骤(三)至步骤(七)，计算Re[detA]和Im[detA]的值。图3和图4给出了在 $\left\{\begin{array}{c}V\∈(0,50)\\ \mathrm{\ω}\∈(0,200\mathrm{\π})\end{array}\right.$ 范围内Re[detA]和Im[detA]的等值线图。图5给出了图3和图4的合成图。从图5中可以发现，在给定范围 $\left\{\begin{array}{c}V\∈(0,50)\\ \mathrm{\ω}\∈(0,200\mathrm{\π})\end{array}\right.$ 内虚线矩形区域约为 $\left\{\begin{array}{c}V\∈(30,40)\\ \mathrm{\ω}\∈(30\mathrm{\π},60\mathrm{\π})\end{array}\right.$ 内存在满足(13)式的空速V和圆频率ω，即机翼颤振速度V_cz在上述范围内。

为了提高计算结果精度，在上述缩小的范围内 $\left\{\begin{array}{c}V\∈(30,40)\\ \mathrm{\ω}\∈(30\mathrm{\π},60\mathrm{\π})\end{array}\right.$ 内，划分更小的步长 $\left\{\begin{array}{c}\mathrm{\Δ}V=0.2\\ \mathrm{\Δ}\mathrm{\ω}=0.4\mathrm{\π}\end{array}\right.,$ 并依据所述步骤(三)到步骤(七)，再计算Re[detA]和Im[detA]的值，获得精度更高的机翼颤振速度值。图6和图7分别给出了在缩小范围 $\left\{\begin{array}{c}V\∈(30,40)\\ \mathrm{\ω}\∈(30\mathrm{\π},60\mathrm{\π})\end{array}\right.$ 内Re[detA]和Im[detA]的等值线图，图8给出了图6和图7的合成图。结合图6至图8，根据(13)式的要求，可从图8中找到画圆的区域内的O点为满足(13)式的空速V和圆频率ω。由图8可得：空速V＝35.3米/秒，此时的空速V＝35.3米/秒即机翼颤振速度，与上述参考文献的计算结果36.6米/秒吻合很好。

Claims (1)

1.一种基于传递函数的机翼颤振速度确定方法，其特征在于所述方法利用机翼弯扭振动微分方程和Theodrosen非定常气动力模型，得出长直机翼颤振微分方程，然后对所述机翼颤振微分方程进行Fourier变换，再运用传递函数法求解机翼颤振速度；

一、所述方法需要依据的计算公式：

a.利用机翼弯扭振动微分方程建立机翼颤振微分方程：

本发明利用机翼弯扭振动微分方程式(1)来描述机翼弯扭振动：

$\left\{\begin{array}{c}EI\frac{{\∂}^{4}h}{\∂y^4}+m\frac{{\∂}^{2}h}{\∂t^2}-{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}-L_h=0\\ GJ\frac{{\∂}^2\mathrm{\α}}{\∂y^2}-I_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}+{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^{2}h}{\∂t^2}+T_{\mathrm{\α}}=0\end{array}\right.---\left(1\right)$

式中，h为机翼弯曲振动位移，单位为米，α为机翼扭转振动转角，单位为弧度，EI为机翼抗弯刚度，单位牛顿·米²，GJ为机翼抗扭刚度，单位为牛顿·米²，m为机翼单位长度质量，单位为千克，I_α为单位长度机翼绕弹性轴的转动惯量，单位为千克·米²，x_α为机翼弹性轴到机翼截面重心的距离，单位为米，L_h为机翼单位长度的升力，单位为牛顿/米，T_α机翼单位长度的扭矩，单位为牛顿，y为机翼展向坐标，单位为米，t为时间，单位为秒，为机翼弯曲振动位移h对机翼展向坐标y的四阶偏导数，为机翼扭转振动转角α对机翼展向坐标y的二阶偏导数，为机翼弯曲振动位移h对时间t的二阶偏导数，为机翼扭转振动转角α对时间t的二阶偏导数；

b.根据Theodrosen非定常气动力模型，得到机翼单位长度的升力L_h和机翼单位长度的扭矩T_α为下式(2)：

$\left\{\begin{array}{c}{L}_h={\mathrm{\π\ρb}}^2(-\frac{{\∂}^{2}h}{\∂t^2}+V\frac{\∂\mathrm{\α}}{\∂t}-b\stackrel{\‾}{a}\frac{{\∂}^2\mathrm{\α}}{\∂t^2})+2\mathrm{\π}\mathrm{\ρ}VbC\left(k\right)(V\mathrm{\α}-\frac{\∂h}{\∂t}+b(\frac{1}{2}-\stackrel{\‾}{a})\frac{\∂\mathrm{\α}}{\∂t})\\ T_{\mathrm{\α}}={\mathrm{\π\ρb}}^2(-b\stackrel{\‾}{a}\frac{{\∂}^{2}h}{\∂t^2}-Vb(\frac{1}{2}-\stackrel{\‾}{a})\frac{\∂\mathrm{\α}}{\∂t}-b^2(\frac{1}{8}+{\stackrel{\‾}{a}}^2)\frac{{\∂}^2\mathrm{\α}}{\∂t^2})\\ +2{\mathrm{\π\ρVb}}^2\left(\frac{1}{2}+\stackrel{\‾}{a}\right)C\left(k\right)(V\mathrm{\α}-\frac{\∂h}{\∂t}+b(\frac{1}{2}-\stackrel{\‾}{a})\frac{\∂\mathrm{\α}}{\∂t})\end{array}\right.---\left(2\right)$

式中，V为空速，单位为米/秒，ρ为空气密度，单位为千克/立方米，b为机翼的半弦长，单位为米，C(k)为Theodrosen函数，为减缩频率，ω为圆频率，单位为弧度/秒，为机翼弹性轴到机翼弦长中点的距离占半弦长的百分比；由于减缩频率k是圆频率ω和空速V的函数，将C(k)写为C(ω,V)，将(2)式代入(1)式得到机翼颤振微分方程式(3)：

$\left\{\begin{array}{c}EI\frac{{\∂}^{4}h}{\∂y^4}+m\frac{{\∂}^{2}h}{\∂t^2}-{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}-{\mathrm{\π\ρb}}^2\left(-\frac{{\∂}^{2}h}{\∂t^2}+V\frac{\∂\mathrm{\α}}{\∂t}-b\stackrel{\‾}{a}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}\right)\\ -2\mathrm{\π}\mathrm{\ρ}VbC(\mathrm{\ω},V)\left(V\mathrm{\α}-\frac{\∂h}{\∂t}+b\left(\frac{1}{2}-\stackrel{\‾}{a}\right)\frac{\∂\mathrm{\α}}{\∂t}\right)=0\\ GJ\frac{{\∂}^{2}h}{\∂y^2}-I_{\mathrm{\α}}\frac{{\∂}^2\mathrm{\α}}{\∂t^2}+{\mathrm{mx}}_{\mathrm{\α}}\frac{{\∂}^{2}h}{\∂t^2}+{\mathrm{\π\ρb}}^2\left(-b\stackrel{\‾}{a}\frac{{\∂}^{2}h}{\∂t^2}-Vb\left(\frac{1}{2}-\stackrel{\‾}{a}\right)\frac{\∂\mathrm{\α}}{\∂t}-b^2\left(\frac{1}{8}+{\stackrel{\‾}{a}}^2\right)\frac{{\∂}^2\mathrm{\α}}{\∂t^2}\right)\\ +2{\mathrm{\π\ρVb}}^2\left(\frac{1}{2}+\stackrel{\‾}{a}\right)C(\mathrm{\ω},V)\left(V\mathrm{\α}-\frac{\∂h}{\∂t}+b\left(\frac{1}{2}-\stackrel{\‾}{a}\right)\frac{\∂\mathrm{\α}}{\∂t}\right)=0\end{array}\right.---\left(3\right)$

c.利用传递函数法求解机翼颤振速度

①对(3)式作Fourier变换，并经过整理可得到下式(4)：

$\left\{\begin{array}{c}\frac{{\∂}^4\mathrm{\α}}{\∂y^4}=A_1(\mathrm{\ω},V)h+B_1(\mathrm{\ω},V)\mathrm{\α}\\ \frac{{\∂}^2\mathrm{\α}}{\∂y^2}=A_2(\mathrm{\ω},V)h+B_2(\mathrm{\ω},V)\mathrm{\α}\end{array}\right.---\left(4\right)$

式中，A₁(ω,V)、A₂(ω,V)、B₁(ω,V)、B₂(ω,V)的具体表达式如下：

$\{\begin{array}{c}{A}_1(\mathrm{\ω},V)=\frac{{\mathrm{\ω}}^2\[(m-{\mathrm{\π\ρb}}^2)+2\mathrm{\π}\mathrm{\ρ}b\frac{V}{\mathrm{\ω}}C(\mathrm{\ω},V)i\]}{EI}\\ B_1(\mathrm{\ω},V)=\frac{{\mathrm{\ω}}^2\[{\mathrm{\π\ρb}}^2(b\stackrel{\‾}{a}+{\mathrm{mx}}_{\mathrm{\α}})+({\mathrm{\π\ρb}}^2\frac{V}{\mathrm{\ω}}+2{\mathrm{\π\ρb}}^2\frac{V}{\mathrm{\ω}}C(\mathrm{\ω},V\left)\right(0.5-\stackrel{\‾}{a}\left)\right)i+2\mathrm{\π}\mathrm{\ρ}b\frac{V^2}{{\mathrm{\ω}}^2}C(\mathrm{\ω},V)\]}{EI}\\ A_2(\mathrm{\ω},V)=\frac{{\mathrm{\ω}}^2\[({\mathrm{mx}}_{\mathrm{\α}}-{\mathrm{\π\ρb}}^3\stackrel{\‾}{a})+2\mathrm{\π}\mathrm{\ρ}\frac{V}{\mathrm{\ω}}{b}^2\left(0.5+\stackrel{\‾}{a}\right)C\left(\mathrm{\ω},V\right)i\]}{GJ}\\ B_2(\mathrm{\ω},V)=\frac{{\mathrm{\ω}}^2\[(I_{\mathrm{\α}}+{\mathrm{\π\ρb}}^4(0.125+{\stackrel{\‾}{a}}^2\left)\right)+\[2\mathrm{\π}\mathrm{\ρ}\frac{V}{\mathrm{\ω}}{b}^3\left(0.5+\stackrel{\‾}{a}\right)C\left(\mathrm{\ω},V\right)\left(0.5-\stackrel{\‾}{a}\right)-{\mathrm{\π\ρb}}^3\frac{V}{\mathrm{\ω}}\left(0.5-\stackrel{\‾}{a}\right)\]i+2{\mathrm{\π\ρb}}^2\frac{V^2}{{\mathrm{\ω}}^2}\left(0.5+\stackrel{\‾}{a}\right)C(\mathrm{\ω},V)\]}{GJ}\end{array}---\left(5\right)$

其中，虚数 $i=\sqrt{-1};$

②采用传递函数法求解上述(4)式，确定机翼的颤振速度；

为了便于应用传递函数理论，定义状态变量向量如下式(6)：

$\mathrm{\η}(y,\mathrm{\ω})={\left[\begin{array}{cccccc}h& \frac{\∂h}{\∂y}& \frac{{\∂}^{2}h}{\∂y^2}& \frac{{\∂}^{3}h}{\∂y^3}& \mathrm{\α}& \frac{\∂\mathrm{\α}}{\∂y}\end{array}\right]}^T---\left(6\right)$

式中，T表示向量转置；

从而，将(4)式写成状态方程的形式如下式(7)：

$\frac{\∂\mathrm{\η}(y,\mathrm{\ω})}{\∂y}=F(\mathrm{\ω},V)\mathrm{\η}(y,\mathrm{\ω})+g(y,\mathrm{\ω})---\left(7\right)$

式中，

$F(\mathrm{\ω},V)=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ A_1(\mathrm{\ω},V)& 0& 0& 0& B_1(\mathrm{\ω},V)& 0\\ 0& 0& 0& 0& 0& 1\\ A_2(\mathrm{\ω},V)& 0& 0& 0& B_2(\mathrm{\ω},V)& 0\end{array}\right]---\left(8\right)$

g(y,ω)＝0

边界条件为：

M_bη(y＝0,ω)+N_bη(y＝l,ω)＝γ(ω)        (9)

式中，η(0,ω)为η(y,ω)在机翼根端y＝0处的值，η(l,ω)为η(y,ω)在机翼梢端y＝l处的值，l为机翼的半展长，单位为米，M_b为机翼根端边界条件选择矩阵，N_b为机翼梢端边界条件选择矩阵，γ(ω)为由边界条件给定的位移或力组成的向量，其表达式分别为：

$M_b=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],\mathrm{\γ}\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}h(y=0,\mathrm{\ω})\\ \frac{\∂h}{\∂x}(y=0,\mathrm{\ω})\\ \mathrm{\α}(y=0,\mathrm{\ω})\\ \frac{{\∂}^{2}h}{\∂x^2}(y=l,\mathrm{\ω})\\ \frac{{\∂}^{3}h}{\∂x^3}(y=l,\mathrm{\ω})\\ \frac{\∂\mathrm{\α}}{\∂x}(y=l,\mathrm{\ω})\end{array}\right\}$

根据传递函数理论，方程(7)式的解为：

$\mathrm{\η}(y,\mathrm{\ω})=H(y,\mathrm{\ω},V)\mathrm{\γ}\left(\mathrm{\ω}\right)+{\∫}_0^{l}G(y,\mathrm{\ξ},\mathrm{\ω},V)g(\mathrm{\ξ},\mathrm{\ω})d\mathrm{\ξ}---\left(10\right)$

式中，G(y,ξ,ω,V)为状态空间方程的域内传递函数，H(y,ω,V)为状态空间方程的边界传递函数，其表达式分别为：

$\left\{\begin{array}{c}G(y,\mathrm{\ξ},\mathrm{\ω},V)=\left\{\begin{array}{c}H(y,\mathrm{\ω},V)M_b{e}^{-F(\mathrm{\ω},V)\mathrm{\ξ}},\mathrm{\ξ}\≤y\\ -H(y,\mathrm{\ω},V)N_b{e}^{-F(\mathrm{\ω},V)(l-\mathrm{\ξ})},\mathrm{\ξ}>y\end{array}\right.\\ H(y,\mathrm{\ω},V)=e^{F(\mathrm{\ω},V)y}{\[M_b+N_b{e}^{F(\mathrm{\ω},V)l}\]}^{-1}\end{array}\right.---\left(11\right)$

式中，变量ξ∈(0,l)，为机翼展向的坐标，单位为米；

机翼颤振时弯曲振动位移h和扭转振动转角α的振幅为非零常数，即(10)式有非零解，根据传递函数理论可知，(10)式有非零解的充分必要条件为：

det[M_b+N_be^F(ω,V)l]^-1＝0           (12)

令A＝[M_b+N_be^F(ω,V)L]^-1，由于A为复矩阵，其行列式值等于零的必要条件为矩阵行列式值的实部与虚部均为零，即

$\left\{\begin{array}{c}\mathrm{Re}\[\det A\]=0\\ \mathrm{Im}\[\det A\]=0\end{array}\right.---\left(13\right)$

求解上述(13)式，可以得到满足方程组的空速V和圆频率ω，分别记为V_cz和ω_cz；

其中，V_cz即为机翼的颤振速度，ω_cz即为机翼的颤振圆频率；

二、所述方法的具体步骤如下：

步骤(一)：测量长直机翼的下列物理参数：

半弦长b，单位为米；

半展长l，单位为米；

单位长度机翼质量m，单位为千克/米；

机翼弹性轴(2)到机翼弦长中点的距离占半弦长的百分比

机翼弹性轴(2)到机翼截面重心的距离x_α，单位为米；

单位长度机翼绕弹性轴的转动惯量I_α，单位为千克·米²；

机翼的抗弯刚度EI，单位为牛顿·米²；

机翼的抗扭刚度GJ，单位为牛顿·米²；

空气密度ρ，单位为千克/立方米；

步骤(二)：确定飞机机翼的空速V和圆频率ω的大致范围，假设为 $\left\{\begin{array}{c}V\∈(V_0,V_1)\\ \mathrm{\ω}\∈({\mathrm{\ω}}_0,{\mathrm{\ω}}_1)\end{array}\right.;$

步骤(三)：在 $\left\{\begin{array}{c}V\∈(V_0,V_1)\\ \mathrm{\ω}\∈({\mathrm{\ω}}_0,{\mathrm{\ω}}_1)\end{array}\right.$ 范围内划分合适步长ΔV和Δω，并进行离散，空速V和圆频率ω的取值为： $\left\{\begin{array}{cc}V=V_0+i\mathrm{\Δ}V& (i=0,1,2,3...)\\ \mathrm{\ω}={\mathrm{\ω}}_0+j\mathrm{\Δ}\mathrm{\ω}& (j=0,1,2,3...)\end{array}\right.;$

步骤(四)：取空速V＝V₀，圆频率ω依次取ω₀+jΔω(j＝0,1,2,3…)，将空速V和圆频率ω的取值与步骤(一)中的机翼各物理参数代入(5)式，得到系数A₁(ω,V)、A₂(ω,V)、B₁(ω,V)、B₂(ω,V)的值；

步骤(五)：将系数A₁(ω,V)、A₂(ω,V)、B₁(ω,V)、B₂(ω,V)的值代入(8)式，得到矩阵F(ω,V)的值；

步骤(六)：将矩阵F(ω,V)的值代入(13)式，计算Re[det A]和Im[det A]的值；

步骤(七)：再依次取空速V＝V₀+jΔV(j＝1,2,3…)，重复步骤(四)至步骤(六)，计算Re[det A]和Im[det A]的值；

步骤(八)：确定同时满足(13)式的空速V的值，即为机翼的颤振速度V_cz。