超音速流中纤维增强FGM梁的气动弹性动力学特性研究

Study on Aeroelastic Dynamic Characteristics of Fiber-Reinforced FGM Beams in Supersonic Flow

摘要: 本文基于经典梁理论(CBT)，研究了在高超音速作用下，活塞气动力理论的非线性效应对纤维增强功能梯度材料(FGM)的颤振特性的影响。首先通过混合率模型来表征纤维增强FGM梁的材料属性，然后通过Hamilton原理推导出只考虑横向振动的纤维增强FGM梁的非线性气动弹性偏微分方程，利用Galerkin方法，把该方程转化为非线性常微分方程，再利用Hurwitz行列式，把该方程的求根问题用以判定Hopf分叉，得到不同温度应力下梁的无量纲临界流速和无量纲临界频率。最后通过Runge-Kutta法得到纤维体积分数和无量纲温升对无量纲临界动压的影响。

Abstract: Based on the classical beam theory (CBT), the influence of the nonlinear effects of piston aerodynamic theory on the flutter characteristics of fiber-reinforced functionally graded materials (FGM) at hypersonic speeds is studied in this paper. First, the material properties of fiber-reinforced FGM beams are characterized by the mixing rate model. Then, the nonlinear aeroelastic partial differen-tial equation of fiber-reinforced FGM beams that only consider the transverse vibration is derived from Hamilton’s principle. The Galerkin method is used to convert the equation into a nonlinear ordinary differential equation, and the Hurwitz determinant is used to determine the Hopf bifurcation, the dimensionless critical velocity and the dimensionless critical frequency of the beam under different temperature stresses are obtained. Finally, fiber volume fraction and dimensionless tem-perature rise on dimensionless critical dynamic pressure are obtained by the Runge-Kutta method.

1. 引言

功能梯度材料(FGM)是两种或者两种以上材料按一定规律组合起来的材料组分按梯度变化的先进材料，因其优异的性能，越来越多的学者对其进行了研究 [1] 。Lai等 [2] 、蹇越傲等 [3] 研究了在了不同载荷下功能梯度梁的大振幅振动问题。滕兆春等 [4] 基于Timoshenko梁变形理论研究了弹性地基上转动功能梯度梁的自由振动问题。现在，增强型复合材料因其更好的刚度和强度，常被应用于工程实际，为此增强功能梯度结构也有了更多的研究。Liew等 [5] 研究了石墨烯薄片(GPL)增强功能梯度多孔板的静态线性弹性、固有频率和屈曲行为。柏冬军等 [6] 研究了将压电功能梯度层合梁中的两节点力–电–热耦合梁单元用于功能梯度层合梁的振动优化。功能梯度材料也因其能很好地减缓热应力集中，在复杂的巡航环境中保证结构的强度，为此在飞行器上有诸多的应用。飞行器结构主要面临气动力和热载荷共同作用下的恶劣环境，现在对于超音速环境下的研究也有很多 [7] 。Mehri等 [8] [9] 研究了碳纳米管增强功能梯度材料板壳的气动弹性问题。Zhang等 [10] 研究了一种由功能梯度梁和压电材料组成的新型结构在超音速中的能量收集问题。Song等 [11] 提出了一种轴向功能梯度材料并研究了其对颤振的抑制效果。Mahmoudkhani等 [12] 对简支功能梯度壳的气动弹性问题做了研究。Fazelzadeh等 [13] 、Sohn等 [14] 研究了各参数对于超音速功能梯度板的散热效果的影响。现在对于各种增强功能梯度材料的研究有很多，对于超音速流中板壳的研究也有很多。但对于超音速流中纤维增强功能梯度梁的研究却很少，所以研究其气动弹性问题具有重要的现实意义。

目前，关于纤维增强FGM梁的研究国内外鲜有文献研究，本文基于经典梁理论，考虑将FGM梁作为基体进行均匀分布的纤维增强，由指数模型和混合率公式计算纤维增强FGM梁的等效物性参数。在超音速和热效应的前提下，通过Hamilton原理得出气动弹性控制微分方程。利用Galerkin方法，转化该方程为非线性常微分方程，再利用Hurwitz行列式，把该方程的求根问题用以判定Hopf分叉，得到不同温度应力下梁的无量纲临界流速和无量纲临界频率。最后通过Runge-Kutta法得到纤维体积分数和无量纲温升对无量纲临界动压的影响。

