考虑热退化和阻尼效应的壁板颤振特性分析

Analysis of flutter characteristics of panel considering thermal degradation and damping effects

摘要

为系统揭示热退化与阻尼效应对壁板颤振边界与响应的作用机制，构建含热退化与阻尼项的壁板气动弹性模型，并基于稳定性判据开展了关键参数的灵敏度分析。首先，基于Von Kármán板大变形理论、Kelvin阻尼模型及一阶活塞理论建立考虑热退化和阻尼效应的二维超声速壁板动力学方程，并通过Galerkin法进行空间离散。然后，基于李雅普诺夫间接法和Routh-Hurwitz准则，得到不同热退化程度下的稳定性区域图，并确定各区域内平衡点的数量及其稳定性。最后，使用四阶Runge-Kutta方法求解非线性常微分方程组，得到壁板非线性气动弹性响应，并通过时间历程图、相轨迹图和分岔图等非线性动力学描述工具进行分析。研究结果表明，热退化和材料阻尼的存在会显著缩小壁板的稳定区域范围且使得壁板厚度成为影响稳定性区域的因素之一。热退化系数的增大不仅会导致壁板振动幅值的增大，还会使得分岔点提前出现，并加速分岔演化过程。并且，热退化效应还会显著提高壁板响应类型的多样性和对系统参数的敏感性。此外，材料阻尼效应能够有效抑制壁板的类周期振动与混沌振动响应。

Abstract

To systematically investigate the mechanisms of thermal degradation and damping effects on panel flutter boundaries and responses， an aeroelastic model incorporating thermal degradation and damping terms was developed， followed by a sensitivity analysis of key parameters based on stability criteria. Firstly， based on the von Kármán plate large deformation theory， the Kelvin damping model， and the first-order piston theory， a two-dimensional supersonic panel dynamic equation was established. This equation considered thermal degradation and damping effects. Spatial discretization was achieved using the Galerkin method. Next， the Lyapunov indirect method and the Routh-Hurwitz criterion were used to obtain stability region diagrams for various degrees of thermal degradation. This identified the number and stability of equilibrium points in each region. Finally， the nonlinear ordinary differential equations were solved using the fourth-order Runge-Kutta method to determine the panel’s nonlinear aeroelastic response. This response was analyzed using nonlinear dynamic tools， such as time history plots， phase trajectory diagrams， and bifurcation diagrams. The results indicate that both thermal degradation and material damping significantly reduce the stability region of the panel. Panel thickness also influences the stability region. An increase in the thermal degradation coefficient not only amplifies the panel’s vibration amplitude but also causes the bifurcation points to appear earlier， accelerating the bifurcation evolution process. Furthermore， thermal degradation significantly increases the diversity of panel response types and their sensitivity to system parameters. In addition， material damping effectively suppresses the panel’s quasi-periodic and chaotic vibrations responses.

Graphical abstract

关键词

热退化 / 材料阻尼 / 壁板 / 稳定性边界 / 气动弹性响应 / 颤振

Key words

thermal degradation / material damping / panel / stability boundary / aeroelastic response / flutter

随着航空航天技术的不断发展，新型材料（如复合材料和高温合金）因其优异的强度、质量比和耐高温性能，在超音速飞机壁板制造中的应用比例迅速增大。然而，超音速飞行环境中复杂的热-力耦合条件使得材料的热退化效应和材料阻尼特性成为影响飞机性能和安全性的关键因素。热退化不仅会导致材料力学性能的衰减，还可能引发结构的稳定性问题；而材料阻尼作为抑制振动和能量耗散的重要机制，直接关系到飞行器的动态性能和疲劳寿命。因此，对新型材料在高温、高速环境下的热退化行为及其阻尼效应进行深入研究，具有重要的理论意义和工程价值。

近几十年来，学者对壁板颤振进行了广泛研究。Dowell^［1-2］在20世纪60年代率先研究了壁板的热颤振，认为几何非线性对热颤振影响显著，并将壁板颤振分析模型分为四类。在此基础上，Gary等^［3］提出了高阶非线性活塞理论，发展出第五类分析模型。随着计算气动力方法精确度的提高，Cheng等^［4］结合前人的结论提出第六类分析模型。夏巍等^［5］和杨智春等^［6］建立了超声速二维受热壁板的气动弹性模型，并分析了其稳定性。Ye等^［7］研究了激波振荡对壁板非线性动力学的影响，发现激波振荡的幅值和频率显著影响极限环振幅及混沌行为。

