气动弹性颤振涉猎

• 1 绪论
• 2 特征值分析的稳定性特征（p法）
• 3 典型翼段的气动弹性分析（p法-定常气动理论）
• 4 经典颤振分析（k法-无阻尼）
  • 4.1 单自由度颤振
  • 4.2 二自由度颤振
• 5 颤振问题的工程求解（k法、p-k法）
  • 5.1 k法
  • 5.2 p-k法
    • 5.2.1 p法和k法的比较
    • 5.2.2 p-k法简介
• 6 非定常气动力
  • 6.1 Theodorsen非定常薄翼理论
  • 6.2 Peter等人的有限状态非定常薄翼理论
• 7 根据假设模态进行颤振预测
• 8 颤振边界特性
• 9 总结

1 绪论

气动弹性颤振是飞行器结构在气动力、弹性力和惯性力三者共同作用下的一种动力学不稳定。颤振是自激振动且具有潜在破坏性的振动不稳定性，它是由在弹性物体上的气动力与它的固有振动模态耦合产生的，有增长振幅的振动运动。

涉及知识：空气动力学、结构动力学、理论力学、材料力学、高等数学、线性代数

2 特征值分析的稳定性特征（p法）

将升力面简化为板模型。令$w(x,y,t)$为在$z$向垂直于机翼平面（即$x-y$平面）的位移，${\varphi }_i(x,y)$为固有振动模态，$w_i$为对应的各阶固有频率，（${\varphi }_i(x,y)$和$w_i$只与模型的结构和边界条件有关，故视为已知量）则典型的结构位移可表示为

$$w(x,y,t)=\sum _{i=1}^n{\xi }_i(t){\varphi }_i(x,y)$$

其中${\xi }_i(t)$为第$i$阶模态的广义坐标。广义运动方程组：

$$M_i({\ddot{\xi }}_i+{\omega }_i^2{\xi }_i)={\Xi }_i(i=0,1,2,\cdots ,n)$$

其中$M_i$为广义质量，与质量分布$m(x,y)$相关，并确定为

$$M_i={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}m(x,y){\varphi }_i^2(x,y)\mathrm{d}x\mathrm{d}y$$

${\Xi }_i(t)$为广义力，与外力载荷$F(x,y,t)$相关，并计算为

$${\Xi }_i={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}F(x,y,t){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y$$

外载荷只考虑气动载荷。压强分布$\Delta p(x,y,t)$可以用广义坐标及其导数的线性方程来表示：

$$\begin{gathered} \Delta p(x,y,t)=\sum _{j=1}^n[a_j(x,y){\xi }_j(t)+b_j(x,y){\dot{\xi }}_j(t)+c_j(x,y){\ddot{\xi }}_j(t)] \\ +\sum _{j=1}^N{d}_j(x,y){\lambda }_j(t) \end{gathered}$$

其中${\lambda }_j$表示流场相关的状态变量。于是广义力表示为

$$\begin{array}{rl}{\Xi }_i(t)& ={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}\Delta p(x,y,t){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y\\ & =\sum _{j=1}^n{\xi }_j(t){\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{a}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y+\\ & \sum _{j=1}^n{\dot{\xi }}_j(t){\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{b}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y+\\ & \sum _{j=1}^n{\ddot{\xi }}_j(t){\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{c}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y+\\ & \sum _{j=1}^N{\lambda }_j(t){\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{d}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y\\ & ={\rho }_{\infty }\frac{U^2}{b^2}[\sum _{j=1}^n(a_{ij}{\xi }_j+\frac{b}{U}{b}_{ij}{\dot{\xi }}_j+\frac{b^2}{U^2}{c}_{ij}{\ddot{\xi }}_j)+\sum _{j=1}^N{d}_{ij}{\lambda }_j]\end{array}$$

其中

$$\begin{gathered} {\rho }_{\infty }\frac{U^2}{b^2}{a}_{ij}={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{a}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y \\ {\rho }_{\infty }\frac{U}{b}{b}_{ij}={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{b}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y \\ {\rho }_{\infty }{c}_{ij}={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{c}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y \\ {\rho }_{\infty }\frac{U^2}{b^2}{d}_{ij}={\iint }_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}}{d}_j(x,y){\varphi }_i(x,y)\mathrm{d}x\mathrm{d}y \end{gathered}$$

将相关量代入广义运动方程：

$$\begin{gathered} \frac{b^2}{U^2}(M_i{\ddot{\xi }}_i+M_i{\omega }_i^2\xi )-{\rho }_{\infty }\frac{b^2}{U^2}\sum _{j=1}^n{c}_{ij}{\ddot{\xi }}_j-{\rho }_{\infty }\frac{b}{U}\sum _{j=1}^n{b}_{ij}{\dot{\xi }}_j- \\ {\rho }_{\infty }\sum _{j=1}^n{a}_{ij}{\xi }_j-{\rho }_{\infty }\sum _{j=1}^n{d}_{ij}{\lambda }_j=0(i=1,2,\cdots ,n)\text{\hspace{1em}(2.1)} \end{gathered}$$

同时由空气动力学理论可以得到$N$个附加方程：

$$\sum _{j=1}^N{A}_{ij}{\dot{\lambda }}_j+\frac{U}{b}({\lambda }_i-\sum _{j=1}^n{E}_{ij}{\xi }_j)=0(i=1,2,\cdots ,N)\text{\hspace{1em}(2.2)}$$

由$(2.1)$和$(2.2)$联立，构成包含了$n+N$个方程的齐次线性常微分方程组，其中有$n+N$个未知函数。于是解的形式可以取作：

$${\xi }_i(t)={\overline{\xi }}_i\mathrm{e}\mathrm{x}\mathrm{p}(vt),{\lambda }_i(t)={\overline{\lambda }}_i\mathrm{e}\mathrm{x}\mathrm{p}(vt)$$

代入方程$(2.1)$和$(2.2)$中，得到关于${\overline{\xi }}_i$和${\overline{\lambda }}_i$的$n+N$阶线性齐次代数方程组，矩阵形式为

$$\begin{array}{c}[p^2\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})+\frac{b^2}{U^2}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M}{\omega }^2)]\overline{\mathbf{\xi }}-{\rho }_{\infty }(p^2\mathbf{c}+p\mathbf{b}+\\ \mathbf{a})\overline{\mathbf{\xi }}-{\rho }_{\infty }\mathbf{d}\overline{\mathbf{\lambda }}=0\\ -\mathbf{E}\overline{\mathbf{\xi }}+(p\mathbf{A}+\mathbf{I})\overline{\mathbf{\lambda }}=0\end{array}\}\text{\hspace{1em}(2.3)}$$

其中$p=bv/U$为未知无量纲特征值，$\mathbf{I}$为单位矩阵。

要求解广义坐标振幅的非平凡解，系数${\overline{\xi }}_i$和${\overline{\lambda }}_i$构成的行列式必须为0。即：

$$\left|\begin{array}{cc}{p}^2\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})+\frac{b^2}{U^2}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M}{\omega }^2)-{\rho }_{\infty }(p^2\mathbf{c}+p\mathbf{b}+\mathbf{a})& {\rho }_{\infty }\mathbf{d}\\ -\mathbf{E}& p\mathbf{A}+\mathbf{I}\end{array}\right|=0\text{\hspace{1em}(2.4)}$$

显然这个行列式是关于$p$的$2n+N$阶多项式。于是得到$2n+N$个根，包括$n_c$对共轭复根和$n_r$个实数根$(2n_c+n_r=2n+N)$。根的典型形式写作

$$v_k=\frac{U{p}_k}{b}={\Gamma }_k\pm \mathrm{i}{\Omega }_k(k=0,1,\cdots ,n_c,\cdots ,n_c+n_r)\text{\hspace{1em}(2.5)}$$

其中，当$n_c\le k\le n_c+n_r$时，${\Omega }_k=0$。
将每一个$p_k$代入$(2.3)$可以解得对应n维列阵${\overline{\mathbf{\xi }}}^{(k)}(i=0,1,\cdots ,n)$和N维列阵${\overline{\mathbf{\lambda }}}^{(\mathbf{k})}(i=0,1,\cdots ,N)$

于是位移场的解为

$$\begin{gathered} w(x,y,t)=\sum _{k=1}^{n_c+n_r}\{w_k(x,y)\mathrm{e}\mathrm{x}\mathrm{p}[({\Gamma }_k+\mathrm{i}{\Omega }_k)t]+ \\ {\overline{w}}_k(x,y)\mathrm{e}\mathrm{x}\mathrm{p}[({\Gamma }_k+\mathrm{i}{\Omega }_k)t]\} \\ {w}_k(x,y)=\sum _{k=1}^n{\overline{\xi }}_i^{(k)}{\varphi }_i(x,y)(k=1,2,\cdots ,n_c+n_r)\text{\hspace{1em}(2.6)} \end{gathered}$$

${\overline{w}}_k$为$w_k$的共轭复数。

当$\Omega \ne 0$时各典型模态幅值的状态如图

则可以得出结论${\Gamma }_k$和${\Omega }_k$不同组合对应的运动状态类型和稳定特性：

${\Gamma }_k$	${\Omega }_k$	运动类型	稳定特征
$<0$	$\ne 0$	衰减运动	稳定
$=0$	$\ne 0$	简谐运动	临界稳定
$>0$	$\ne 0$	发散运动	不稳定
$<0$	$=0$	收敛	稳定
$=0$	$=0$	与时间相关	临界稳定
$>0$	$=0$	发散	不稳定