2. 基本方程

建立一个梁长为L，宽为b，高为h，弹性模量为E，密度为 $\rho$ 的纤维增强FGM梁在超音速气流中的力学模型，如图1所示，假设温度在厚度方向均匀分布，温度变化为 $\Delta T$ ，在其上表面沿x方向上存在超音速气流，流速为 $U_{\infty }$ ，马赫数为 $M_{\infty }$ ，空气密度为 ${\rho }_{\infty }$ 。

Figure 1. Structure and coordinate diagram of fiber reinforced FGM beam

图1. 纤维增强FGM梁的结构和坐标示意图

纤维增强FGM梁的物性参数可以由下面的混合率公式 [15] 表示：

$P\left(z\right)=P_F\left(1-V_f\right)+P_f{V}_f$ (1)

基体FGM的物性参数由下面的指数模型 [16] 表示：

$P_F\left(z\right)=P_m{\text{e}}^{\beta \left(z+h/2\right)}$ (2)

其中，V为体积分数，下标F表示FGM的基体，c表示陶瓷，m表示金属，f表示纤维， $\beta =\left(\ln \left(P_c/P_m\right)\right)/h$ ，P为弹性模量和密度。

根据经典梁理论，假定温度场是准定常的，只考虑梁的横向振动，则应变位移关系为：

${\epsilon }_x=\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2-z\frac{{\partial }^{2}w}{\partial x^2}$ (3)

应力应变关系为：

${\sigma }_x=E\left(z\right)\left[{\epsilon }_x-\alpha \left(z\right)T\right]$ (4)

纤维增强FGM梁的内力位移关系为：

$N_x={\displaystyle {\int }_A{\sigma }_x\text{d}A}=A_x{\epsilon }_x-N_T=A_x\left[\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\right]-N_T$ (5a)

$M_x={\displaystyle {\int }_{A}z{\sigma }_x\text{d}A}=-D_x{\epsilon }_x-M_T=-D_x\frac{{\partial }^{2}w}{\partial x^2}-M_T$ (5b)

式中： $N_x$ 和 $M_x$ 分别为梁截面上的轴力和弯矩。 $N_T$ 和 $M_T$ 为热薄膜力和热弯矩， $A_x$ 和 $D_x$ 为刚度系数，分别表示如下：

$\left(N_T,M_T\right)={\displaystyle {\int }_{A}E\left(z\right)\alpha T\left(1,z\right)\text{d}A},$ (6a)

$\left(A_x,D_x\right)={\displaystyle {\int }_{A}E\left(z\right)\left(1,z^2\right)\text{d}A},$ (6b)

$A_x=E_{m}bh{g}_1$ , $D_x=E_{m}b{h}^3{g}_2/12$ , $N_T=\alpha T{E}_{m}bh{g}_1$ , $M_T=\alpha T{E}_{m}b{h}^3{g}_2/12$

$g_1=\left(1-V_f\right)\frac{E_c/E_m-1}{\ln \left(E_c/E_m\right)}+\frac{E_f}{E_m}{V}_f$

$g_2=\left(1-V_f\right)\left[\frac{3\left(E_c/E_m-1\right)}{\ln \left(E_c/E_m\right)}-\frac{12\left(E_c/E_m+1\right)}{{\ln }^2\left(E_c/E_m\right)}\right]-\left[\left(1-V_f\right)\frac{24\left(1-E_c/E_m\right)}{{\ln }^3\left(E_c/E_m\right)}-\frac{E_f}{E_m}{V}_f\right]$

利用Hamilton原理推出非线性经典梁下的运动方程：

$\frac{\partial N_x}{\partial x}-\left(I_1\frac{{\partial }^{2}u}{\partial t^2}-I_2\frac{{\partial }^{3}w}{\partial x\partial t^2}\right)=0$ (7a)

$\frac{\partial }{\partial x}\left(N_x\frac{\partial w}{\partial x}\right)+\frac{{\partial }^2{M}_x}{\partial x^2}-I_1\frac{{\partial }^{2}w}{\partial t^2}+\left(-I_2\frac{{\partial }^{3}u}{\partial x\partial t^2}+I_3\frac{{\partial }^{4}w}{\partial x^2\partial t^2}\right)+\Delta p_a=0$ (7b)