壁板在超音速气流中不可避免地受到气动加热效应的显著影响^［8］。研究显示，热效应通过两种方式影响壁板气动弹性响应：1）材料因热膨胀产生面内热应力；2）材料受热导致特征属性的改变，即热退化现象。当前文献多集中于第一种情况对壁板颤振的影响。Yuan等^［9］分析了加筋复合材料板在超声速气流和热效应下的颤振特性。Xie等^［10］提出了一种基于正交分解（proper orthogonal decomposition，POD）方法的非线性壁板颤振分析方法。Chai等^［11-12］分析了复合层合板在超音速气流中热-气动弹性特性及其主动颤振控制。Miller等^［13］对超音速气流中层合板的气动弹性响应进行了研究。

实际上，除热效应外，壁板的材料阻尼也显著影响颤振边界和气动弹性响应特性。Pourtakdoust等^［14］提出材料阻尼对壁板稳定性边界有重要影响。Wang等^［15］建立了同时受到上下表面气动载荷作用的加热黏弹性壁板模型，发现高温升下黏弹性阻尼对壁板稳定性和后颤振混沌运动起到稳定作用。

综上，本文建立了一个系统的理论模型，研究超声速流中的壁板在同时考虑热效应的两种影响方式（热应力与热退化）及材料阻尼下的稳定性边界与气动弹性响应特性。

1 计算模型

1.1 运动微分方程

超音速流中二维简支壁板模型如图1所示。其沿x方向长度为a，沿y方向长度为b\to \mathrm{\infty }，厚度为h，壁板上表面作用沿x坐标轴正方向的超音速流，其速度为U_{\mathrm{\infty }}，密度为{\rho }_{\mathrm{\infty }}，马赫数为Ma， w 为壁板沿x方向的横向位移。

基于Von Kármán板大变形理论和Kirchhoff板理论，结合简支壁板的边界条件，得到二维壁板的运动微分方程为

D\frac{{\partial }^{\mathrm{4}}w}{\partial x^{\mathrm{4}}}-\left(N_x-N_x^{\left(T\right)}\right)\frac{{\partial }^{\mathrm{2}}w}{\partial x^{\mathrm{2}}}+{\rho }_{m}h\frac{{\partial }^{\mathrm{2}}w}{\partial t^{\mathrm{2}}}-q_a=\mathrm{0}（1）

式中：D=Eh/\left[\mathrm{12}\left(\mathrm{1}-{\mu }^{\mathrm{2}}\right)\right]为壁板弯曲刚度；N_x=Eh/\left[\mathrm{2}a\left(\mathrm{1}-{\mu }^{\mathrm{2}}\right)\right]{\int }_{\mathrm{0}}^a{\left(\partial w/\partial x\right)}^{\mathrm{2}}\mathrm{d}x为壁板几何非线性力；N_x^{\left(T\right)}={\int }_{-h/\mathrm{2}}^{h/\mathrm{2}}E{\alpha }_{r}T/\left(\mathrm{1}-\mu \right)\mathrm{d}z为面内热应力载荷；{\rho }_m为壁板的密度；q_a为气动载荷；E为弹性模量；{\alpha }_r为热膨胀系数；\mu为泊松比；T为壁板的温升。

1.2 气动载荷

根据修正后的一阶活塞理论，可得气动力为

\begin{array}{l}{q}_a=P_{\mathrm{\infty }}-P_p=\\ -\frac{\mathrm{2}{q}_{\mathrm{\infty }}}{\beta }\left(\frac{\partial w}{\partial x}+\frac{\mathrm{1}}{U_{\mathrm{\infty }}}·\frac{M{a}^{\mathrm{2}}-\mathrm{2}}{M{a}^{\mathrm{2}}-\mathrm{1}}·\frac{\partial w}{\partial t}\right)\end{array}（2）

式中：\beta =\sqrt[]{M{a}^{\mathrm{2}}-\mathrm{1}}，为修正因子。