可知，要求临界颤振条件只需要令${\Gamma }_k=0$即可。

$-{\Gamma }_k$称为第$k$阶模态的模态阻尼，${\Omega }_k$称为$k$阶模态的模态频率。

3 典型翼段的气动弹性分析（p法-定常气动理论）

建立如图简单的、弹簧支承的刚硬机翼模型

该模型所关注的是P、C、Q和T，它们分别是参考基准点（即沉浮位移$h$的测量点），机翼质心、焦心（亚声速薄翼理论中假定为1/4弦长点）和3/4弦长点（薄翼理论中的一个重要弦向位置）。无量纲参数$e$、$a$（$-1\le e\le 1$、$-1\le a\le 1$）决定了C、P的位置。$b$为翼型半弦长。记$x_{\theta }=e-a$为静不稳定参数。模型的刚体沉浮、俯仰运动由线性轻质弹簧约束，弹簧常数分别为$k_h$和$k_{\theta }$。

系统势能（不考虑重力势能）为

$$P=\frac{1}{2}{k}_h{h}^2+\frac{1}{2}{k}_{\theta }{\theta }^2\text{\hspace{1em}(3.1)}$$

质心C的速度为

$${\mathbf{v}}_C={\mathbf{v}}_P+\dot{\theta }{\mathbf{b}}_3\times b[(1+a)-(1+e)]{\mathbf{b}}_1$$

其中，参考点P的惯性速度为

$${\mathbf{v}}_P=-\dot{h}{\mathbf{i}}_2$$

则：

$${\mathbf{v}}_C=-\dot{h}{\mathbf{i}}_2-b{x}_{\theta }\dot{\theta }{\mathbf{b}}_2$$

动能为

$$K=\frac{1}{2}m{\mathbf{v}}_C\cdot {\mathbf{v}}_C+\frac{1}{2}{I}_C{\dot{\theta }}^2$$

假设$\theta$为小量，得

$$\begin{array}{rl}K& =\frac{1}{2}m({\dot{h}}^2+b^2{x}_{\theta }^2{\dot{\theta }}^2+2b{x}_{\theta }\dot{h}\dot{\theta })+\frac{1}{2}{I}_C{\dot{\theta }}^2\\ & =\frac{1}{2}m({\dot{h}}^2+2b{x}_{\theta }\dot{h}\dot{\theta })+\frac{1}{2}{I}_P{\dot{\theta }}^2\end{array}\text{\hspace{1em}(3.2)}$$

Q点的虚位移为

$$\delta {\mathbf{p}}_Q=-\delta h{\mathbf{i}}_2+b\delta \theta (\frac{1}{2}+a){\mathbf{b}}_2$$

机翼的虚拟转角为$\delta \theta {\mathbf{b}}_3$
于是虚功为

$$\overline{\delta W}=L[-\delta h+b(\frac{1}{2}+a)\delta \theta ]+M_{\frac{1}{4}}\delta \theta$$

因此广义力为

$$\begin{array}{c}{Q}_h=-L\\ Q_{\theta }=M_{\frac{1}{4}}+b(\frac{1}{2}+a)L\end{array}\}\text{\hspace{1em}(3.3)}$$

将$(3.1)$、$(3.2)$、$(3.3)$代入拉格朗日方程中得

$$\begin{array}{c}m(\ddot{h}+b{x}_{\theta }\ddot{\theta })+k_{h}h=-L\\ I_P\ddot{\theta }+mb{x}_{\theta }\ddot{h}+k_{\theta }\theta =M_{\frac{1}{4}}+b(\frac{1}{2}+a)L\end{array}\}\text{\hspace{1em}(3.4)}$$

这里采用定常气动理论和薄翼理论：

$$\begin{array}{c}L=2\pi {\rho }_{\infty }b{U}^2\theta \\ M_{\frac{1}{4}}=0\end{array}\}$$

引入零空速时的非耦合结构固有频率：

$${\omega }_h=\sqrt{\frac{k_h}{m}},{\omega }_{\theta }=\sqrt{\frac{k_{\theta }}{m}}$$

将相关量代入$(3.4)$，写作矩阵形式：