其中：

$\left(I_1,I_2,I_3\right)={\displaystyle {\int }_A\rho \left(z\right)\left(1,z,z^2\right)\text{d}A},$ (8)

$I_1={\rho }_{m}bh{G}_1,I_2={\rho }_{m}b{h}^2{G}_2,I_3={\rho }_{m}b{h}^3{G}_3/12$ ,

$G_1=\left(1-V_f\right)\frac{{\rho }_c/{\rho }_m-1}{\ln \left({\rho }_c/{\rho }_m\right)}+\frac{{\rho }_f}{{\rho }_m}{V}_f$

$G_2=\left(1-V_f\right)\left[\frac{{\rho }_c/{\rho }_m-1}{2\ln \left({\rho }_c/{\rho }_m\right)}+\frac{1-{\rho }_c/{\rho }_m}{{\ln }^2\left({\rho }_c/{\rho }_m\right)}\right]$

$G_3=\left(1-V_f\right)\left[\frac{3\left({\rho }_c/{\rho }_m-1\right)}{\ln \left({\rho }_c/{\rho }_m\right)}-\frac{12\left({\rho }_c/{\rho }_m+1\right)}{{\ln }^2\left({\rho }_c/{\rho }_m\right)}\right]-\left[\left(1-V_f\right)\frac{24\left(1-{\rho }_c/{\rho }_m\right)}{{\ln }^3\left({\rho }_c/{\rho }_m\right)}-\frac{{\rho }_f}{{\rho }_m}{V}_f\right]$

$\Delta p_a=-\frac{{\rho }_{\infty }{U}_{\infty }^2}{\sqrt{M_{\infty }^2-1}}\left(\frac{M_{\infty }^2-2}{M_{\infty }^2-1}\frac{1}{U_{\infty }}\frac{\partial w}{\partial t}+\frac{\partial w}{\partial x}\right)$ (9)

对于纤维增强FGM梁的边界采取对边简支，轴向位移边界条件为： $x=0,L$ 时， $u=0$ 。略去面内惯性、转动惯性以及耦合惯性的影响，并消去轴向位移，可以得到气动弹性微分方程，用 $w^{\prime }$ 表示 $\partial w/\partial x$ ， $\stackrel{\dot{}}{w}$ 表示 $\partial w/\partial t$ ，则其无量纲形式可表示如下：

${\bar{w}}^{\left(4\right)}-\left(6f_1{\displaystyle {\int }_0^1{\left({\bar{w}}^{\prime }\right)}^2\text{d}\bar{x}}-R_T\right){\bar{w}}^{″}+f_2\ddot{\stackrel{¯}{w}}+p_1{\bar{w}}^{\prime }+p_2\stackrel{\dot{}}{\bar{w}}=0$ (10)

式中： $f_1=\frac{g_1}{g_2}$ ， $R_T=\frac{N_T{L}^2}{D_x}=P{f}_1$ ， $f_2=\frac{I_1{L}^4}{D_x{T}^2}$ ， $\bar{w}=\frac{w}{h}$ ， $\bar{x}=\frac{x}{L}$ ， $\tau =\frac{t}{T}$ 。

令 $\ddot{\stackrel{¯}{w}}$ 的系数 $f_2$ 为1时，气动刚性系数和气动阻尼系数 $p_1$ 和 $p_2$ 分别为：

$p_1=\frac{L^3{\rho }_{\infty }{U}_{\infty }^2}{D_x\sqrt{M_{\infty }^2-1}}$ , $p_2=\frac{L^2{\rho }_{\infty }{U}_{\infty }\left(M_{\infty }^2-2\right)}{\sqrt{I_{1}D}{\left(M_{\infty }^2-1\right)}^{\frac{3}{2}}}$

纤维增强FGM梁对边简支横向边界条件如下：

$\bar{w}\left(0,\tau \right)=\bar{w}\left(L,\tau \right)=0,\frac{{\partial }^2\bar{w}\left(0,\tau \right)}{\partial {\bar{x}}^2}=\frac{{\partial }^2\bar{w}\left(L,\tau \right)}{\partial {\bar{x}}^2}=0$ (11)