1.3 引入热退化及材料阻尼并无量纲化

温度对于壁板气动弹性的影响不仅体现在热应力方面，还会对壁板的材料属性产生影响，较为突出的就是弹性模量与热膨胀系数。通过查阅材料手册，可知弹性模量和热膨胀系数大多与温度呈线性关系，因此有

E=E_{\mathrm{0}}\left(\mathrm{1}+eT\right)（3）

{\alpha }_r={\alpha }_{\mathrm{0}}\left(\mathrm{1}+\alpha T\right)（4）

式中：E_{\mathrm{0}}为升温前材料的弹性模量；e为弹性模量的线性增量系数；{\alpha }_{\mathrm{0}}为升温前材料的热膨胀系数；\alpha为热膨胀系数的线性增量系数。

采用应用较为广泛的Kelvin阻尼模型，将E_{\mathrm{0}}用E_{\mathrm{0}}\left(\mathrm{1}+g\cdot \partial /\partial t\right)代换，则式（3）可写为

E=E_{\mathrm{0}}\left(\mathrm{1}+eT\right)\left(\mathrm{1}+g\frac{\partial }{\partial t}\right)（5）

将式（2）—（5）代入式（1），并采取合适的无量纲化形式，可得考虑热退化于材料阻尼效应的二维超声速壁板运动微分方程。

\begin{array}{l}-\mathrm{6}{C}_e\left(\mathrm{1}+\zeta \frac{\partial }{\partial \tau }\right)\left[{\int }_{\mathrm{0}}^{\mathrm{1}}{\left(\frac{\partial W}{\partial \xi }\right)}^{\mathrm{2}}\mathrm{d}\xi \right]\frac{{\partial }^{\mathrm{2}}W}{\partial {\xi }^{\mathrm{2}}}+\\ C_e{C}_{\alpha }{R}_x\left(\mathrm{1}+\zeta \frac{\partial }{\partial \tau }\right)\frac{{\partial }^{\mathrm{2}}W}{\partial {\xi }^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}W}{\partial {\tau }^{\mathrm{2}}}+\lambda \frac{\partial W}{\partial \xi }+\\ C_e\frac{{\partial }^{\mathrm{4}}W}{\partial {\xi }^{\mathrm{4}}}+\sqrt[]{\lambda R_m}\frac{\partial W}{\partial \tau }=\mathrm{0}\end{array}（6）

式中：\xi =\frac{x}{a}；W=\frac{w}{h}；\tau =t\sqrt[]{\frac{D_{\mathrm{0}}}{{\rho }_{m}h{a}^{\mathrm{4}}}}；\lambda =\frac{\mathrm{2}{q}_{\mathrm{\infty }}{a}^{\mathrm{3}}}{\beta D_{\mathrm{0}}}；\nu =\frac{{\rho }_{\mathrm{\infty }}a}{{\rho }_{m}h}；R_x=T/\overline{T}；\overline{T}=\frac{h^{\mathrm{2}}}{\mathrm{12}{\alpha }_{\mathrm{0}}\left(\mathrm{1}+\mu \right)a^{\mathrm{2}}}；D_{\mathrm{0}}=\frac{E_{\mathrm{0}}{h}^{\mathrm{3}}}{\mathrm{12}\left(\mathrm{1}-{\mu }^{\mathrm{2}}\right)}；R_m=\left(\frac{M{a}^{\mathrm{2}}-\mathrm{2}}{M{a}^{\mathrm{2}}-\mathrm{1}}\right)\frac{\nu }{\beta }；\zeta =g\sqrt[]{\frac{D_{\mathrm{0}}}{{\rho }_{m}h{a}^{\mathrm{4}}}}；C_e=\mathrm{1}+eT；C_{\alpha }=\mathrm{1}+\alpha T。

1.4 Galerkin离散

为满足壁板的简支边界条件，设位移函数为

W\left(\xi ,\tau \right)={\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{\infty }}}{a}_m\left(\tau \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\mathrm{\pi }\xi \right)（7）