$$\left[\begin{array}{cc}m{b}^2& m{b}^2{x}_{\theta }\\ m{b}^2{x}_{\theta }& I_P\end{array}\right]\left[\begin{array}{c}\frac{\ddot{h}}{b}\\ \ddot{\theta }\end{array}\right]+\left[\begin{array}{cc}m{b}^2{\omega }_h^2& 2\pi {\rho }_{\infty }{b}^2{U}^2\\ 0& I_P{\omega }_{\theta }^2-2(\frac{1}{2}+a)\pi {\rho }_{\infty }{b}^2{U}^2\end{array}\right]\left[\begin{array}{c}\frac{h}{b}\\ \theta \end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

设此齐次线性常微分方程组的解结构为$h=\overline{h}\mathrm{e}\mathrm{x}\mathrm{p}(vt)$和$\theta =\overline{\theta }\mathrm{e}\mathrm{x}\mathrm{p}(vt)$，并设$v=\Gamma \pm \mathrm{i}\mathrm{\Omega }$。得到

$$\left[\begin{array}{cc}m{b}^2{v}^2+m{b}^2{\omega }_h^2& m{b}^2{v}^2{x}_{\theta }+2\pi {\rho }_{\infty }{b}^2{U}^2\\ m{b}^2{v}^2{x}_{\theta }& I_P{\omega }_{\theta }^2\theta +I_P{v}^2-2(a+\frac{1}{2})\pi {\rho }_{\infty }{b}^2{U}^2\end{array}\right]\left[\begin{array}{c}\frac{\overline{h}}{b}\\ \overline{\theta }\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

引入无量纲参数：

$$p=\frac{vb}{U},r^2=\frac{I_P}{m{b}^2},\sigma =\frac{{\omega }_h}{{\omega }_{\theta }},\mu =\frac{m}{{\rho }_{\infty }\pi b^2},V=\frac{U}{b{\omega }_{\theta }}$$

其中$p$和$V$是未知的无量纲复特征值。$V$表示无量纲空气的来流速度，也称为“减缩速度”。于是：

$$\left[\begin{array}{cc}{p}^2+\frac{{\sigma }^2}{V^2}& x_{\theta }{p}^2+\frac{2}{\mu }\\ x_{\theta }{p}^2& r^2{p}^2+\frac{r^2}{V^2}-\frac{2}{\mu }(a+\frac{1}{2})\end{array}\right]\left[\begin{array}{c}\frac{\overline{h}}{b}\\ \overline{\theta }\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]\text{\hspace{1em}(3.5)}$$

若存在非平凡解，则系数矩阵的行列式必须等于零：

$$\left|\begin{array}{cc}{p}^2+\frac{{\sigma }^2}{V^2}& x_{\theta }{p}^2+\frac{2}{\mu }\\ x_{\theta }{p}^2& r^2{p}^2+\frac{r^2}{V^2}-\frac{2}{\mu }(a+\frac{1}{2})\end{array}\right|=0\text{\hspace{1em}(3.6)}$$

从而得到关于$p$和$V$的方程。

将$V$视为函数自变量，$p$视为函数因变量。上面方程可以得到两组$p$的共轭复根：

$$\begin{array}{c}{p}_1=\frac{b{v}_1}{U}=\frac{b}{U}({\Gamma }_1\pm \mathrm{i}{\mathrm{\Omega }}_1)=p_1(V)\\ p_2=\frac{b{v}_2}{U}=\frac{b}{U}({\Gamma }_2\pm \mathrm{i}{\mathrm{\Omega }}_2)=p_2(V)\end{array}\}$$

它们分别对应俯仰($\theta$）频率和沉浮（$h$）频率。
乘以减缩速度$V$得

$$\begin{array}{c}V{p}_1=\frac{b}{U}({\Gamma }_1\pm \mathrm{i}{\mathrm{\Omega }}_1)\frac{U}{b{\omega }_{\theta }}=\frac{{\Gamma }_1}{{\omega }_{\theta }}\pm \mathrm{i}\frac{{\Omega }_1}{{\omega }_{\theta }}\\ V{p}_2=\frac{b}{U}({\Gamma }_2\pm \mathrm{i}{\mathrm{\Omega }}_2)\frac{U}{b{\omega }_{\theta }}=\frac{{\Gamma }_2}{{\omega }_{\theta }}\pm \mathrm{i}\frac{{\Omega }_2}{{\omega }_{\theta }}\end{array}\}$$

由第2节的结论知，当$\Gamma =0\Omega =0$即$p=0$时，其对应的$V=V_{\mathrm{D}}=U_{\mathrm{D}}/(b{\omega }_{\theta })$为发散速度的减缩速度。$U_{\mathrm{D}}$是发散速度。将$p=0$代入$(3.6)$式得：

$$V_{\mathrm{D}}=\frac{U_{\mathrm{D}}}{b{\omega }_{\theta }}=r\sqrt{\frac{\mu }{1+2a}}$$

同样由第2节的结论知，当$\Gamma =0\Omega \ne 0$时，其对应的$V=V_{\mathrm{F}}=U_{\mathrm{F}}/(b{\omega }_{\theta })$为颤振速度的减缩速度。$U_{\mathrm{F}}$是颤振速度。
当$a=-1/5e=-1/10\mu =20r^2=6/25\sigma =2/5$时，方程$(3.6)$对应的V-p函数图像如图所示：


可以看出无量纲颤振速度为$V_{\mathrm{F}}=U_{\mathrm{F}}/(b{\omega }_{\theta })=1.843$，颤振频率为${\Omega }_{\mathrm{F}}/{\omega }_{\theta }=0.5568$。

此外，当$V$逼近$V_{\mathrm{F}}$时俯仰与沉浮频率靠近。这条性质在颤振分析中常作为判断依据。

然而仅凭上述分析无法准确得出颤振速度，而且在低于颤振速度时，模态阻尼$\Gamma =0$不符合实际情况。同时，在颤振点上的两个根完全相等使得定常理论呈现出重合的特征。这与实验和观测均不相符合。

出现这些缺陷的原因在于使用了定常气动理论。然而非定常气动理论加入将会使得问题变为复杂。但如果预先假设运动为简谐运动，问题则可以大大简化。

4 经典颤振分析（k法-无阻尼）

根据第2节的结论知颤振边界对应的条件为：具有一个依赖于时间的简谐运动的模态。即某一阶的模态阻尼${\Gamma }_k=0$。于是在求颤振边界时，可以设第3节的模型中：

$$\begin{gathered} h(t)=\overline{h}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ \theta (t)=\overline{\theta }\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \end{gathered}$$

即，运动为简谐运动。
因为这种简谐运动的振幅比较小，于是采用线性的空气动力学理论（非定常气动理论），所以得到的升力$L$和绕$P$点的俯仰力矩$M$也是简谐的：

$$\begin{gathered} L(t)=\overline{L}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ M(t)=\overline{M}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \end{gathered}$$

其中，气动载荷的振幅可以按照运动振幅的复线性函数计算得：

$$\begin{gathered} \overline{L}=-\pi {\rho }_{\infty }{b}^3{\omega }^2[l_h(k,M{a}_{\infty })\frac{\overline{h}}{b}+l_{\theta }(k,M{a}_{\infty })\overline{\theta }] \\ \overline{M}=\pi {\rho }_{\infty }{b}^4{\omega }^2[m_h(k,M{a}_{\infty })\frac{\overline{h}}{b}+m_{\theta }(k,M{a}_{\infty })\overline{\theta }] \end{gathered}$$

其中，$l_h,l_{\theta },m_h,m_{\theta }$四个复数函数表示由沉浮运动和俯仰运动引起的无量纲气动升力系数和力矩系数。这些系数通常是参数$k$和$M{a}_{\infty }$的函数。且

$$\begin{gathered} k=\frac{b\omega }{U}(减缩频率) \\ M{a}_{\infty }=\frac{U}{c_{\infty }}(来流马赫数) \end{gathered}$$

4.1 单自由度颤振

考虑一个单自由度气动弹性系统，由刚硬的二维机翼组成，允许它绕指定的参考点俯仰运动。如图所示：

根据前面的分析：

$$\begin{gathered} \theta =\overline{\theta }\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ M=\overline{M}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ \overline{M}=\pi {\rho }_{\infty }{b}^4{\omega }^2{m}_{\theta }(k,M{a}_{\infty })\overline{\theta } \end{gathered}$$

系统的运动方程可以写作：

$$I_P\ddot{\theta }+k_{\theta }\theta =M$$

代入相关值：

$$k_{\theta }-{\omega }^2{I}_P=\pi {\rho }_{\infty }{b}^4{\omega }^2{m}_{\theta }(k,M{a}_{\infty })$$

引入系统固有频率：

$${\omega }_{\theta }=\sqrt{\frac{k_{\theta }}{I_P}}$$

可以得到颤振边界飞行条件的最终方程：

$$\frac{I_P}{\pi {\rho }_{\infty }{b}^4}[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+m_{\theta }(k,M{a}_{\infty })=0\text{\hspace{1em}(4.1.1)}$$

上述方程中$m_{\theta }(k,M{a}_{\infty })$为复数，所以可以写作：

$$\frac{I_P}{\pi {\rho }_{\infty }{b}^4}[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+\Re [m_{\theta }(k,M{a}_{\infty })]+\mathrm{i}\Im [m_{\theta }(k,M{a}_{\infty })]=0\text{\hspace{1em}(4.1.2)}$$

上述方程中有未知参数$\omega$、${\rho }_{\infty }$、$k$和$M{a}_{\infty }$，考察四个变量的相互关系，看出有$\omega$、${\rho }_{\infty }$、$U$和$c_{\infty }$四个独立变量。而方程有实部和虚部对应相等两个关系。于是只要给定其中两个量，其他的量也可以求出。

下面介绍计算方法：
首先给定高度，确定来流空气密度${\rho }_{\infty }$。将来流马赫数$M{a}_{\infty }$暂定为零（即假设来流声速$c_{\infty }$相对于空速$U$极大，设为不可压缩流）。代入$(4.1.2)$中：

$$\frac{I_P}{\pi {\rho }_{\infty }{b}^4}[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+\Re [m_{\theta }(k,0)]+\mathrm{i}\Im [m_{\theta }(k,0)]=0$$

因此可以得到左边虚部为零，实部为零，即

$$\begin{gathered} \Im [m_{\theta }(k_{\mathrm{F}},0)]=0 \\ (\frac{{\omega }_{\theta }}{{\omega }_{\mathrm{F}}}{)}^2=1+\frac{\pi {\rho }_{\infty }{b}^4\Re [m_{\theta }(k_{\mathrm{F}},0)]}{I_P} \end{gathered}$$

上面第一式可以求得颤振边界的减缩频率$k_{\mathrm{F}}$。再代入第二式可以求得颤振边界的振动频率${\omega }_{\mathrm{F}}$
由减缩频率定义式可以计算颤振速度：

$$U_{\mathrm{F}}=\frac{b{\omega }_{\mathrm{F}}}{k_{\mathrm{F}}}$$

这时给出来流声速$c_{\infty }$，进而计算得到颤振马赫数：

$$M{a}_{\mathrm{F}}=\frac{U_{\mathrm{F}}}{c_{\infty }}$$

如果此时计算得到的颤振马赫数接近之前假设的马赫数（第一次迭代为接近零），则可以取此结果。若与之前假设差距过大。则将此颤振马赫数代入$(4.1.2)$中，重复接下来的步骤。这个迭代方案将收敛于一个颤振边界上的飞行条件。最终得出此高度下的颤振速度$U_{\mathrm{F}}$和振动频率${\omega }_{\mathrm{F}}$。

4.2 二自由度颤振

这里讨论第3节的模型，如图所示：

与第3节不同的是，本节采用非定常气动理论。则其运动方程类似于$(3.4)$，为：

$$\begin{array}{c}m(\ddot{h}+b{x}_{\theta }\ddot{\theta })+k_{h}h=-L\\ I_P\ddot{\theta }+mb{x}_{\theta }\ddot{h}+k_{\theta }\theta =M\end{array}\}\text{\hspace{1em}(4.2.1)}$$

根据第4节开始的分析：

$$\begin{gathered} h(t)=\overline{h}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ \theta (t)=\overline{\theta }\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ L(t)=\overline{L}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ M(t)=\overline{M}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ \overline{L}=-\pi {\rho }_{\infty }{b}^3{\omega }^2[l_h(k,M{a}_{\infty })\frac{\overline{h}}{b}+l_{\theta }(k,M{a}_{\infty })\overline{\theta }] \\ \overline{M}=\pi {\rho }_{\infty }{b}^4{\omega }^2[m_h(k,M{a}_{\infty })\frac{\overline{h}}{b}+m_{\theta }(k,M{a}_{\infty })\overline{\theta }] \end{gathered}$$

将上面的函数代入到运动方程$(4.2.1)$中，得：

$$\begin{gathered} \{\frac{m}{\pi {\rho }_{\infty }{b}^2}[1-(\frac{{\omega }_h}{\omega }{)}^2]+l_h(k,M{a}_{\infty })\}\frac{\overline{h}}{b}+[\frac{m{x}_{\theta }}{\pi {\rho }_{\infty }{b}^2}+l_{\theta }(k,M{a}_{\infty })]\overline{\theta }=0 \\ [\frac{m{x}_{\theta }}{\pi {\rho }_{\infty }{b}^2}+m_h(k,M{a}_{\infty })]\frac{\overline{h}}{b}+\{\frac{I_P}{\pi {\rho }_{\infty }{b}^4}[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+m_{\theta }(k,M{a}_{\infty })\}\overline{\theta }=0 \end{gathered}$$

引入无量纲参数：

$$\begin{gathered} \mu =\frac{m}{\pi {\rho }_{\infty }{b}^2}(质量比) \\ r=\sqrt{\frac{I_P}{m{b}^2}}(关于P点的回转质量半径) \end{gathered}$$

可以简化为：

$$\begin{array}{c}\{\mu [1-(\frac{{\omega }_h}{\omega }{)}^2]+l_h\}\frac{\overline{h}}{b}+(\mu x_{\theta }+l_{\theta })\overline{\theta }=0\\ (\mu x_{\theta }+m_h)\frac{\overline{h}}{b}+\{\mu r^2[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+m_{\theta }\}\overline{\theta }=0\end{array}\}\text{\hspace{1em}(4.2.2)}$$

要使方程有非平凡解，系数行列式为零。同时引入参数$\sigma ={\omega }_h/{\omega }_{\theta }$和$\lambda =({\omega }_{\theta }/\omega {)}^2$，即：

$$\left|\begin{array}{cc}\mu (1-{\sigma }^2\lambda )+l_h& \mu x_{\theta }+l_{\theta }\\ \mu x_{\theta }+m_h& \mu r^2(1-\lambda )+m_{\theta }\end{array}\right|=0\text{\hspace{1em}(4.2.3)}$$

这个行列式被称为“颤振行列式”。其中未知量有$\lambda =({\omega }_{\theta }/\omega {)}^2$、$\mu =m/(\pi {\rho }_{\infty }{b}^2)$、$M{a}_{\infty }=U/c_{\infty }$和$k=b\omega /U$。同样有$\omega$、${\rho }_{\infty }$、$U$和$c_{\infty }$四个独立变量，且上式有虚部、实部分别等于零两个关系。于是只需要确定其中两个变量，其它变量也可以求得。

下面介绍求解颤振边界的主要步骤：

1. 给定一个高度（即给定来流密度${\rho }_{\infty }$），则参数$\mu$确定。
2. 给定$M{a}_{\infty }$的初始值，例如零（即假设来流声速$c_{\infty }$相对于空速$U$极大，设为不可压缩流）。
3. 令颤振行列式等于0，得到关于$\lambda$的二次方程。由虚部为零得$k$值，即$k_{\mathrm{F}}$。
4. 将$k=k_{\mathrm{F}}$代入上面二次方程的实部。得到$\lambda$的值，即${\lambda }_{\mathrm{F}}=({\omega }_{\theta }/{\omega }_{\mathrm{F}}{)}^2$。
5. 给定来流声速$c_{\infty }$。确定$U_{\mathrm{F}}=b{\omega }_{\mathrm{F}}/k_{\mathrm{F}}$和$M{a}_{\infty \mathrm{F}}=U_{\mathrm{F}}/c_{\infty }$。
6. 使用步骤5计算的数值$M{a}_{\infty \mathrm{F}}$重复步骤3~5，直至在给定的$\mu$值下，$M{a}_{\infty \mathrm{F}}$、$k_{\mathrm{F}}$和$U_{\mathrm{F}}$的值收敛为止。
7. 取不同的$\mu$值（即不同高度），重复整个过程确定颤振边界，即给出高度与$M{a}_{\infty \mathrm{F}}$、$k_{\mathrm{F}}$及$U_{\mathrm{F}}$的变化曲线。

5 颤振问题的工程求解（k法、p-k法）

工程实践中除了求出颤振边界，还有两个着重考虑的问题。第一个，确定在颤振边界附近飞行条件下的稳定裕度；第二个，确定导致不稳定的物理机制。

5.1 k法

在4节的基础上，再引入表征结构阻尼效应的参数。实验观测，在简谐振动中，每个周期消耗的能量与振幅的平方近似成比例，而与频率无关，故可以通过一个与位移成比例、与速度同相位的阻尼力来体现这个特征。

引入阻尼力后运动方程$(4.2.1)$修正为：

$$\begin{array}{c}m(\ddot{h}+b{x}_{\theta }\ddot{\theta })+k_{h}h=-L+D_h\\ I_P\ddot{\theta }+mb{x}_{\theta }\ddot{h}+k_{\theta }\theta =M+D_{\theta }\end{array}\}$$

其中耗散阻尼项为：

$$\begin{gathered} D_h=\overline{D_h}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t)=-\mathrm{i}{g}_{h}m{\omega }_h^2\overline{h}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \\ {D}_{\theta }=\overline{D_{\theta }}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t)=-\mathrm{i}{g}_{\theta }{I}_P{\omega }_{\theta }^2\overline{\theta }\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t) \end{gathered}$$