3. Galerkin法求解屈曲和静态分叉问题

纤维增强FGM梁在超音速气流中会产生气动热，在热载荷的作用下引起静气动弹性变形。为研究该模型在对边简支情况下的静气动弹性变形，需将上述控制方程去掉阻尼项、强迫项和时间相关项，用 ${\rm K}\left(\bar{x}\right)$ 表示去掉时间相关项后的 $\bar{w}\left(\bar{x},\tau \right)$ ，作为该模型的屈曲构型，由无量纲微分方程得到：

${{\rm K}}^{\left(4\right)}-\left(6f_1{\displaystyle {\int }_0^L{\left({{\rm K}}^{\prime }\right)}^2\text{d}\bar{x}}-P{f}_1\right){{\rm K}}^{″}=0$ (12)

该模型两端简支的边界条件为：

${{\rm K}|}_{x=0,1}=0$ , ${{{\rm K}}^{″}|}_{x=0,1}=0$ (13)

假设 ${\lambda }^2=P{f}_1-6f_1{\displaystyle {\int }_0^L{\left({{\rm K}}^{\prime }\right)}^2\text{d}\bar{x}}$ ，式(12)可写为：

${{\rm K}}^{\left(4\right)}+{\lambda }^2{{\rm K}}^{″}=0$ (14)

上式通解为：

${\rm K}\left(\bar{x}\right)=C_1\sin \left(\lambda \bar{x}\right)+C_2\cos \left(\lambda \bar{x}\right)+C_3\bar{x}+C_4$ (15)

把边界条件(13)带入式(15)可得：