将式（7）代入式（6）后，在等式两侧分别乘以\mathrm{s}\mathrm{i}\mathrm{n}\left(r\mathrm{\pi }\xi \right)，r=\mathrm{1,2},\cdots ,\mathrm{\infty }，并沿无量纲化板长积分，可得到二阶非线性常微分方程组。壁板的变形主要以前N阶模态为主，令{\dot{a}}_r=a_{N+r}，代入二阶N维非线性常微分方程组，可到一阶\mathrm{2}N维非线性常微分方程组为

\begin{array}{l}{\dot{a}}_r=a_{N+r}\\ {\dot{a}}_{N+r}=-\left\{\begin{array}{l}{C}_e{(r\mathrm{\pi })}^{\mathrm{4}}\\ +\mathrm{3}{C}_e\left[{\displaystyle \sum _{m=\mathrm{1}}^N}{a}_m^{\mathrm{2}}{(m\mathrm{\pi })}^{\mathrm{2}}\right]{(r\mathrm{\pi })}^{\mathrm{2}}\\ -C_e{C}_{\alpha }{R}_x{(r\mathrm{\pi })}^{\mathrm{2}}\end{array}\right\}{a}_r-\end{array}

\lambda {\displaystyle \sum _{\begin{array}{l}m=\mathrm{1}\\ m\ne r\end{array}}^N}\frac{\mathrm{2}mr}{r^{\mathrm{2}}-m^{\mathrm{2}}}\left[\mathrm{1}-{(-\mathrm{1})}^{r+m}\right]a_m-

\begin{array}{l}\sqrt[]{\lambda R_m}{a}_{N+r}-\\ \zeta C_e\left(\begin{array}{l}{(r\mathrm{\pi })}^{\mathrm{4}}-C_{\alpha }{R}_x{(r\mathrm{\pi })}^{\mathrm{2}}\\ +\mathrm{3}\left[{\displaystyle \sum _{m=\mathrm{1}}^N}{a}_m^{\mathrm{2}}{(m\mathrm{\pi })}^{\mathrm{2}}\right]{(r\mathrm{\pi })}^{\mathrm{2}}\end{array}\right)a_{N+r},\\ r=\mathrm{1,2},\cdots ,N\end{array}（8）

2 壁板稳定性边界

由于材料属性基本已经被包含在无量纲化参数中，改变材料导致材料属性的变化基本不会对分析过程产生影响，只是不同材料对应着不同的热退化系数及材料阻尼系数。因此，假设壁板尺寸、属性及所处流体属性分别为a=\mathrm{1} m，h=\mathrm{0.005} m，E_{\mathrm{0}}=\mathrm{1.1}\times {\mathrm{10}}^{\mathrm{11}} Pa，{\alpha }_{\mathrm{0}}=\mathrm{8.71}\times {\mathrm{10}}^{-\mathrm{6}}/℃，\mu =\mathrm{0.3}，{\rho }_m=\mathrm{4}\mathrm{429} kg/m³，Ma=\mathrm{3}，{\rho }_{\mathrm{\infty }}=\mathrm{0.364} kg/m³。

对于式（8），暂取N=\mathrm{2}，可得到一阶四维非线性常微分方程。当{\dot{a}}_k=\mathrm{0},k=\mathrm{1,2},\cdots ,\mathrm{2}N时，可以得到该二维壁板热气弹问题的定态方程为