于是$(4.2.2)$修正为：

$$\begin{array}{c}\{\mu [1-(\frac{{\omega }_h}{\omega }{)}^2(1+\mathrm{i}{g}_h)]+l_h\}\frac{\overline{h}}{b}+(\mu x_{\theta }+l_{\theta })\overline{\theta }=0\\ (\mu x_{\theta }+m_h)\frac{\overline{h}}{b}+\{\mu r^2[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2(1+\mathrm{i}{g}_{\theta })]+m_{\theta }\}\overline{\theta }=0\end{array}\}\text{\hspace{1em}(5.1.1)}$$

注意$g_h$和$g_{\theta }$为人为引入修正参数。可以凭经验预先指定，典型的取值范围为0.01~0.05。这里令$g_h=g_{\theta }=g$，并将$\omega$和$g$合并为一个参数$Z$：

$$Z=(\frac{{\omega }_{\theta }}{\omega }{)}^2(1+\mathrm{i}g)$$

同样记$\sigma ={\omega }_h/{\omega }_{\theta }$。可得颤振行列式：

$$\left|\begin{array}{cc}\mu (1-{\sigma }^{2}Z)+l_h& \mu x_{\theta }+l_{\theta }\\ \mu x_{\theta }+m_h& \mu r^2(1-Z)+m_{\theta }\end{array}\right|=0\text{\hspace{1em}(5.1.2)}$$

此方程中有独立未知变量$\omega$、$g$、${\rho }_{\infty }$、$U$和$c_{\infty }$。而有实部和虚部两个等式关系。可以给定${\rho }_{\infty }$（即高度$H$）和$c_{\infty }$的值，得出函数$g=g(U)$和$\omega =\omega (U)$，绘制出函数图像U-g和U-ω。

其中，U-g图像显示出颤振边界附近的飞行条件下的稳定裕度，当$g=0$时为颤振边界。结合U-ω图像可以得出导致颤振不稳定的物理机理。沿着$U=0$轴上的频率值对应于初始结构动力学系统耦合模态。随着速度的增加，这些根的个体特征或者他们之间的相互影响体现了能量从一个模态到另一个模态的传递。

5.2 p-k法

5.2.1 p法和k法的比较

现在研究第2节的模型，即：将升力面简化为板模型。

将$(2.3)$式中的$\overline{\mathbf{\lambda }}$消去，得：

$$[p^2\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})+\frac{b^2}{U^2}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\omega }}^2)-{\rho }_{\infty }\mathcal{A}(p)]\overline{\mathbf{\xi }}=0\text{\hspace{1em}(5.2.1.1)}$$

其中$\mathcal{A}(p)$为非定常气动算子矩阵，可用$(2.3)$式中其他矩阵（即$\mathbf{a}$、$\mathbf{b}$、$\mathbf{c}$、$\mathbf{d}$、$\mathbf{A}$和$\mathbf{E}$）表示。

广义振幅要有非平凡解，上式的系数行列式必须为零，即

$$\left|\begin{array}{c}{p}^2\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})+\frac{b^2}{U^2}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\omega }}^2)-{\rho }_{\infty }\mathcal{A}(p)\end{array}\right|=0\text{\hspace{1em}(5.2.1.2)}$$

由于$(2.3)$和$(5.2.1.1)$是等价的，所以$(5.2.1.2)$和$(2.4)$同解。所以$(5.2.1.2)$的根为：

$$p_l=\frac{bv}{U}=\frac{b}{U}({\Gamma }_l\pm \mathrm{i}{\Omega }_l)(l=0,1,\cdots ,n_c,\cdots ,n_c+n_r)\text{\hspace{1em}(5.2.1.3)}$$

而第4节和第5节的k法本质上将$v$设为$\mathrm{i}\omega$，且$k=b\omega /U$。而k法首先是假定${\Gamma }_l=0$，所以$k=\Omega b/U$。

当$\Omega \ne 0$时各典型模态幅值的状态如图

注意，图中$\mathrm{e}\mathrm{x}\mathrm{p}[({\Gamma }_k+\mathrm{i}{\Omega }_k)t]$表示的是由两个共轭复根$\mathrm{e}\mathrm{x}\mathrm{p}[({\Gamma }_k\pm \mathrm{i}{\Omega }_k)t]$导出的一对等价实根$\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t)$和$\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t)$中的前者。即：

$$\begin{gathered} \mathrm{e}\mathrm{x}\mathrm{p}[({\Gamma }_k\pm \mathrm{i}{\Omega }_k)t]=\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)[\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t)\pm \mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t)] \\ \frac{\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)[\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t)+\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t)]+\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)[\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t)-\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t)]}{2}=\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t) \\ \frac{\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)[\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t)+\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t)]-\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)[\mathrm{c}\mathrm{o}\mathrm{s}({\Omega }_{k}t)-\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t)]}{2\mathrm{i}}=\mathrm{e}\mathrm{x}\mathrm{p}({\Gamma }_{k}t)\mathrm{s}\mathrm{i}\mathrm{n}({\Omega }_{k}t) \end{gathered}$$

假设$t_0$和$t_0+T$为图中某两个相邻波峰的时间（$\mathrm{c}\mathrm{o}\mathrm{s}(\Omega t_0)=1$、$\mathrm{c}\mathrm{o}\mathrm{s}[\Omega (t_0+T)]=1$）。设$a_m$和$a_{m+1}$为连续周期的振幅。则：

$$\begin{gathered} a_m=\mathrm{e}\mathrm{x}\mathrm{p}(\Gamma t_0)\mathrm{c}\mathrm{o}\mathrm{s}(\Omega t_0)=\mathrm{e}\mathrm{x}\mathrm{p}(\Gamma t_0) \\ {a}_{m+1}=\mathrm{e}\mathrm{x}\mathrm{p}[\Gamma (t_0+T)]\mathrm{c}\mathrm{o}\mathrm{s}[\Omega (t_0+T)]=\mathrm{e}\mathrm{x}\mathrm{p}[\Gamma (t_0+T)] \end{gathered}$$

则：