$\{\begin{array}{l}{C}_2+C_4=0\\ C_1\sin \lambda +C_2\cos \lambda +C_3+C_4=0\\ {\lambda }^2{C}_2=0\\ {\lambda }^2{C}_2\cos \lambda +{\lambda }^2{C}_1\sin \lambda =0\end{array}$ (16)

${\rm K}\left(\bar{x}\right)$ 应具有非零解，所以使齐次线性方程的系数行列式为零，以此得到特征方程：

$C\sin \lambda =0$ 或 $\lambda =m\pi$ ( $m=1,2,3,\cdots$ ) (17)

得到 ${\rm K}=C\sin \left(m\pi \bar{x}\right)$ ，将它带回关于 $\lambda$ 的表达式可得：

$C=\pm \frac{1}{\lambda }\sqrt{\frac{P{f}_1-{\lambda }^2}{3f_1}}$ (18)

从而得到：

${\rm K}=\pm \frac{1}{\lambda }\sqrt{\frac{P{f}_1-{\lambda }^2}{3f_1}}\sin \left(\lambda \bar{x}\right)$ (19)

令纤维增强FGM梁的各项参数为： $L=0.3\text{m}$ ， $b=0.03\text{m}$ ， $h=0.004\text{m}$ ，取其四分之一长度处进行计算，得到在热应力影响下的前三阶弯曲构型的静态分叉图(见图2)。

(a) $V_f=0$ (b) $V_f=0.1$(c) $V_f=0.3$

Figure 2. Static bifurcation diagrams of the first three bending configurations with different fiber volume fractions

图2. 不同的纤维体积分数前三阶弯曲构型的静态分叉图

Table 1. Static bifurcation values of the first three bending configurations with different fiber volume fractions

表1. 不同的纤维体积分数前三阶弯曲构型的静态分叉数值

由图2和表1可知，在 $P<P_1$ 时，即未达到一阶临界屈曲载荷时，系统是稳定的，随着轴向力的增加，达到 $P_1$ 以后，系统产生分岔，发生屈曲，随后变为两种平衡：静平衡构型和屈曲构型。当轴向力大于 $P_2$ 和 $P_3$ ，即达到二阶临界屈曲载荷和三阶临界屈曲载荷时，又会分别增加一个新的屈曲构型。当 $P>P_3$ 时，系统会存在一个静平衡构型和三个屈曲构型。随着纤维的体积分数增大，系统的前三阶弯曲构型的临界屈曲载荷随之减小，这是由于功能梯度材料中增加纤维，改变了梁中的温度应力。

4. Galerkin法求解振动问题

利用Galerkin法，简支边界的位移函数表示为正弦函数的线性组合：

$\bar{w}\left(x,t\right)={\displaystyle \underset{i}{\overset{n}{\sum }}{q}_i\left(\tau \right)}\sin \left(i\pi \bar{x}\right)$ (20)

取二阶Galerkin截断阶数，则 $q=\left(\begin{array}{c}{q}_1\left(\tau \right)\\ q_2\left(\tau \right)\end{array}\right)$ ， $\phi =\left(\begin{array}{c}{\phi }_1\left(\bar{x}\right)\\ {\phi }_2\left(\bar{x}\right)\end{array}\right)$ ， $\bar{w}\left(\bar{x},\tau \right)={\phi }^{\text{T}}q$ ，代入无量纲微分方程得到：

${\left[{\phi }^4\right]}^{\text{T}}q-\left[6f_1{\displaystyle {\int }_0^1{\left({\left[{\phi }^{\prime }\right]}^{\text{T}}q\right)}^2\text{d}\bar{x}}-N_T{f}_2\right]{\left[{\phi }^{″}\right]}^{\text{T}}q+f_3{\phi }^{\text{T}}\ddot{q}+p_1{\left[{\phi }^{\prime }\right]}^{\text{T}}q+p_2{\phi }^{\text{T}}\stackrel{\dot{}}{q}=0$ (21)

方程两端各项左乘 $\phi =\left(\begin{array}{c}{\phi }_1\left(\bar{x}\right)\\ {\phi }_2\left(\bar{x}\right)\end{array}\right)$ ，并在气动弦长区间 $\left[0,1\right]$ 上积分，从而得到：

$f_2\left(\begin{array}{c}{\ddot{q}}_1\\ {\ddot{q}}_2\end{array}\right)+p_2\left(\begin{array}{c}{\stackrel{\dot{}}{q}}_1\\ {\stackrel{\dot{}}{q}}_2\end{array}\right)+\left[\begin{array}{cc}{\pi }^4-R_T{\pi }^2& -\frac{8}{3}{p}_1\\ \frac{8}{3}{p}_1& 16{\pi }^4-4R_T{\pi }^2\end{array}\right]\left(\begin{array}{c}{q}_1\\ q_2\end{array}\right)+f_1\left(\begin{array}{c}{F}_1\\ F_2\end{array}\right)=0$ (22)

$\left(\begin{array}{c}{F}_1\\ F_2\end{array}\right)=\left(\begin{array}{c}3{\pi }^4{q}_1^3+12{\pi }^4{q}_1{q}_2^2\\ 12{\pi }^4{q}_1^2{q}_2+48{\pi }^4{q}_2^3\end{array}\right)$ (23)

下面分析超音速流中受热壁板系统由平衡点失稳而产生分岔的解析边界条件。令 ${\stackrel{\dot{}}{q}}_1=q_3$ ， ${\stackrel{\dot{}}{q}}_2=q_4$ ，则 ${\ddot{q}}_1={\stackrel{\dot{}}{q}}_3$ ， ${\ddot{q}}_2={\stackrel{\dot{}}{q}}_4$ 。于是，二维壁板系统的控制方程化为四维一阶非线性常微分方程组，写成矩阵形式为：

$\left(\begin{array}{c}{\stackrel{\dot{}}{q}}_1\\ {\stackrel{\dot{}}{q}}_2\\ {\stackrel{\dot{}}{q}}_3\\ {\stackrel{\dot{}}{q}}_4\end{array}\right)=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -\frac{{\pi }^4-R_T{\pi }^2}{f_2}& \frac{8}{3f_2}{p}_1& -\frac{1}{f_2}{p}_2& 0\\ -\frac{8}{3f_2}{p}_1& -\frac{16{\pi }^4-4R_T{\pi }^2}{f_2}& 0& -\frac{1}{f_2}{p}_2\end{array}\right)\left(\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\end{array}\right)-\left(\begin{array}{c}0\\ 0\\ \frac{f_1}{f_2}\left(3{\pi }^4{q}_1^3+12{\pi }^4{q}_1{q}_2^2\right)\\ \frac{f_1}{f_2}\left(12{\pi }^4{q}_1^2{q}_2+48{\pi }^4{q}_2^3\right)\end{array}\right)$ (24)