\left\{\begin{array}{l}{C}_e{\mathrm{\pi }}^{\mathrm{4}}\left[\mathrm{1}+\mathrm{3}\left(a_{\mathrm{1}}^{\mathrm{2}}+\mathrm{4}{a}_{\mathrm{2}}^{\mathrm{2}}\right)-C_{\alpha }{R}_x/{\mathrm{\pi }}^{\mathrm{2}}\right]a_{\mathrm{1}}=\frac{\mathrm{8}}{\mathrm{3}}\lambda a_{\mathrm{2}}\\ \mathrm{4}{C}_e{\mathrm{\pi }}^{\mathrm{4}}\left\{\mathrm{4}+\mathrm{3}\left[a_{\mathrm{1}}^{\mathrm{2}}+\mathrm{4}{a}_{\mathrm{2}}^{\mathrm{2}}\right]-C_{\alpha }{R}_x/{\mathrm{\pi }}^{\mathrm{2}}\right\}{a}_{\mathrm{2}}=-\frac{\mathrm{8}}{\mathrm{3}}\lambda a_{\mathrm{1}}\end{array}\right.（9）

显然，a_{\mathrm{1}}=a_{\mathrm{2}}=\mathrm{0}是该定态方程的一个解，这个解代表平板的初始静力平衡状态。定态方程的解对应的状态点称为不动点，在该不动点（记为p_{\mathrm{0}}）附近的稳定性通过李雅普诺夫第二方法来进行判断。在该不动点邻域内的伴随线性方程为\dot{\mathit{\epsilon }}=\mathit{J}\mathit{\epsilon }，\mathit{J}为该不动点处的Jacobi矩阵，可以得到伴随线性系统的特征方程为\left|\mathit{J}-{\lambda }^J\mathit{I}\right|=\mathrm{0}。通过分析伴随线性系统的稳定性，进而推断原非线性系统的稳定性。

在不动点p_{\mathrm{0}}处的Jacobi矩阵为

\mathit{J}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ j_{\mathrm{31}}& j_{\mathrm{32}}& j_{\mathrm{33}}& j_{\mathrm{34}}\\ j_{\mathrm{41}}& j_{\mathrm{42}}& j_{\mathrm{43}}& j_{\mathrm{44}}\end{array}\right]（10）

伴随系统特征方程为

\left|\mathit{J}-{\lambda }^J\mathit{I}\right|=({\lambda }^J{)}^{\mathrm{4}}+j_{\mathrm{3}}({\lambda }^J{)}^{\mathrm{3}}+j_{\mathrm{2}}({\lambda }^J{)}^{\mathrm{2}}+j_{\mathrm{1}}{\lambda }^J+j_{\mathrm{0}}=\mathrm{0}（11）

根据特征方程，列出Routh表，基于Routh-Hurwitz稳定判据，该特征方程所有特征值的实部为负，即系统稳定的条件为

\begin{array}{l}{j}_{\mathrm{0}}>\mathrm{0},j_{\mathrm{1}}>\mathrm{0},j_{\mathrm{3}}>\mathrm{0}\\ j_{\mathrm{3}}{j}_{\mathrm{2}}{j}_{\mathrm{1}}-j_{\mathrm{1}}^{\mathrm{2}}-j_{\mathrm{3}}^{\mathrm{2}}{j}_{\mathrm{0}}>\mathrm{0}\end{array}（12）

令\mathrm{\Delta }T=R_x/{\mathrm{\pi }}^{\mathrm{2}}=T/\left(\overline{T}{\mathrm{\pi }}^{\mathrm{2}}\right)，记为无量纲化温度增量，无量纲动压\lambda >\mathrm{0}，得到该系统在不动点p_{\mathrm{0}}稳定时，\lambda和\mathrm{\Delta }T满足关系

（13）

（14）

{\lambda }^{\mathrm{2}}>\frac{\mathrm{9}}{\mathrm{16}}{\mathrm{\pi }}^{\mathrm{8}}{C}_e^{\mathrm{2}}\left(\mathrm{1}-\mathrm{\Delta }T{C}_{\alpha }\right)\left(\mathrm{\Delta }T{C}_{\alpha }-\mathrm{4}\right)（15）

（16）

据上述分析和不等式可知，定态方程的零解（不动点p_{\mathrm{0}}）一定存在。在满足式（13）—（16）所述不等式边界，p_{\mathrm{0}}为稳定不动点。当超过该边界时，系统的不动点p_{\mathrm{0}}将变为不稳定不动点。

下面考虑定态方程的非零解，即壁板发生热屈曲变形的情况。将式（9）写为矩阵方程形式。

\mathit{A}\mathit{a}=\mathrm{0}（17）

令M=\mathrm{3}\left(a_{\mathrm{1}}^{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}+\mathrm{4}{a}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}\right)，根据方程有非零解条件：\left|\mathit{A}\right|=\mathrm{0}，可求解得到