$$\mathrm{l}\mathrm{n}(\frac{a_{m+1}}{a_m})=\Gamma T$$

又周期$T=2\pi /\Omega$，所以

$$\Gamma =\frac{1}{2\pi }\mathrm{l}\mathrm{n}(\frac{a_{m+1}}{a_m})\Omega$$

令$\gamma =\mathrm{l}\mathrm{n}(a_{m+1}/a_m)/(2\pi )$，命名为衰减率。则

$$p_l=\frac{bv}{U}=\frac{b}{U}({\Gamma }_l\pm \mathrm{i}{\Omega }_l)=\frac{b}{U}{\Gamma }_l\pm \mathrm{i}\frac{b}{U}{\Omega }_l=\gamma \frac{b}{U}\Omega \pm \mathrm{i}\frac{b}{U}{\Omega }_l$$

已推得$k=\Omega b/U$。所以：

$$p=\gamma k\pm \mathrm{i}k\text{\hspace{1em}(5.2.1.4)}$$

于是令$(5.2.1.3)$中$p=\mathrm{i}k$，可得颤振行列式：

$$\left|\begin{array}{c}-k^2\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})+\frac{b^2}{U^2}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{M})\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{\omega }}^2)-{\rho }_{\infty }\mathcal{A}(\mathrm{i}k)\end{array}\right|=0$$

其中独立未知变量有$\Omega$、$U$、${\rho }_{\infty }$、$g$（k法会引入一个阻尼参数）和$c_{\infty }$，而有实部和虚部两个关系式。于是在给定${\rho }_{\infty }$（高度）和$c_{\infty }$时得出函数$g=g(U)$和$\Omega =\Omega (U)$，绘制出函数图像U-g和U-Ω。

经过分析，阻尼参数$g$和衰减率$\gamma$关系为：

$$g\cong 2\gamma$$

可以看出，k法在分析颤振附近的飞行条件时，将原本$p$中的阻尼项$\gamma k$略去，而人为加上了一个阻尼参数$g$进行修正。然而这个修正并不能完全符合原方程关系。如图所示为p法和k法计算相同的模型得出的曲线。

k法不仅得出错误的模态耦合，而且更重要的是，得出了错误的不稳定模态。这两种分析唯一一致有效的结果是在$g=\gamma =0$点处所对应的颤振速度。

尽管k法给出了不一致的模态耦合，但是它可以使用非定常气动力的简谐模型。简谐运动载荷的计算精度优于瞬态运动气动载荷的计算精度。于是产生了两种方法的折中。

5.2.2 p-k法简介

p-k法就是这样一种折中。它基于p法进行颤振分析，但使用简谐运动的非定常气动力矩阵。使用$k$值来计算矩阵$\mathcal{A}(\mathrm{i}k)$，颤振行列式为

$$\left|\begin{array}{c}{p}^2\mathbf{M}+\frac{b^2}{U^2}\mathbf{M}{\mathbf{\omega }}^2-{\rho }_{\infty }\mathcal{A}(\mathrm{i}k)\end{array}\right|=0$$

求解方法如下：

1. 给定一组$k$的初始估计值，比如对于第$i$个根，$k_0=b{\omega }_i/U$。
2. 将$k_0$代入方程中，解出$p$。
3. 由于$p$是简单出现在方程中的，所以以上方程可以看作关于$p$的标准特征值问题。典型解是一组共轭复根对。选择其中的一个根，初始解记为：

$$k_1=|\Im (p)|,{\gamma }_1=\frac{\Re (p)}{k_1}$$

1. 将新解得的$k_1$代入方程中，重复步骤2~3，直到每一个根收敛。

对于低阶问题的求根过程，可以取$k=|\Im (p)|$时的行列式为零。

用p法和p-k法计算5.2.1节的模型，得曲线如图所示：

p-k法得到了于p法几乎相同的结果。p-k法的最大优势在于它可以使用简谐运动形式的气动载荷。

将k法和p-k法用在水平安定面/升降舵构型的颤振分析中。这是一个强耦合的系统。分析结果如图所示：

这里也得到了大相径庭的模态耦合结论。p-k法除了对强耦合系统很容易分析频率-速度和阻尼-速度曲线外，它的第二个优势在于计算耗费。在密度不变的情况下，k法需要大量的运算，以确保马赫数与速度和高度匹配，而p-k法则不需要。

6 非定常气动力

在第3节使用定常气动理论进行颤振分析时，升力和俯仰力矩仅是瞬时俯仰角$\theta$的函数。然而深入探究，很容易看出攻角并不简单地就等于$\theta$。例如，翼型参考基准点以速度$\dot{h}$做沉浮运动，可以论证至少对于小角度，应修正攻角以变引入沉浮影响，即

$$\alpha =\theta +\frac{\dot{h}}{U}$$

这源于考虑飞机滚转对机翼攻角的影响。但必须小心做此特定的推论，因为可能还有其他同等量级的影响没有考虑。

在非定常流场中，升力和俯仰力矩由两部分构成——非环量效应和环量效应。
环量效应通常对机翼更加重要。事实上，在定常飞行中，是环量产生的升力使飞机在空中飞行。回顾Helmholtz定理：在任意包围特定流体微团的封闭曲线内，总涡量始终为零。因此，如果一个顺时针涡量随机翼发展变化，那么一定有相同强度的逆时针漩涡脱落到流动中。随着脱落涡运动、发展，它会诱导产生非定常流动，从而改变流场，反过来影响机翼。在估算这些脱落涡的影响时，可以假设脱落涡随来流一起运动，简化运算。
非环量效应，也称为视在质量和惯性效应，重要性次于环量效应。当机翼有加速度时会产生非环量效应，而且会因此带动部分周围的空气随之运动，这部分空气具有确定的质量，会产生与其加速度对应的惯性力。

总而言之，非定常气动理论需要考虑至少以下三个独立的物理现象：

1. 风速方向的变化改变了有效攻角，从而改变了升力。
2. 机翼的运动对气流有扰动，并使漩涡自后缘脱落。脱落涡产生的下洗，反过来会改变机翼的绕流。这个非定常下洗改变了有效攻角，从而改变了升力。
3. 机翼的运动会使翼表面附近的空气微团加速运动，因此需要考虑由此产生的惯性力。

6.1 Theodorsen非定常薄翼理论

Theodorsen（1934年）推导了小简谐振动情况下不可压流的薄翼（平板）非定常气动理论。

根据Theodorsen理论，升力和俯仰力矩可表示为

$$\begin{array}{c}L=2\pi {\rho }_{\infty }UbC(k)[\dot{h}+U\theta +b(\frac{1}{2}-a)\dot{\theta }]+\pi {\rho }_{\infty }{b}^2(\ddot{h}+U\dot{\theta }-ba\ddot{\theta })\\ M_{\frac{1}{4}}=-\pi {\rho }_{\infty }{b}^3[\frac{1}{2}\ddot{h}+U\dot{\theta }+b(\frac{1}{8}-\frac{a}{2})\ddot{\theta }]\end{array}\}\text{\hspace{1em}(6.1.1)}$$

函数$C(k)$是减缩频率$k$的复数函数，可表示为

$$C(k)=\frac{H_1^{(1)}(k)}{H_1^{(2)}(k)+\mathrm{i}{H}_0^{(2)}(k)}$$

其中$H_n^{(2)}(k)$是第二类Hankel函数，可以分别由第一类和第二类Bessel函数表示为

$$H_n^{(2)}(k)={\mathrm{J}}_n(k)-\mathrm{i}{Y}_n(k)$$

函数$C(k)=F(k)+\mathrm{i}G(k)$称为Theodorsen函数。$k$从0变化到1时$C(k)$的实部和虚部曲线如图所示：

注意对于定常情况（即$k=0$），而且等于$1$。随着$k$增加，虚部也增加，实部减小；当$k$趋于无穷大时，$C(k)$接近于$1/2$。但在实际情况中，$k$不会超过$1$，所以图中只绘制了直到$k=1$的情况。任何一个调和函数乘以$C(k)$后，数值都会减小，同时都会产生相位滞后。

下面对升力和俯仰力矩表达式简要分析介绍：
Theodorsen理论推导基于线性势流理论，升力线斜率等于$2\pi$。
升力表达式的第一项式环量升力（不考虑脱落效应）乘以$C(k)$，这是为了考虑脱落涡的影响。
非环量项（即升力表达式第二项和全部俯仰力矩）取决于机翼的加速度和角加速度。注意到，升力中$\ddot{h}$的系数是以$b$为半径的单位长度的圆柱（无限长圆柱）所含的空气质量，它反映了因机翼运动使多少空气获得加速度。
对于定常流，环量升力与攻角呈线性关系。对于非定常流和简谐运动，可以从Theodorsen理论中推出有效攻角：

$$\alpha =C(k)[\theta +\frac{\dot{h}}{U}+\frac{b}{U}(\frac{1}{2}-a)\dot{\theta }]$$

通过与有限状态气动力模型比较表明，这是3/4弦线处的攻角。

注意：一般而言Theodorsen理论只在简谐振动条件下才有效，是频域微分方程。
但当$C(k)$等于1时，Theodorsen理论的近似称为“准定常”薄翼理论。这样的近似仅当$k$取值非常小时才适用。对于缓慢的简谐运动或缓慢的非简谐的变化运动，准定常理论可以用于时域。

Theodorsen理论可以用于k法和p-k法。

6.2 Peter等人的有限状态非定常薄翼理论