将Hopf分叉转化为非线性方程的求根问题，用以Hopf分叉的代数判断和计算。假设系统的一个平衡点为 $X_0\left(0,0,0,0\right)$ ，该平衡点处的Jacobi矩阵对应的特征方程为：

$\left|\begin{array}{cccc}-\lambda & 0& 1& 0\\ 0& -\lambda & 0& 1\\ -\frac{{\pi }^4-R_T{\pi }^2}{f_2}& \frac{8}{3f_2}{p}_1& -\frac{1}{f_2}{p}_2-\lambda & 0\\ -\frac{8}{3f_2}{p}_1& -\frac{16{\pi }^4-4R_T{\pi }^2}{f_2}& 0& -\frac{1}{f_2}{p}_2-\lambda \end{array}\right|=0$ (25)

从上得到关于 $\lambda$ 的四次实系数方程：

$a_0{\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0$ (26)

$\begin{array}{l}{a}_0=1\\ a_1=\frac{2p_2}{f_2}\\ a_2=\frac{17{\pi }^4{f}_2-5{\pi }^2{f}_2{R}_T+p_2^2}{f_2^2}\\ a_3=\frac{\left(17{\pi }^4-5{\pi }^2{R}_T\right)p_2}{f_2^2}\\ a_4=\frac{16{\pi }^8-20{\pi }^6{R}_T+4{\pi }^4{R}_T^2}{f_2^2}+\frac{64p_1^2}{9f_2^2}\end{array}$ (27)

计算相应的各阶Hurwitz行列式：

$\begin{array}{l}{\Delta }_1=a_1=2p_2\\ {\Delta }_2=\left|\begin{array}{cc}{a}_1& 1\\ a_3& a_2\end{array}\right|=a_1{a}_2-a_3\\ {\Delta }_3=\left|\begin{array}{ccc}{a}_1& 1& 0\\ a_3& a_2& a_1\\ 0& a_4& a_3\end{array}\right|=a_1{a}_2{a}_3-a_1^2{a}_4-a_3^2\\ =p_2^2\left(\frac{225f_2{\pi }^8-90f_2{\pi }^6{R}_T+9{\pi }^4{f}_2{R}_T^2}{f_2^5}+\frac{34{\pi }^4{p}_2^2-10{\pi }^2{R}_T{p}_2^2}{f_2^5}-\frac{256p_1^2}{9f_2^5}\right)\end{array}$ (28)

${\Delta }_1>0$ 且 ${\Delta }_2>0$ ，当 $M_{\infty }\gg 0$ 时， $\sqrt{M_{\infty }^2-1}\to M_{\infty }$

$p_1^2=\frac{L^6{\rho }_{\infty }^2{U}_{\infty }^4}{D_x^2{M}^2}$ , $p_2^2=\frac{L^4{\rho }_{{}_{\infty }}^2{U}_{{}_{\infty }}^2}{I_1{D}_x{M}_{\infty }^2}$ (29)

${\Delta }_3=0$ 可得临界流速为：

$U_{cr}^2=\frac{9}{256}\frac{D_x}{I_1{L}^2}\frac{17{\pi }^4-5{\pi }^2{R}_T}{f_2}+\frac{9}{256}\frac{{\pi }^2{D}_x{M}_{\infty }}{{\rho }_{\infty }{L}^3}\sqrt{\frac{{\rho }_{\infty }^2{L}^2}{{\left(I_1{M}_{\infty }\right)}^2}{\left(\frac{17{\pi }^2-5R_T}{f_2}\right)}^2+{\left(80{\pi }^2+16R_T\right)}^2}$ (30)

$p_1=P_{1cr}=\frac{L^3{\rho }_{\infty }}{D_x{M}_{\infty }}{U}_{cr}^2$ 时，方程有一对纯虚根，可得颤振频率：

${\omega }^2=\frac{{\Delta }_1}{{\Delta }_2}{a}_4$ (31)

设向量 $X=\left(x_1,x_2,x_3,x_4\right)$ 和 $Y={\left(y_1,y_2,y_3,y_4\right)}^{\text{T}}$ 分别是矩阵A的属于特征值 $i\omega$ 的归一化的左、右特征向量。根据 $XA=i\omega A,YA=i\omega A,XA=1$ ，可得：

$X=v{\left(\begin{array}{c}\frac{3f_2}{8p_1}\left[\left(\frac{16{\pi }^4-4R_T{\pi }^2}{f_2}\right)\left(i\omega +\frac{p_2}{f_2}\right)+i\omega {\left(i\omega +\frac{p_2}{f_2}\right)}^2\right]\\ i\omega +\frac{p_2}{f_2}\\ \frac{3f_2}{8p_1}\left[\left(\frac{16{\pi }^4-4R_T{\pi }^2}{f_2}\right)+i\omega \left(i\omega +\frac{p_2}{f_2}\right)\right]\\ 1\end{array}\right)}^{\text{T}}$ (32)

$Y=\left(\begin{array}{c}1\\ \frac{3f_2}{8p_1}\left[\left(\frac{{\pi }^4-R_T{\pi }^2}{f_2}\right)-\left({\omega }^2+\frac{i\omega p_2}{f_2}\right)\right]\\ i\omega \\ \frac{3f_2}{8p_1}\left[i\omega \left(\frac{{\pi }^4-R_T{\pi }^2}{f_2}\right)-\left(i{\omega }^3+{\omega }^2\frac{p_2}{f_2}\right)\right]\end{array}\right)$ (33)

$v=\frac{8p_1{f}_2}{3p_2\left(17{\pi }^4-5{\pi }^2{R}_T-\left(5f_2+1\right){\omega }^2\right)+6i\omega p_2^2}+\frac{8p_1{f}_2}{3i\omega \left(33{\pi }^4{f}_2+{\pi }^4-9{\pi }^2{f}_2{R}_T-{\pi }^2{R}_T-3{\omega }^2{f}_2^2-{\omega }^2{f}_2\right)}$ (34)

因为：

$B={\frac{\text{d}A\left(P_1\right)}{\text{d}{P}_1}|}_{P_1=P_{1}cr}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& \frac{8}{3}& 0& 0\\ -\frac{8}{3}& 0& 0& 0\end{array}\right)$ (35)