\begin{array}{l}M=-\frac{\mathrm{5}}{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}+R_x{C}_{\alpha }\pm \\ \frac{\mathrm{1}}{\mathrm{6}{C}_e{\mathrm{\pi }}^{\mathrm{2}}}\sqrt[]{\mathrm{81}{C}_e^{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{8}}-\mathrm{64}{\lambda }^{\mathrm{2}}}\end{array}（18）

若要使M为实数，则需

\lambda \le \frac{\mathrm{9}}{\mathrm{8}}{C}_e{\mathrm{\pi }}^{\mathrm{4}}（19）

若定态方程存在非零解，即壁板发生热屈曲变形，需要满足式（19）所述不等式边界，当超过该边界时，系统将不存在屈曲不动点。

根据\lambda的范围，令

\lambda =\frac{\mathrm{9}}{\mathrm{8}}{C}_e{\mathrm{\pi }}^{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\theta \right),\mathrm{0}\le \theta \le \frac{\mathrm{\pi }}{\mathrm{4}}（20）

将式（20）代入式（18）可得

M=-\frac{\mathrm{5}}{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}+R_x{C}_{\alpha }\pm \frac{\mathrm{3}}{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\theta \right)（21）

当M=-{\mathrm{\pi }}^{\mathrm{2}}+R_x{C}_{\alpha }-\mathrm{3}{\mathrm{\pi }}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta时，将其与式（20）代入式（17），整理后可得

\mathit{A}\mathit{a}=\left[\begin{array}{cc}-\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta & -\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\theta )\\ \mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\theta )& \mathrm{4}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta \end{array}\right]\left[\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]=\mathrm{0}（22）

求解该定态方程可得一组解为

\left[\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]=\pm \sqrt[]{\frac{M}{\mathrm{12}{\mathrm{\pi }}^{\mathrm{2}}}}\left[\begin{array}{c}\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ -\mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right]（23）

将这组解对应的不动点分别记为p_{\mathrm{1}}和p_{\mathrm{2}}，可以发现这组不动点关于中性面对称。由于M>\mathrm{0}，即-{\mathrm{\pi }}^{\mathrm{2}}+R_x{C}_{\alpha }-\mathrm{3}{\mathrm{\pi }}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta >\mathrm{0}，并根据式（20），可得\mathrm{\Delta }T-\lambda的不等式。

\mathrm{\Delta }T>\frac{\mathrm{1}}{\mathrm{2}{C}_{\alpha }}\left[\mathrm{5}-\mathrm{3}\sqrt[]{\mathrm{1}-{\left(\frac{\mathrm{8}\lambda }{\mathrm{9}{C}_e{\mathrm{\pi }}^{\mathrm{4}}}\right)}^{\mathrm{2}}}\right]（24）

根据上述不等式可知，若系统存在不动点p_{\mathrm{1}}、p_{\mathrm{2}}，需要满足式（19）和（24）所述不等式边界。

当M=-\mathrm{4}{\mathrm{\pi }}^{\mathrm{2}}+R_x{C}_{\alpha }+\mathrm{3}{\mathrm{\pi }}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta时，同理可得

\left[\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\end{array}\right]=\pm \sqrt[]{\frac{M}{\mathrm{12}{\mathrm{\pi }}^{\mathrm{2}}}}\left[\begin{array}{c}\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ -\mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right]（25）

将这组解对应的不动点记为p_{\mathrm{3}}和p_{\mathrm{4}}，同样可以发现，这组不动点同样关于中性面对称。同理由M>\mathrm{0}可得

\mathrm{\Delta }T>\frac{\mathrm{1}}{\mathrm{2}{C}_{\alpha }}\left[\mathrm{5}+\mathrm{3}\sqrt[]{\mathrm{1}-{\left(\frac{\mathrm{8}\lambda }{\mathrm{9}{C}_e{\mathrm{\pi }}^{\mathrm{4}}}\right)}^{\mathrm{2}}}\right]（26）

同理，若系统存在不动点p_{\mathrm{3}}、p_{\mathrm{4}}，需要满足式（19）和（26）所述不等式边界。

将不动点p_{\mathrm{1}}、p_{\mathrm{2}}和p_{\mathrm{3}}、p_{\mathrm{4}}对应的解代入对应的Jacobi矩阵并求解矩阵特征值可知，p_{\mathrm{1}}、p_{\mathrm{2}}始终稳定；p_{\mathrm{3}}、p_{\mathrm{4}}始终不稳定。