在进行颤振分析时，通常需要计算亚临界飞行条件下的模态阻尼；在颤振主动控制设计中，控制器的设计要求系统可表示为状态空间形式。为了满足这些要求，需要用时域微分方程来表示真实的气动载荷（Theodorsen的理论中为频域）。
有限状态理论是在工程精度范围内对真实的无线状态的空气动力模型进行近似。其中一个方法是Peter等人（1955年）提出的针对无粘不可压流的有限状态诱导流理论。

考虑一个刚硬、对称机翼的典型剖面，如图所示。

附加矢量方向定义。如图所示：

其中有三组单位矢量：

1. 固定于惯性系的一组矢量${\mathbf{i}}_1$和${\mathbf{i}}_2$，则空气流动速度为$-U{\mathbf{i}}_1$；
2. 固定于机翼的一组矢量${\mathbf{b}}_1$和${\mathbf{b}}_2$，其中${\mathbf{b}}_1$沿着零升力线指向前缘，${\mathbf{b}}_2$垂直于${\mathbf{b}}_1$；
3. 与3/4弦线处（T点）的当地相对风速矢量相关的一组矢量${\mathbf{a}}_1$和${\mathbf{a}}_2$，${\mathbf{a}}_1$沿着相对风速矢量的方向，${\mathbf{a}}_1$和${\mathbf{a}}_2$垂直，沿着假定的升力方向。

这些单位矢量之间的关系可以表示为：

$$\begin{gathered} \left[\begin{array}{c}{\mathbf{b}}_1\\ {\mathbf{b}}_2\end{array}\right]=\left[\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{s}\mathrm{i}\mathrm{n}\theta \\ -\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right]\left[\begin{array}{c}{\mathbf{i}}_1\\ {\mathbf{i}}_2\end{array}\right],\left[\begin{array}{c}{\mathbf{a}}_1\\ {\mathbf{a}}_2\end{array}\right]=\left[\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\theta & -\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ \mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right]\left[\begin{array}{c}{\mathbf{b}}_1\\ {\mathbf{b}}_2\end{array}\right], \\ {\mathbf{i}}_3={\mathbf{a}}_3={\mathbf{b}}_3={\mathbf{b}}_1\times {\mathbf{b}}_2 \end{gathered}$$

诱导流理论以脱落涡引起的翼型附近流场的改变为基础，对脱落涡的影响进行了近似。因此，翼型附近的速度场包含自由流速度和用来描述诱导流的额外速度分量。并将诱导流近似为弦线方向上的平均值。以此当地的惯性风速可以写为$-U{\mathbf{i}}_1-{\lambda }_0{\mathbf{b}}_2$，这里${\lambda }_0$是平均诱导流速度（垂直于翼型的零升力线）。根据经典薄翼理论，应该使用T点的瞬时相对风速来计算攻角$\alpha$。T点相对风速$W{\mathbf{a}}_1$可以写作：

$$W{\mathbf{a}}_1={\mathbf{v}}_T-(-U{\mathbf{i}}_1-{\lambda }_0{\mathbf{b}}_2)={\mathbf{v}}_T+U{\mathbf{i}}_1+{\lambda }_0{\mathbf{b}}_2$$

其中${\mathbf{v}}_T$是3/4弦长处的惯性速度，表达式为

$${\mathbf{v}}_T={\mathbf{v}}_P+\dot{\theta }{\mathbf{b}}_3\times {\mathbf{r}}_{PT}$$

${\mathbf{r}}_{PT}$是从P点到T点位置矢量。从图中可以看出：

$${\mathbf{r}}_{PT}=[\frac{b}{2}+(1+a)b-2b]{\mathbf{b}}_1=b(a-\frac{1}{2}){\mathbf{b}}_1$$

${\mathbf{v}}_P$是P点的惯性速度：

$${\mathbf{v}}_P=-\dot{h}{\mathbf{i}}_2$$

$\dot{\theta }{\mathbf{b}}_3$为机翼的惯性角速度。
所以

$${\mathbf{v}}_T=-\dot{h}{\mathbf{i}}_2+b\dot{\theta }(a-\frac{1}{2}){\mathbf{b}}_2$$

代入相对风速的表达式：

$$W{\mathbf{a}}_1=U{\mathbf{i}}_1-\dot{h}{\mathbf{i}}_2+[b\dot{\theta }(a-\frac{1}{2})+{\lambda }_0]{\mathbf{b}}_2$$

于是：

$$\begin{array}{c}W{\mathbf{a}}_1\cdot {\mathbf{b}}_1=U\mathrm{c}\mathrm{o}\mathrm{s}\theta -\dot{h}\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ W{\mathbf{a}}_1\cdot {\mathbf{b}}_2=-U\mathrm{s}\mathrm{i}\mathrm{n}\theta -\dot{h}\mathrm{c}\mathrm{o}\mathrm{s}\theta +b\dot{\theta }(a-\frac{1}{2})+{\lambda }_0\end{array}\}\text{\hspace{1em}(6.2.1)}$$

另外,由矢量关系：

$$W{\mathbf{a}}_1=W\mathrm{c}\mathrm{o}\mathrm{s}\alpha {\mathbf{b}}_1-W\mathrm{s}\mathrm{i}\mathrm{n}\alpha {\mathbf{b}}_2$$

于是：

$$\begin{array}{c}W{\mathbf{a}}_1\cdot {\mathbf{b}}_1=W\mathrm{c}\mathrm{o}\mathrm{s}\alpha \\ W{\mathbf{a}}_1\cdot {\mathbf{b}}_2=-W\mathrm{s}\mathrm{i}\mathrm{n}\alpha \end{array}\}\text{\hspace{1em}(6.2.2)}$$

对比$(6.2.1)$和$(6.2.2)$得：

$$\begin{gathered} W\mathrm{c}\mathrm{o}\mathrm{s}\alpha =U\mathrm{c}\mathrm{o}\mathrm{s}\theta -\dot{h}\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ -W\mathrm{s}\mathrm{i}\mathrm{n}\alpha =-U\mathrm{s}\mathrm{i}\mathrm{n}\theta -\dot{h}\mathrm{c}\mathrm{o}\mathrm{s}\theta +b\dot{\theta }(a-\frac{1}{2})+{\lambda }_0 \end{gathered}$$

$\theta$和$\alpha$为小角度，则可设：

$$\begin{gathered} \mathrm{s}\mathrm{i}\mathrm{n}\theta =\theta \\ \mathrm{s}\mathrm{i}\mathrm{n}\alpha =\alpha \\ \mathrm{c}\mathrm{o}\mathrm{s}\theta =1 \\ \mathrm{s}\mathrm{i}\mathrm{n}\alpha =1 \\ \alpha \theta =0 \end{gathered}$$

上式可变为：

$$\begin{gathered} W=U-\dot{h}\theta \\ -W\alpha =-U\theta -\dot{h}+b\dot{\theta }(a-\frac{1}{2})+{\lambda }_0 \end{gathered}$$

解上面方程，注意$\alpha \theta =0$：

$$\begin{array}{c}\alpha =\theta +\frac{\dot{h}}{U}+\frac{b}{U}(\frac{1}{2}-a)\dot{\theta }-\frac{{\lambda }_0}{U}\\ W=U-\dot{h}\theta \end{array}\}\text{\hspace{1em}(6.2.3)}$$

可见，攻角$\alpha$取受俯仰角速度、沉浮速度和诱导流影响。

上述考虑了风速方向改变和漩涡脱落两个因素（即环量效应的因素）改变了有效攻角从而引起的非定常特性。除此之外还有非环量效应因素引起的非定常特性。如果将三者共同考虑，将会升力和力矩表达式：

$$\begin{array}{c}L=2\pi {\rho }_{\infty }Ub[\dot{h}+U\theta +b(\frac{1}{2}-a)\dot{\theta }-{\lambda }_0]+\pi {\rho }_{\infty }{b}^2(\ddot{h}+U\dot{\theta }-ba\ddot{\theta })\\ M_{\frac{1}{4}}=-\pi {\rho }_{\infty }{b}^3[\frac{1}{2}\ddot{h}+U\dot{\theta }+b(\frac{1}{8}-\frac{a}{2})\ddot{\theta }]\end{array}\}\text{\hspace{1em}(6.2.4)}$$