在参数 $P_1=P_{1cr}=\frac{L^3{\rho }_{\infty }}{D_x{M}_{\infty }}{U}_{cr}^2$ 附近，纯虚根 $i\omega$ 可以写为：

${\zeta }_{1,2}\left(p_1\right)=\alpha \left(p_1\right)\pm i\omega \left(p_1\right)$ (36)

其中 $\alpha \left(p_1\right)=0,\omega \left(p_1\right)>0$ ，于是 $X\left(p_1\right)A\left(p_1\right)Y\left(p_1\right)=\alpha \left(p_1\right)\pm i\omega \left(p_1\right)$ ，得到：

${\alpha }^{\prime }\left(p_1\right)=\mathrm{Re}\left(X\left(p_{1cr}\right)A\left(p_{1cr}\right)Y\left(p_{1cr}\right)\right)$ (37)

通过计算得到， ${\alpha }^{\prime }\left(p_1\right)$ 在 $p_1=p_{1cr}=\frac{L^3{\rho }_{\infty }}{M_{\infty }}{U}_{cr}^2\ne 0$ ，并且满足：

$\begin{array}{l}{a}_1>0,a_2>0,\cdots ,a_n>0\\ {\Delta }_{n-1}=0,{\Delta }_i>0\left(i=n-3,n-5,\cdots \right)\end{array}$ (38)

此时系统发生Hopf分岔，即在 $U=U_{cr}$ 时，超音速流中的纤维增强FGM梁会发生颤振。

4.1. 数值算例

假设飞行器的飞行高度为11 km，此时空气密度为0.364 kg/m³，音速为295.065 m/s，马赫数为5。碳纤维、陶瓷和金属的物性参数如下： $E_f=400\text{GPa}$ ， ${\rho }_f=1800\text{kg}/{\text{m}}^{\text{3}}$ ， $E_c=3800\text{GPa}$ ， ${\rho }_c=3800\text{kg}/{\text{m}}^{\text{3}}$ ， $E_m=70\text{GPa}$ ， ${\rho }_m=2780\text{kg}/{\text{m}}^{\text{3}}$ 。取相同的材料物性参数，得到不同温度应力、不同纤维体积分数情况下的无量纲临界流速和无量纲临界频率。