暂不考虑热退化及材料阻尼效应时，根据式（13）—（16）、（19）、（24）和（26）所述不等式边界，绘制\mathrm{\Delta }T-\lambda的不动点稳定性区域图，并根据文献中给出的表达式绘制稳定性区域图^［5］，与本文结果进行对比，如图2所示。

图2中：AG称为颤振失稳边界，BFGC称为屈曲失稳边界，FD称为屈曲-颤振边界。在这种条件下，颤振失稳边界和屈曲-颤振边界均与屈曲失稳边界相切。稳定性边界将图2划分为5个区域。稳定性区域不动点及其稳定性情况如表1所示。

表1中：Ⅰ区域称为基本稳定区；Ⅱ、Ⅲ区域称为热屈曲区；Ⅳ区域称为基本颤振区；Ⅴ区域称为过渡区。壁板在两种情况下会发生动态失稳：

1）当区域内不存在稳定不动点，则壁板必然发生颤振，如区域Ⅳ。本文主要讨论该区域发生的动态失稳现象，即发生颤振。

2）当系统既存在渐近稳定的不动点又存在不稳定的不动点时，则壁板可能静态稳定在某一个稳定不动点，也有可能动态失稳而发生颤振，如区域Ⅱ、Ⅲ、Ⅴ，这种颤振称为二次失稳型颤振。

通过表2中不动点对应特征值实部的正负可以判断不动点在该区域的稳定性，当4个特征值实部全为负值时，不动点为稳定不动点；当特征值实部出现正值时，不动点为不稳定不动点。可以看出表2与表1中的结论相互验证，确保了前文理论分析的正确性。

考虑不同方面热退化效应下的稳定性边界，如图3所示。通过图3可以发现，只考虑热膨胀系数的热退化效应时，边界FD不受影响；只考虑弹性模量的热退化效应时，3个边界都受到影响；而完整考虑热退化效应时，边界的变化情况结合了前两者的特点。除基本颤振区域外，其余区域均向坐标系左下角缩减，这意味着考虑热退化效应后，壁板更容易出现动态失稳状态。

在上述基础上，绘制不同无量纲材料阻尼系数下的稳定性边界图，如图4所示。可以发现，随材料阻尼增大，颤振失稳边界向坐标系左下方移动，并且与屈曲失稳边界的相切关系不复存在；过渡区域Ⅴ逐渐减小，最终消失。

不同壁板厚度下，考虑热退化与材料阻尼效应壁板的稳定性区域如图5所示。对比图5中不同厚度下的稳定性区域图可以发现，壁板厚度增大时，热退化和材料阻尼效应对壁板稳定性边界的影响愈加显著。这是由于随壁板厚度的增加，受热退化影响的材料体积也在增加，热退化对壁板的影响进一步加强。

3 壁板气动弹性响应特性分析

使用RK4数值积分方法求解非线性微分方程并得到壁板位移响应。截至模态数大于4时即可较为准确地描述壁板实际响应^［6］，本文选择N=\mathrm{6}。为更清晰地看出热退化及材料阻尼的影响，选择壁板厚度h=\mathrm{0.01}。选择特征点\xi =\mathrm{0.75}进行位移响应观测。得到不同热退化程度与材料阻尼下系统分叉参数λ的分岔图，如图6所示。

通过对比图6a—6c可知，随热退化程度的加深，在较低动压下（\lambda \le \mathrm{200}）分岔形式逐渐复杂。从图6a—6c中的分岔图中可以发现，有分叉点的出现，根据分岔图的定义可知出现了类周期振动。在\lambda =\mathrm{409,605}时出现类周期振动，并且随热退化的加深，类周期振动出现时对应的动压降低、幅值增大且范围扩大。而对比图6d—6i可知，随材料阻尼的增大，壁板高动压下的混沌振动被极大程度地抑制，并且前文提到的类周期振动也被抑制。但是在较低动压时出现了周期倍化现象，并提高了振动的幅值。

4 结论