但其中有一个未确定量${\lambda }_0$，即需要将${\lambda }_0$用机翼运动表示。Peter等人的诱导理论的做法是：用$N$个诱导流速度${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_N$来表示平均诱导速度${\lambda }_0$。如：

$${\lambda }_0\approx \frac{1}{2}\sum _{n=1}^N{p}_n{\lambda }_n$$

其中$p_n$由最小二乘法求得。$N$是由计算者给定。
假设脱落涡保持在机翼平面，并且以同样的速度随气流顺流而下，由此可以推导出诱导流的动力学方程：

$$\mathbf{A}\dot{\mathbf{\lambda }}+\frac{U}{b}\mathbf{\lambda }=\mathbf{c}[\ddot{h}+U\dot{\theta }+b(\frac{1}{2}-a)\ddot{\theta }]$$

其中：

$$\begin{gathered} \mathbf{\lambda }={\left(\begin{array}{cccc}{\lambda }_1& {\lambda }_2& \cdots & {\lambda }_N\end{array}\right)}^{\mathrm{T}} \\ \mathbf{A}=\mathbf{D}+\mathbf{d}{\mathbf{b}}^{\mathrm{T}}+\mathbf{c}{\mathbf{d}}^{\mathrm{T}}+\frac{1}{2}\mathbf{c}{\mathbf{b}}^{\mathrm{T}} \\ {D}_{nm}=\{\begin{array}{l}\frac{1}{2n},(n=m+1)\\ -\frac{1}{2n},(n=m+1)\\ 0,(n\ne m\pm 1)\end{array} \\ {b}_n=\{\begin{array}{l}(-1{)}^{n-1}\frac{(N+n-1)!}{(N-n-1)!}{\frac{1}{n!}}^2,(n\ne N)\\ (-1{)}^{n-1},(n=N)\end{array} \\ {d}_n=\{\begin{array}{l}\frac{1}{2},(n=1)\\ 0,(n\ne 1)\end{array} \\ {c}_n=\frac{2}{n} \end{gathered}$$

求解上面的一阶常微分方程就可以得到${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_N$的表达式，进而求得${\lambda }_0$的表达式。

将本小节的有限状态非定常薄翼理论引入第3节的模型。即将$(6.2.4)$式代入$(3.4)$式中，对第3节中的算例（$a=-1/5e=-1/10\mu =20r^2=6/25\sigma =2/5$）进行结果分析，分析时选取$N=6$个诱导流速度。频率和阻尼分析结果如图所示：


与之前类似，在接近不稳定时，可以观察到有频率相互靠近的趋势，但是颤振条件式以某一个根的实部穿越过横轴成为正值来确定的，从而得到颤振速度$V_{\mathrm{F}}=U_{\mathrm{F}}/(b{\omega }_{\theta })=2.165$，颤振频率${\Omega }_{\mathrm{F}}/{\omega }_{\theta }=0.6545$。

7 根据假设模态进行颤振预测

考虑一个安装在风洞壁上的非后掠机翼，用长为$l$的均匀悬臂梁模拟，如图所示：

因此对于一个具有弯曲刚度为$\overline{EI}$和扭转刚度$\overline{GJ}$的梁，其应变能为：

$$U=\frac{1}{2}{\int }_0^l[\overline{EI}(\frac{{\partial }^{2}w}{\partial y^2}{)}^2+\overline{GJ}(\frac{\delta \theta }{\delta y}{)}^2]\mathrm{d}y\text{\hspace{1em}(7.1)}$$

其翼剖面如图所示：

材料密度为$\rho$，横截面内典型点的速度为

$$\mathbf{v}=z\frac{\partial \theta }{\partial t}\mathbf{i}+(\frac{\partial w}{\partial t}-x\frac{\partial \theta }{\partial t})\mathbf{k}$$

$\mathbf{i}$和$\mathbf{k}$分别是x和y方向上的单位矢量，则动能可表示为

$$\begin{array}{rl}K& =\frac{1}{2}{\int }_0^l{\iint }_A\rho [(\frac{\partial w}{\partial t}-x\frac{\partial \theta }{\partial t}{)}^2+z^2(\frac{\partial \theta }{\partial t}{)}^2]\mathrm{d}x\mathrm{d}y\mathrm{d}z\\ & =\frac{1}{2}{\int }_0^l[m(\frac{\partial w}{\partial t}{)}^2+2md\frac{\partial w}{\partial t}\frac{\partial \theta }{\partial t}+m{b}^2{r}^2(\frac{\partial \theta }{\partial t}{)}^2]\mathrm{d}y\end{array}$$

其中，$m$是单位长度质量，$d$是质心与弹性轴之间的距离（当质心朝向前缘为正），$b$是半弦长，$br$是截面质量绕弹性轴的回转半径。

气动力所做的虚功为：

$$\overline{\delta W}={\int }_0^l[L^{\prime }\delta w+(M_{ac}^{\prime }+e{L}^{\prime })\delta \theta ]\mathrm{d}y$$

上面使用的符号是静气动力学的符号习惯。现在改为非定常气动力学符号习惯。

$$\begin{gathered} d\to -b{x}_{\theta } \\ e\to (\frac{1}{2}+a)b \\ {L}^{\prime }\to L^{\prime } \\ {M}_{ac}^{\prime }\to M_{\frac{1}{4}}^{\prime } \end{gathered}$$

动能变为：

$$K=\frac{1}{2}{\int }_0^l[m(\frac{\partial w}{\partial t}{)}^2-2mb{x}_{\theta }\frac{\partial w}{\partial t}\frac{\partial \theta }{\partial t}+m{b}^2{r}^2(\frac{\partial \theta }{\partial t}{)}^2]\mathrm{d}y\text{\hspace{1em}(7.2)}$$

虚功变为：

$$\overline{\delta W}={\int }_0^l\{L^{\prime }\delta w+[M_{\frac{1}{4}}^{\prime }+(\frac{1}{2}+a)b{L}^{\prime }]\delta \theta \}\mathrm{d}y\text{\hspace{1em}(7.3)}$$

合理的假设模态是选择一组不耦合的悬臂梁的自由弯曲和和扭转模态：

$$\begin{array}{c}w(y,t)={\sum }_{i=1}^{N_w}{\eta }_i(t){\Psi }_i(y)\\ \theta (y,t)={\sum }_{i=1}^{N_{\theta }}{\varphi }_i(t){\Theta }_i(y)\end{array}\}\text{\hspace{1em}(7.4)}$$

其中，$N_w$和$N_{\theta }$分别表示弯曲模态和扭转模态的个数，${\eta }_i$和${\varphi }_i$分别是与弯曲和扭转相关的广义坐标，${\Psi }_i$和${\Theta }_i$分别是弯曲和扭转模态振型。且：

$$\begin{gathered} {\Theta }_i=\sqrt{2}\mathrm{s}\mathrm{i}\mathrm{n}({\gamma }_{i}y) \\ {\Psi }_i=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}({\alpha }_{i}y)-\mathrm{c}\mathrm{o}\mathrm{s}({\alpha }_{i}y)-{\beta }_i[\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}({\alpha }_{i}y)-\mathrm{s}\mathrm{i}\mathrm{n}({\alpha }_{i}y)] \end{gathered}$$

其中，${\gamma }_i=\pi (i-1/2)/l$，而${\alpha }_i$和${\beta }_i$满足：

$$\begin{gathered} \mathrm{c}\mathrm{o}\mathrm{s}(\alpha l)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(\alpha l)+1=0 \\ {\beta }_i=\frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}({\alpha }_{i}l)+\mathrm{c}\mathrm{o}\mathrm{s}({\alpha }_{i}l)}{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}({\alpha }_{i}l)+\mathrm{s}\mathrm{i}\mathrm{n}({\alpha }_{i}l)} \end{gathered}$$

将$(7.4)$代入$(7.1)$和$(7.2)$中，利用弯曲和扭转二者模态存在正交性可简化得：

$$U=\frac{1}{2}[\frac{\overline{EI}}{l^3}\sum _{i=1}^{N_w}({\alpha }_{i}l{)}^4{\dot{\eta }}_i^2+\frac{\overline{GJ}}{l}\sum _{i=1}^{N_{\theta }}({\gamma }_{i}l{)}^2{\varphi }_i^2]\text{\hspace{1em}(7.5)}$$

$$K=\frac{ml}{2}(\sum _{i=1}^{N_w}{\dot{\eta }}_i^2+b^2{r}^2\sum _{i=1}^{N_{\theta }}{\dot{\varphi }}_i^2-2b{x}_{\theta }\sum _{i=1}^{N_{\theta }}\sum _{i=1}^{N_w}{A}_{ij}\dot{{\varphi }_i}\dot{{\eta }_i})\text{\hspace{1em}(7.6)}$$

其中

$$A_{ij}=\frac{1}{l}{\int }_0^l{\Theta }_i{\Psi }_j\mathrm{d}y(i=1,2,\cdots ,N_{\theta };j=1,2,\cdots ,N_w)$$

将$(7.4)$代入$(7.3)$得

$$\overline{\delta W}=\sum _{i=1}^{N_w}\delta {\eta }_i{\int }_0^l{\Psi }_i{L}^{\prime }\mathrm{d}y+\sum _{i=1}^{N_{\theta }}\delta {\varphi }_i{\int }_0^l{\Theta }_i[M_{\frac{1}{4}}^{\prime }+(\frac{1}{2}+a)b{L}^{\prime }]\mathrm{d}y$$

于是广义力为

$$\begin{array}{c}{\Xi }_{w_i}={\int }_0^l{\Psi }_i{L}^{\prime }\mathrm{d}y\\ {\Xi }_{{\theta }_i}={\int }_0^l{\Theta }_i[M_{\frac{1}{4}}^{\prime }+(\frac{1}{2}+a)b{L}^{\prime }]\mathrm{d}y\end{array}\}$$

由Theodorsen理论$(6.1.1)$式（注意$h=-w$）得

$$\begin{array}{c}{L}^{\prime }=2\pi {\rho }_{\infty }UbC(k)[U\theta -\frac{\partial w}{\partial t}+b(\frac{1}{2}-a)\frac{\partial \theta }{\partial t}]+\pi {\rho }_{\infty }{b}^2(U\frac{\partial \theta }{\partial t}-\frac{{\partial }^{2}w}{\partial t^2}-ba\frac{{\partial }^2\theta }{\partial t^2})\\ M_{\frac{1}{4}}^{\prime }=-\pi {\rho }_{\infty }{b}^3[U\frac{\partial \theta }{\partial t}-\frac{1}{2}\frac{{\partial }^{2}w}{\partial t^2}+b(\frac{1}{8}-\frac{a}{2})\frac{{\partial }^2\theta }{\partial t^2}]\end{array}\}$$