(a) (b)

Figure 3. Effect of fiber volume fraction on dimensionless critical frequency and dimensionless critical flow rate under different temperature stresses

图3. 不同温度应力下，纤维体积分数对无量纲临界频率和无量纲临界流速的影响

图3反映了当温度应力 $R_T={\pi }^2/4$ 、 $R_T=3{\pi }^2/4$ 、 $R_T=5{\pi }^2/4$ 三种情况时，不同纤维体积分数情况下，无量纲临界频率和无量纲临界流速的情况。结果显示；温度应力一定时，随着纤维体积分数 $V_f$ 的增加，梁系统发生颤振的无量纲临界流速呈现上升趋势；而无量纲临界频率仅以非常微小的变化呈现上升趋势时，影响并不明显。取纤维体积分数 $V_f$ 为一个定值时，随着温度应力的增加，无量纲临界流速和无量纲临界频率都在明显变小。

4.2. 稳定性验证

当 $R_T=3{\pi }^2/4$ 、 $p_1=235.5$ 时，初值取(0.005, 0.005, 0.01, 0.1)，从而得到 $V_f$ 分别为0、0.1、0.3时系统的时间历程图和相位图(见图4)。

Table 2. When R T = 3 π 2 / 4 , the eigenvalue of Jacobi matrix for different fiber volume fractions

表2. 当 $R_T=3{\pi }^2/4$ 时，不同的纤维体积分数的Jacobi矩阵特征值

(a) V_f = 0时，系统的时间历程图和相位图 (b) V_f = 0.1时，系统的时间历程图和相位图 (c) V_f = 0.3时，系统的时间历程图和相位图

Figure 4. Time history diagrams and phase diagrams of the system with different fiber volume fractions

图4. 不同纤维体积分数时系统的时间历程图和相位图

表2为在温度应力 $R_T=3{\pi }^2/4$ 时，不同纤维体积分数下所求得的Jacobi矩阵特征值。 $p_1=235.41$ 、 $p_1=235.53$ 、 $p_1=235.81$ 分别为三种纤维体积分数情况下梁发生颤振时的无量纲临界来流动压，此时特征值存在一对纯虚根，并且其余特征值都具有负实部，系统处于稳定状态。而在一对特征值存在负实部，另一对特征值实部为正的情况下，将会在平衡点附近产生稳定极限环，一直振动，从而影响系统寿命。由图3可以看出，在 $p_1=235.5$ 时， $V_f=0.1$ 和 $V_f=0.3$ 的情况下，时间历程图能够趋于一个值，处于稳定状态，而 $V_f=0$ 时处于一直抖动的状态，这刚好验证了前面得到的特征值状况。在温度应力一定时，纤维增强FGM梁在 $V_f=0$ 的情况下，发生颤振的无量纲临界来流动压最小，说明功能梯度材料中增加一定的纤维增强，更加不容易发生颤振。

5. 结论

本文基于经典梁理论，研究了纤维增强FGM梁在超音速流中的稳定性问题。推导得到了纤维增强FGM梁在超音速流和热载荷共同作用下的气动弹性微分方程，讨论了在不同温度应力、不同纤维体积分数的情况下，无量纲临界频率和无量纲临界流速的变化规律以及不同纤维体积分数下梁系统的颤振情况，得到如下结论：

1) 在两边简支的边界条件下，纤维增强FGM梁的无量纲临界载荷随着纤维体积分数的增大呈现变小趋势。

2) 在温度应力一定时，无量纲临界流速随着纤维体积分数的增大而增大，无量纲临界频率虽然也随着纤维体积分数的增大而增大，不过影响非常小。当纤维体积分数一定时，随着温度应力的增加，无量纲临界流速和无量纲临界频率都表现为下降趋势。

3) 在不同纤维体积分数下，系统发生颤振的无量纲来流动压不同。在功能梯度材料中存在纤维增强时要比不存在纤维增强时，能够承受更大的无量纲来流动压。在相同情况下，纤维增强FGM梁中所含纤维体积分数越大，越不容易发生颤振。