研究了考虑热退化的超声速阻尼壁板稳定性边界及非线性气动弹性响应特性。基于李雅普诺夫间接法和Routh-Hurwitz准则，推导了稳定性边界的理论解，并得到稳定性边界图。通过数值积分方法求解微分方程，结合时间历程图、相轨迹图及分岔图，分析了热退化后的阻尼壁板的气动弹性响应特性。得到以下结论：1）随热退化效应的加深，静稳定及屈曲区域范围均有所缩小，而颤振区域范围显著扩大。2）随阻尼系数的增大，壁板的静稳定区域范围显著缩小，颤振区域范围显著扩大，而屈曲区域几乎不受影响。3）壁板厚度增大时，热退化和材料阻尼效应对壁板稳定性边界的影响愈加显著。这说明热退化和材料阻尼效应对壁板的影响机制与壁板厚度密切相关。4）热退化系数的增大不仅会导致壁板振动幅值的增大，还会使分岔点提前出现，并加速分岔演化过程。5）材料阻尼可以有效地抑制壁板类周期振动和混沌振动响应。

参考文献

[1] Dowell E H.Nonlinear oscillations of a fluttering plate［J］.AIAA Journal，1966，4（7）：1267-1275.

[2] Dowell E H.Panel flutter：a review of the aeroelastic stability of plates and shells［J］.AIAA Journal，1970，8（3）：385-399.

[3] Gray C， Mei C.Large-amplitude finite element flu-tter analysis of composite panels in hypersonic flow［J］.AIAA Journal，1993，31（6）：1090-1099.

[4] Cheng G F， Mei C.Finite element modal formulation for hypersonic panel flutter analysis with thermal effects［J］.AIAA Journal，2004，42（4）：687-695.

[5] 夏巍， 杨智春. 超音速气流中受热壁板的稳定性分析［J］. 力学学报， 2007（5）： 602-609.

[6] 杨智春， 夏巍. 超音速气流中二维壁板的非线性热颤振响应分析［J］. 振动工程学报， 2009， 22（3）： 221-226.

[7] Ye L Q， Ye Z Y， Ye K，et al.Aeroelastic stability and nonlinear flutter analysis of viscoelastic heated panel in shock-dominated flows［J］.Aerospace Scien- ce and Technology，2021，117：106909.

[8] Ye K， Ye Z Y， Li C N，et al.Effects of the aerothermoelastic deformation on the performance of the three-dimensional hypersonic inlet［J］.Aerospace Science and Technology，2019，84：747-762.

[9] Yuan K H， Qiu Z P.Nonlinear flutter analysis of stiffened composite panels in supersonic flow［J］.Science China Physics，Mechanics and Astro-nomy，2010，53（2）：336-344.

[10] Xie D， Xu M， Dai H H，et al.Proper orthogonal decomposition method for analysis of nonlinear panel flutter with thermal effects in supersonic flow［J］.Journal of Sound and Vibration，2015，337：263-283.

[11] Chai Y Y， Song Z G， Li F M.Active aerother-moelastic flutter suppression of composite laminated panels with time-dependent boundaries［J］.Composite Structures，2017，179：61-76.

[12] Chai Y Y， Li F M， Song Z G，et al.Aerothermoelastic flutter analysis and active vibration suppression of nonlinear composite laminated panels with time-dependent boundary conditions in supersonic airflow［J］.Journal of Intelligent Material Systems and Structures，2018，29（4）：653-668.

[13] Miller B A， Mcnamara J J.Impact of aerothermodynamic model fidelity on aerothermoelastic response of a skin panel［J］.AIAA Journal，2018，56（12）：5028-5032.

[14] Pourtakdoust S H， Fazelzadeh S A.Chaotic analysis of nonlinear viscoelastic panel flutter in supersonic flow［J］.Nonlinear Dynamics，2003，32（4）：387-404.

[15] Wang X C， Yang Z C， Wang W，et al.Nonlinear viscoelastic heated panel flutter with aerodynamic loading exerted on both surfaces［J］.Journal of Sound and Vibration，2017，409：306-317.

基金资助

辽宁省教育厅基本科研面上项目(JYTMS20230253)

辽宁省教育厅基本科研面上项目(JYTMS20230224)

辽宁省科技厅联合计划面上项目(2025-MSLH-565)