将上式和$(7.4)$式代入前面一个广义力表达式中，写成矩阵的形式：

$$\begin{gathered} \left[\begin{array}{c}{\mathbf{\Xi }}_{\mathbf{w}}\\ {\mathbf{\Xi }}_{\mathbf{\theta }}\end{array}\right]=-\pi {\rho }_{\infty }{b}^{2}l\left[\begin{array}{cc}\mathbf{E}& ba{\mathbf{A}}^{\mathrm{T}}\\ ba\mathbf{A}& b^2(a^2+\frac{1}{8})\mathbf{E}\end{array}\right]\left[\begin{array}{c}\ddot{\mathbf{\eta }}\\ \ddot{\mathbf{\varphi }}\end{array}\right]- \\ \pi {\rho }_{\infty }bUl\left[\begin{array}{cc}2C(k)\mathbf{E}& -b[1+2(\frac{1}{2}-a)C(k)]{\mathbf{A}}^{\mathrm{T}}\\ 2b(\frac{1}{2}+a)C(k)\mathbf{A}& b^2(\frac{1}{2}-a)[1-2(\frac{1}{2}+a)C(k)]\mathbf{E}\end{array}\right]\left[\begin{array}{c}\dot{\mathbf{\eta }}\\ \dot{\mathbf{\varphi }}\end{array}\right]- \\ \pi {\rho }_{\infty }b{U}^{2}l\left[\begin{array}{cc}0& -2C(k){\mathbf{A}}^{\mathrm{T}}\\ 0& -b(1+2a)C(k)\mathbf{E}\end{array}\right]\left[\begin{array}{c}\mathbf{\eta }\\ \mathbf{\varphi }\end{array}\right]\text{\hspace{1em}(7.7)} \end{gathered}$$

其中$\mathbf{A}$由$A_{ij}$组成，$E$表示单位矩阵，$0$表示零矩阵。

将$(7.5)$、$(7.6)$和$(7.7)$代入拉格朗日方程得到广义运动方程，写成矩阵的形式：

$$ml\left[\begin{array}{cc}\mathbf{E}& -b{x}_{\theta }{\mathbf{A}}^{\mathrm{T}}\\ -b{x}_{\theta }\mathbf{A}& b^2{r}^2\mathbf{E}\end{array}\right]\left[\begin{array}{c}\ddot{\mathbf{\eta }}\\ \ddot{\mathbf{\varphi }}\end{array}\right]+\left[\begin{array}{cc}\frac{EI}{l^3}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{B})& 0\\ 0& \frac{GJ}{l}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{T})\end{array}\right]\left[\begin{array}{c}\mathbf{\eta }\\ \mathbf{\varphi }\end{array}\right]=\left[\begin{array}{c}{\mathbf{\Xi }}_{\mathbf{w}}\\ {\mathbf{\Xi }}_{\mathbf{\theta }}\end{array}\right]\text{\hspace{1em}(7.8)}$$

其中对角矩阵$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{B})$和$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathbf{T})$的元素为

$$\begin{gathered} B_{ij}=({\alpha }_{i}l{)}^4 \\ {T}_{ij}=({\gamma }_{i}l{)}^2 \end{gathered}$$

由于Theodorsn理论的局限性，这个运动方程金对简谐运动有效。则可以用于第4节经典颤振分析的方法，即令：

$$\begin{array}{c}\dot{\eta }(t)=\overline{\eta }\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t)\\ \dot{\varphi }(t)=\overline{\varphi }\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{i}\omega t)\end{array}\}$$

其中$\omega$是简谐运动的频率。这样就可以导出颤振行列式，并求得颤振速度。

下面考虑$N_w=N_{\theta }=1$的情况，类似于第4节，写出一个类似与$(4.2.2)$式：

$$\begin{array}{c}\{\mu [1-(\frac{{\omega }_w}{\omega }{)}^2]+l_w\}\frac{{\overline{\eta }}_1}{b}+(-\mu x_{\theta }+l_{\theta })A_{11}{\overline{\varphi }}_1=0\\ (-\mu x_{\theta }+m_w)\frac{{\overline{\eta }}_1}{b}+\{\mu r^2[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+m_{\theta }\}{\overline{\varphi }}_1=0\end{array}\}$$

其中：

$$\begin{gathered} \mu =\frac{m}{\pi {\rho }_{\infty }{b}^2} \\ {l}_w=1-\frac{2\mathrm{i}C(k)}{k} \\ {l}_{\theta }=a+\frac{\mathrm{i}}{k}[1+2(\frac{1}{2}-a)C(k)]+\frac{2C(k)}{k^2} \\ {m}_w=a-\frac{2\mathrm{i}(\frac{1}{2}+a)C(k)}{k} \\ {m}_{\theta }=a^2+\frac{1}{8}-\frac{(\frac{1}{2}-a)[1-2(\frac{1}{2}+a)C(k)]\mathrm{i}}{k}+\frac{2(\frac{1}{2}+a)C(k)}{k^2} \\ {\omega }_w=({\alpha }_{1}l{)}^2\sqrt{\frac{\overline{EI}}{m{l}^4}} \\ {\omega }_{\theta }=\frac{\pi }{2}\sqrt{\frac{\overline{GJ}}{m{b}^2{r}^2{l}^2}} \\ {A}_{11}=0.958641 \end{gathered}$$

则可以写出颤振行列式：

$$\left|\begin{array}{cc}\mu [1-(\frac{{\omega }_w}{\omega }{)}^2]+l_w& (-\mu x_{\theta }+l_{\theta })A_{11}\\ (-\mu x_{\theta }+m_w)A_{11}& \mu r^2[1-(\frac{{\omega }_{\theta }}{\omega }{)}^2]+m_{\theta }\end{array}\right|=0$$

由此可以按照第4节的方法求得颤振边界。

8 颤振边界特性

由前面的章节介绍了以高度、速度和马赫数表示的颤振边界。给定一个飞行器，其颤振速度常与高度、来流声速有关。颤振速度通常用无量纲减缩速度$V=U/(b{\omega }_{\theta })$表示，称为减缩颤振速度$V_{\mathrm{F}}=U_{\mathrm{F}}/(b{\omega }_{\theta })$。高度通常由质量比$\mu =m/(\pi {\rho }_{\infty }{b}^2)$表示（来流密度${\rho }_{\infty }$是高度的函数）。来流声速通常用来流马赫数$M{a}_{\infty }=U/c_{\infty }$表示。
来流马赫数一定时，$V_{\mathrm{F}}$和$\mu$的关系如图所示（$\sigma =1/\sqrt{10},r=1/2,x_{\theta }=0,a=-3/10$）：

可以看出，颤振速度随着质量比的增加几乎线性增加。意味着飞行器在低空时要比高空时更容易颤振。

对于不同的飞行器，颤振边界也对结构参数非常敏感。

升力面单位展长的质量$m$影响着质量比$\mu =m/(\pi {\rho }_{\infty }{b}^2)$。下表为大气密度为从海平面到10000英尺之间飞行器的构型及其典型质量比数值：

飞行器类型	$\mu$
滑翔机和轻型飞机	5~15
普通航空飞机	10~15
商用运输机	15~30
战斗机	25~55
直升机叶片	65~110

可见滑翔机和轻型飞机更易发生颤振。

在质量比非常小的情况下，颤振速度随频率比$\sigma ={\omega }_h/{\omega }_{\theta }$的变换十分剧烈。如图所示无量纲颤振速度随频率比的变化曲线（$\mu =3,r=1/2,a=-1/5$）：

当$x_{\theta }=0.2$时，$\sigma =1.4$附近的颤振速度陡然降至极小，这一点具有最重要的实际意义。受其他参数影响，在某些频率比处，颤振速度会变得很低。这类凹坑可以大多数高性能飞机机翼的颤振速度随频率比的变化曲线看到，这类飞机具有相对较大的质量比$\mu$和正的静不平衡$x_{\theta }$。

弦向的无量纲距离参数（见第3节的模型）同样对颤振速度有很强的影响，如图所示无量纲颤振速度随$e$的变化曲线（$\mu =10,\sigma =1/\sqrt{2},r=1/2$）：

质心位置$e$和无量纲回转半径$r$必然同时改变，但是前者小百分比变化对颤振速度造成的影响要大于后者同样百分比变化对颤振速度造成的影响。如果质心相对于参考点前移，那么颤振速度通常会提高。

最后，注意到弯扭耦合型颤振频率介于${\omega }_h$和${\omega }_{\theta }$之间，即通常$\sigma <1$，不过也有可能出现颤振频率高于${\omega }_{\theta }$的情况。

值得注意的是，目前的简化理论表明，在弦向偏移参数$e$和$a$的某些组合下不会发生颤振。当$\sigma$取小值时，若$e\le -1/2$则不会发生颤振；当$\sigma$取大值时，若$a\ge -1/2$则不会发生颤振。当质心、弹性轴和气动中心三者全部重合（即$e=a=-1/2$）时，对于任意的$\sigma$和$r$组合似乎都不会发生颤振。

颤振边界一种可能的表示方式是与飞行包线一起给出。一架马赫数为2的攻击机的典型飞行包线如图所示：

图中记为No.1和No.2的曲线表示颤振边界。可以看出此飞行器容易发生颤振的条件：低高度，跨声速流，大动压。

9 总结

气动弹性颤振的分析方法有：p法、经典颤振分析方法、k法和p-k法。

---

本文章使用limfx的vsocde插件快速发布