功能梯度材料(Functionally gradient materials，FGM)是由两种或两种以上材料组成的复合材料^[1]。它可以被设计成一种非均匀材料，其结构沿优选方向连续变化，其物理性质的分布根据体积分数定律进行分级。同传统复合材料相比，一个显著的区别是功能梯度材料可以减少应力集中^[2-3]，增加粘结强度^[4]，提高断裂韧性^[5]，并提高耐腐蚀性^[6]。此外，FGM所需的机械特性，包括材料密度、泊松比和杨氏模量，可以通过调整沿优选方向的体积分数变化来调整^[7-8]。随着材料制造技术的进步，复合材料的制造方式也越来越多样化，其中烧结法是制备功能梯度材料的一种灵活而合适的方法。然而，在使用无压烧结技术^[9]或多步骤顺序填充技术^[10]制造FGM的过程中，材料内部不可避免地会出现孔隙和微孔。在这里，孔隙度对管道的影响有两个方面：一方面是作为一个缺点，由于制造FGM的制造技术问题或质量差，它可能会削弱结构的力学性能；另一方面是可以优化孔隙度的适当分布，以实现改进的结构和机械性能，包括有效的能量耗散和所需的强度和刚度^[11]。因此，含孔隙流体输送功能梯度材料管道在核反应堆、飞机工程等工程领域有着广泛的应用。

近几十年来，功能梯度材料在热环境中的工程应用让人类越来越重视功能梯度材料的物理性质。在热环境中，材料的力学性能(如弹性模量、热膨胀系数)并不稳定，会随着温度的改变而改变，所以热载荷下的热屈曲问题的研究十分重要。黄怀伟等^[12]对不同梯度温度场下的功能梯度薄壁管进行了热屈曲分析，分析了温度场参数、材料参数对临界温度的影响。Li等^[13]利用边界层理论对温度场中承受外压和轴压耦合载荷的层合圆柱壳的屈曲进行了研究，考虑非线性初始屈曲及几何缺陷，采用奇异摄动法获得数值解。徐凯^[14]对功能梯度材料厚壁圆柱壳进行了热屈曲分析，基于Reddy高阶剪切理论，对不同参数条件的圆柱壳进行了热屈曲分析。李亚杯等^[15]利用Rayleigh-Ritz法对复合材料薄壁管的热屈曲和振动问题进行了研究。侯子豪等^[16]确定了碳纳米管在湿热环境中的后屈曲路径。聂聪聪等^[17]研究了纤维起始角和终止角的变化对管的屈曲临界荷载的影响。刘子赫等^[18]研究了管道保温层变化对管道屈曲能力的影响。武行等^[19]研究了不同缺陷位置及尺寸对管道屈曲压力的影响。王小龙等^[20]研究了含缺陷管道的屈曲失稳临界外压的拐点方法。李学松等^[21]研究了热环境下对圆弧拱的非线性屈曲的影响。龚顺风等^[22]研究了外压作用下含缺陷管道的屈曲失稳，利用模拟实验，测得屈曲失稳时的压力和形变。诸多学者对结构稳定性的兴趣导致屈曲失稳成为连续结构中的一个吸引人的话题。同时，结构元件的后屈曲强度对于设计和分析各种工程结构具有重要意义。

对于多孔的功能梯度材料，目前研究人员主要研究的是关于线性和非线性响应以及稳定性。而且对于功能梯度材料在温度场中的热稳定性、屈曲响应均有研究。但是，在温度场中，对存在缺陷的功能梯度管道，如几何缺陷、制造缺陷等研究还十分匮乏。在工程中，研究在温度场中存在缺陷的功能梯度管道的振动特性，对工程指导具有重要的工程意义。为了更好的预测管道系统的力学性能，需要对这一问题进行深入的讨论。

本文考虑了由中性面拉伸而引起几何非线性的现象，故控制方程存在非线性项。基于Euler-Bernouilli理论，建立了在温度场中的多孔功能梯度管道的数学模型。首先研究了孔隙度体积分数、幂律指数和温度对屈曲的静态平衡的影响；其次研究了伪非线性动力学模型的动态响应；最后通过数值解，讨论了孔隙度体积分数、幂律指数和温度等重要参数对模型振动分析的影响。

1 动力模型的建立

图 1所示为多孔功能梯度管道热-流-固耦合系统的示意图，该管道的长度为L，内径为R_i，外径为R_o的直管。在均匀流速U内，管道两端铰接。假设材料成分的微观结构形态仅在厚度方向上连续变化，但是由于加工缺陷，横截面内存在均匀的空隙。同温度相关的材料特性如表 1和图 1所示。

建立笛卡尔坐标系，引入孔隙度α(α≤1)，空隙均匀分布在陶瓷相和金属相之间的FGM管道的材料的修正系数可以表示为^[23]：

$ P_{\mathrm{f}}=P_{\mathrm{m}}\left(V_{\mathrm{m}}-\frac{\alpha}{2}\right)+P_{\mathrm{c}}\left(V_{\mathrm{c}}-\frac{\alpha}{2}\right) 。$

(1)

式中：P_m，P_c是金属和陶瓷的材料特性；α是孔隙度；V_m和V_c分别是碱金属和陶瓷的体积分数，成分表示为^[24]:

$ V_{\mathrm{c}}+V_{\mathrm{m}}=1。$

(2)

陶瓷相的体积分数可以表示为：

$ V_{\mathrm{c}}=\left(\frac{R_{\mathrm{o}}-r}{R_{\mathrm{o}}-R_{\mathrm{i}}}\right)^n。$

(3)

式中：r为管道横截面任意一点半径；n为功能梯度材料幂律指数。

金属和陶瓷材料的与温度相关的材料特性(杨氏模量E、热膨胀系数δ、导热系数k、密度ρ)可通过Touloukian模型得到^[25]：

$ P_ {\rm{f}}(T)=P_0\left(P_{-1} T^{-1}+1+P_1 T^1+P_2 T^2+P_3 T^3\right) 。$

(4)

式中：P_f代表功能梯度材料的杨氏模量E、热膨胀系数δ、传导率K以及密度ρ；P_－1、P₀、P₁、P₂、P₃是温度系数；T是以开尔文为单位的温度。

基于von-Kármán非线性理论^[26], 假设管道大振动的应变-位移关系可以写成

$ \varepsilon_{x x}=\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2。$

(5)

式中：ε_xx为应变；u和w分别为管道中间平面的轴向积横向位移。

基于Euler-Bernoulli梁理论，功能梯度管道的位移场定义如下：

$ u_1(x, z, t)=u(x, t)-z \frac{\partial w(x, t)}{\partial x}, $

(6)

$ u_2(x, z, t)=0, $

(7)

$ u_3(x, z, t)=w(x, t)。$

(8)

根据Euler-Bernoulli理论，将式(6)~(8)代入式(5)中可以得到：

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

(9)

在温度场中，温度产生的热应变为：

$ \varepsilon_T=\delta(r, T)\left(T(r)-T_0\right)。$

(10)

对于温度升高的功能梯度材料，根据式(4)应力-应变可以表示为：

$ \sigma_{x x}=E(r, T) \varepsilon(x)-E(r, T) \delta(r, T)\left(T(r)-T_0\right) 。$

(11)

功能梯度管道的应变能u₁可以表示为：

$ U_1=\int_V \sigma_{x x}\left(\varepsilon_{x x}-\varepsilon_T\right) \mathrm{d} V 。$

(12)

将式(12)代入，可以得到：

$ \begin{gathered} U_1=\frac{1}{2} \int_0^L\left\{A_{11}\left[\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2\right]^2+D_{11}\left(\frac{\partial^2 w}{\partial x^2}\right)^2-\right. \\ \left.2 N_{\mathrm{T}}\left[\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2\right]+N_{0, \mathrm{~T}}\right\} \mathrm{d} x 。\end{gathered} $

(13)

式中：

$ A_{11}=\int_0^{2 \pi} \int_{R_{\mathrm{i}}}^{R_{\mathrm{o}}} E(r, T) r \mathrm{~d} r \mathrm{~d} \theta, $

(14)

$ D_{11}=\int_0^{2 \pi} \int_{R_{\mathrm{i}}}^{R_ {\rm{o}}} E(r, T) r^2 \sin ^2(\theta) r \mathrm{~d} r \mathrm{~d} \theta, $

(15)

$ N_{\mathrm{T}}=\int_0^{2 \pi} \int_{R_{\mathrm{i}}}^{R_{\mathrm{o}}} E(r, T) \delta(r, T)\left(T(r)-T_0\right) r \mathrm{~d} r \mathrm{~d} \theta, $

(16)

$ N_{0, \mathrm{~T}}=\int_0^{2 \pi} \int_{R_{\mathrm{i}}}^{R_{\rm{o}}} E(r, T)\left(\delta(r, T)\left(T(r)-T_0\right)\right)^2 r \mathrm{~d} r \mathrm{~d} \theta \text { 。} $

(17)

在该系统发生振动时，管道发生弯曲变形而伸长，会有附加的轴向力。因弯曲而诱导发生的应变能u₂可以表示为：

$ \begin{array}{*{20}{c}} {U_2=\int_0^L \frac{P_0}{A}\left(\int_A \varepsilon_{x x} \mathrm{~d} A\right) \mathrm{d} x=}\\ {\frac{1}{2} \int_0^L P_0\left[2 \frac{\partial u}{\partial x}+\left(\frac{\partial w}{\partial x}\right)^2\right] \mathrm{d} x 。} \end{array} $

(18)

式中P₀为初始轴向力。

系统总的应变能U可以表示为：

$ U=U_1+U_2 \text { 。} $

(19)

功能梯度输流管的总动能可以表示为：

$ T=T_{\mathrm{f}}+T_{\mathrm{p}}。$

(20)

式中：T_f为管道内流体产生的动能；T_p为管道产生的动能。具体为：

$ \begin{aligned} T_{\mathrm{f}}= & \frac{1}{2} m_{\mathrm{f}} \int_0^L\left(v+\frac{\partial u}{\partial t}+v \frac{\partial u}{\partial x}\right)^2+ \\ & \left.\left(\frac{\partial w}{\partial t}+v \frac{\partial w}{\partial x}\right)^2\right] \mathrm{d} x, \end{aligned} $

(21)

$ T_{\mathrm{p}}=\frac{1}{2} m_{\mathrm{p}} \int_0^L\left[\left(\frac{\partial u}{\partial t}\right)^2+\left(\frac{\partial w}{\partial t}\right)^2\right] \mathrm{d} x 。$

(22)

式中：v表示管道内流体的流速；m_f、m_p分别表示管道和流体单位长度的质量，且

$ m_{\mathrm{f}}=\rho_{\mathrm{f}} A_{\mathrm{f}}, $

(23)

$ m_{\mathrm{p}}=\int_A \rho(r) \mathrm{d} {A}。$

(24)

根据Hamilton原理：

$ \delta \int_{t_1}^{t_2}(T-U) \mathrm{d} t=0 。$

(25)

将式(20)和(21)带入到式(23)和(24)中，然后进行分部积分，就可以得到控制方程：

$ \begin{gathered} \delta u:-\frac{\partial}{\partial x}\left[A_{11}\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2\right)\right]-\frac{\partial P_0}{\partial x}+m_{\mathrm{f}} v^2 \frac{\partial^2 u}{\partial x^2}+ \\ 2 m_{\mathrm{f}} v \frac{\partial^2 u}{\partial x \partial t}+\left(m_{\mathrm{f}}+m_{\mathrm{p}}\right) \frac{\partial^2 u}{\partial t^2}=0, \end{gathered} $

(26)

$ \begin{gathered} \delta w: D_{11} \frac{\partial^4 w}{\partial x^4}+\left(m_{\mathrm{f}} v^2+N^{\mathrm{T}}\right) \frac{\partial^2 w}{\partial x^2}- \\ \frac{\partial}{\partial x}\left[\left(A_{11}\left(\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)\right)^2\right)+P_0\right) \frac{\partial w}{\partial x}\right]+ \\ C_{\mathrm{d}} \frac{\partial w}{\partial t}+2 m_{\mathrm{f}} v \frac{\partial^2 w}{\partial x \partial t}+\left(m_{\mathrm{f}}+m_{\mathrm{p}}\right) \frac{\partial^2 w}{\partial t^2}=0 。\end{gathered} $

(27)

边界条件在(x=0, L)可以表示为：

$ \begin{gathered} \delta u: A_{11}\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2\right)-P_0+ \\ m_{\mathrm{f}} v^2\left(1+\frac{\partial u}{\partial x}\right)+2 m_{\mathrm{f}} v \frac{\partial u}{\partial t}=0 。\end{gathered} $

(28)

$ \begin{gathered} \delta w: D_{11} \frac{\partial^3 w}{\partial x^3}+\left(m_{\mathrm{f}} v^2+N^T\right) \frac{\partial w}{\partial x}- \\ A_{11}\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2\right) \frac{\partial w}{\partial x}- \\ P_0 \frac{\partial w}{\partial x}+2 m_{\mathrm{f}} v \frac{\partial w}{\partial t}=0 。\end{gathered} $

(29)

对于具有2个简单支承的FGM管道，满足的边界条件为：

$ u(0)=u(L)=0, $

(30)

$ w(0)=w(L)=0, $

(31)

$ w^{\prime \prime}(0)=w^{\prime \prime}(L)=0 。$

(32)

在本文的研究中，假设功能梯度管的振动是一个有限拉伸的问题，只考虑横向运动时，可以得到：

$ \frac{\partial}{\partial x}\left[\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial w}{\partial x}\right)^2\right]=0。$

(33)

对式(33)进行2次积分可以得到纵向位移的表达式：

$ u=x C_1(t)-\frac{1}{2} \int_0^x\left(\frac{\partial w}{\partial x}\right)^2 \mathrm{~d} x+C_2(t) 。$

(34)

根据简支输流管的边界条件，可以得到：

$ C_1(t)=\frac{1}{2 L} \int_0^L\left(\frac{\partial w}{\partial x}\right)^2 \mathrm{~d} x, C_2(t)=0 。$

(35)

$ u=\frac{x}{2 L} \int_0^L\left(\frac{\partial w}{\partial x}\right)^2 \mathrm{~d} x-\frac{1}{2} \int_0^x\left(\frac{\partial w}{\partial x}\right)^2 \mathrm{~d} x 。$

(36)

将式(36)代入到式(29)中可以得到横向振动的控制方程为：

$ \begin{gathered} \left(m_{\mathrm{f}}+m_{\mathrm{p}}\right) \frac{\partial^2 w}{\partial t^2}+2 m_{\mathrm{f}} v \frac{\partial^2 w}{\partial x \partial t}+m_{\mathrm{f}} v^2 \frac{\partial^2 w}{\partial x^2}+ \\ c \frac{\partial w}{\partial t}+D_{11} \frac{\partial^4 w}{\partial x^4}+\left(N^T-P_0\right) \frac{\partial^2 w}{\partial x^2}- \\ \frac{A_{11}}{2 L} \frac{\partial^2 w}{\partial x^2} \int_0^L\left(\frac{\partial w}{\partial x}\right)^2 \mathrm{~d} x=0 。\end{gathered} $

(37)

将控制方程无量纲化，引入无量纲参数：

$ \begin{gathered} \xi=\frac{\dot{x}}{L}, \eta=\frac{w}{L}, \tau=t \sqrt{\frac{D_{11}}{\left(m_{\mathrm{f}}+m_{\mathrm{p}}\right) L^4}}, \\ \alpha=\left(\frac{D_{11}}{m}\right)^{\frac{1}{2}} \frac{\mu}{L^2}, \Pi=\frac{A_{11} L^2}{D_{11}}, u=\sqrt{\frac{m_{\mathrm{f}}}{D_{11}}} v L, \\ {\mathit{\Lambda}}=\frac{N_{\mathrm{T}} L^2}{D_{11}}, P=\frac{P_0 L^2}{D_{11}}, \beta=\sqrt{\frac{m_{\mathrm{f}}}{m_{\mathrm{f}}+m_{\mathrm{p}}}}。\end{gathered} $

(38)

通过引入无量纲参数，无量纲控制方程可以表示为：

$ \begin{gathered} \frac{\partial^2 \eta}{\partial \tau^2}+2 u \sqrt{\beta} \frac{\partial^2 \eta}{\partial \xi \partial \tau}+u^2 \frac{\partial^2 \eta}{\partial \xi^2}+\frac{\partial^4 \eta}{\partial \xi^4}+ \\ ({\mathit{\Lambda}}-P) \frac{\partial^2 \eta}{\partial \xi^2}-\frac{\Pi}{2} \frac{\partial^2 \eta}{\partial \xi^2} \int_0^L\left(\frac{\partial \eta}{\partial \xi}\right)^2 \mathrm{~d} \xi=0 。\end{gathered} $

(39)

2 传热模型

考虑内表面温度为T_a，外表面温度为T_o的非理想功能梯度材料管。通过求解一维热传导问题，得到了功能梯度材料管截面的温度分布。FGM管稳态热传导问题的热边界条件和微分方程可表示为^[27]：

$ \frac{1}{r} \frac{\mathrm{d}}{\mathrm{d} r}\left(r k(r, T) \frac{\mathrm{d} T}{\mathrm{~d} r}\right)=0, $

(40)

$ T\left(r_{\mathrm{a}}\right)=T_i, T\left(r_{\mathrm{b}}\right)=T_{\mathrm{o}}。$

(41)

式中k(r, T)由式(4)确定。将边界条件(41)代入式(40)中可以得到：

$ \left\{\begin{array}{l} T\left(r_{\mathrm{a}}\right)=C_1 \int \frac{\mathrm{d} r}{r k\left(r_{\mathrm{a}}, T\right)}+C_2 。\\ T\left(r_{\mathrm{b}}\right)=C_1 \int \frac{\mathrm{d} r}{r k\left(r_{\mathrm{b}}, T\right)}+C_2 。\end{array}\right. $

(42)

式中：r表示管道横截面任意一点的半径; r_a和r_b分别表示材料管的内、外半径; C₁和C₂是2个积分常数，由FGM管内、外表面的热边界条件确定。因此，根据热边界条件，可获得温度分布函数：

$ T(r)=T_{\mathrm{a}}+\frac{T_{\mathrm{o}}-T_{\mathrm{i}}}{\int_{R_{\mathrm{i}}}^{R_{\mathrm{o}}} \int \frac{\mathrm{d} r}{r k(r, T)}} \int_{R_{\mathrm{i}}}^r \frac{1}{r k(r, T)} \mathrm{d} r。$

(43)

进而可以得到热合力的表达式为：

$ \begin{array}{*{20}{c}} {N_{\mathrm{T}}=\frac{2 \pi\left(T_{\mathrm{o}}-T_{\mathrm{i}}\right)}{\int_{R_{\mathrm{i}}}^{R_{\mathrm{o}}} \frac{\mathrm{d} r}{r k(r, T)}} \int_{R_{\mathrm{i}}}^{R_{\mathrm{o}}} {\mathit{\Psi}} \mathrm{d} r+}\\ {2 \pi\left(T_{\mathrm{i}}-T_0\right) \int_{R_{\mathrm{i}}}^{R_{\mathrm{o}}} r E(r, T) \delta(r, T) \mathrm{d} r, } \end{array} $

(44)

$ {\mathit{\Psi}}=r E(r, T)\left[\int_{R_{\mathrm{i}}}^r \frac{\mathrm{d} r}{r k(r, T)}\right] \delta(r, T) 。$

(45)

3 后屈曲静态特性

众所周知，如果流体速度接近临界值，则输送流体的简单支承管道会出现发散不稳定或屈曲^[28]。为了确定后屈曲行为，通过删除式(39)中同时间相关的项来分析静态问题。因此，结果可以写成：

$ \begin{gathered} u^2 \frac{\partial^2 \eta}{\partial \xi^2}+\frac{\partial^4 \eta}{\partial \xi^4}+({\mathit{\Lambda}}-P) \frac{\partial^2 \eta}{\partial \xi^2}- \\ \frac{\Pi}{2} \frac{\partial^2 \eta}{\partial \xi^2} \int_0^1\left(\frac{\partial \eta}{\partial \xi}\right)^2 \mathrm{~d} \xi=0 。\end{gathered} $

(46)

对于给定的η(ξ)，$\int_0^1\left(\frac{\partial \eta}{\partial \xi}\right)^2 d \xi$是一个定值。故可令：

$ \vartheta=\frac{\Pi}{2} \int_0^1\left(\frac{\partial \eta}{\partial \xi}\right)^2 \mathrm{~d} \xi \text { 。} $

(47)

屈曲方程可以表示为：

$ \frac{\partial^4 \eta}{\partial \xi^4}+\lambda^2 \frac{\partial^2 \eta}{\partial \xi^2}=0 。$

(48)

式中：

$ \lambda^2=u^2+{\mathit{\Lambda}}-P-\vartheta 。$

(49)

式(4)的通解形式可以写为：

$ \eta(\xi)=C_1+C_2 \xi+C_3 \cos (\lambda \xi)+C_4 \sin (\lambda \xi) 。$

(50)

边界条件可以表示为：

$ \begin{gathered} \eta(0, t)=\eta(1, t)=0, \\ \frac{\partial^2 \eta(0, t)}{\partial \xi^2}=\frac{\partial^2 \eta(1, t)}{\partial \xi^2}=0 。\end{gathered} $

(51)

将边界条件代入到方程(48)中，可以的得到：

$ \begin{array}{*{20}{c}} {C_1=C_2=C_3=0, }\\ {C_4 \sin (\lambda \xi)=0, \lambda=n \pi, n=1, 2, 3 \cdots \cdots}。\end{array} $

(52)

将式(9)带回到式(4)中可以得到：

$ C_4= \pm 2\left[\frac{\left(u^2+{\mathit{\Lambda}}-P\right)-(n \pi)^2}{\Pi(n \pi)^2}\right]^{0.5} 。$

(53)

进而得到后屈曲的精确表达式：

$ \eta(\xi)= \pm 2\left[\frac{\left(u^2+{\mathit{\Lambda}}-P\right)-(n \pi)^2}{\Pi(n \pi)^2}\right]^{0.5} \sin (n \pi x) 。$

(54)

4 屈曲后的动力学行为

为研究屈曲条件下的动力学响应，假设位移的表达式为：

$ \varphi(\xi, t)=\eta(\xi)+\varphi(\xi, t) 。$

(55)

式中φ(ξ, t)表示在屈曲位形η(ξ)附近的扰动。将式(55)带入到方程(39)中，可以得到输流管绕非平衡位置的控制方程：

$ \begin{gathered} \frac{\partial^2 \varphi}{\partial t^2}+2 \sqrt{\beta} u \frac{\partial^2 \varphi}{\partial \xi \partial t}+\left(u^2+{\mathit{\Lambda}}-P\right) \frac{\partial^2 \varphi}{\partial \xi^2}+\frac{\partial^4 \varphi}{\partial \xi^4}- \\ \frac{\Pi}{2} \frac{\partial^2 \varphi}{\partial \xi^2} \int_0^1\left(\frac{\partial \eta}{\partial \xi}\right)^2 \mathrm{~d} \xi-\Pi \frac{\partial^2 \eta}{\partial \xi^2} \int_0^1\left(\frac{\partial \varphi}{\partial \xi} \frac{\partial \eta}{\partial \xi}\right) \mathrm{d} \xi- \\ \frac{\Pi}{2} \frac{\partial^2 \varphi}{\partial \xi^2} \int_0^1\left(\frac{\partial \varphi}{\partial \xi}\right)^2 \mathrm{~d} \xi-\frac{\Pi}{2} \frac{\partial^2 \eta}{\partial \xi^2} \int_0^1\left(\frac{\partial \varphi}{\partial \xi}\right)^2 \mathrm{~d} \xi- \\ \Pi \frac{\partial^2 \varphi}{\partial \xi^2} \int_0^1\left(\frac{\partial \varphi}{\partial \xi} \frac{\partial \eta}{\partial \xi}\right) \mathrm{d} \xi 。\end{gathered} $

(56)

利用式(56)去掉方程中的非线性项阻尼项，可以得到后屈曲伪非线性控制方程:

$ \begin{gathered} \frac{\partial^2 \varphi}{\partial t^2}+2 \sqrt{\beta} u \frac{\partial^2 \varphi}{\partial \xi \partial t}+\left(u^2+{\mathit{\Lambda}}-P\right) \frac{\partial^2 \varphi}{\partial \xi^2}+\frac{\partial^4 \varphi}{\partial \xi^4}- \\ \frac{\Pi}{2} \frac{\partial^2 \varphi}{\partial \xi^2} \int_0^1\left(\frac{\partial \eta}{\partial \xi}\right)^2 \mathrm{~d} \xi-\Pi \frac{\partial^2 \eta}{\partial \xi^2} \int_0^1\left(\frac{\partial \varphi}{\partial \xi} \frac{\partial \eta}{\partial \xi}\right) \mathrm{d} \xi=0 \end{gathered} $

(57)

由文献[29]可知，简支梁屈曲后的模态函数与屈曲前是相同的。采用伽辽金截断法来离散控制方程，假设：

$ \varphi(\xi, t)=\sum\limits_{i=1}^n \varphi_i(\xi) \psi_i(t), $

(58)

$ \eta(\xi)=\mu_1 \varphi_1(\xi) 。$

(59)

5 数值分析与讨论

为了验证该方法的正确性与准确性，本文作者将其同现有文献进行了对比验证。同Dai等^[30]的工作相比，图 2给出了流体速度变化时各向同性管道的前四阶屈曲前和屈曲后频率。可以看出，目前的结果同文献[31]中报告的结果非常一致，这接受了目前模型的准确性。此外，临界流体速度等于3.14, 这同文献[32]的结果一致。

图 3显示了管道的前3个屈曲配置的静态分叉图。根据无量纲流速绘制四分之一跨度处管道的无量纲横向位移。在没有达到临界流速时，管道在其原始静态平衡位置稳定。当轴向速度达到π时，第一次达到屈曲的不稳定性临界值通过分支点分叉发生，导致挠度位移突然上升。超过这一点，可以看到横向位移的2个可能的稳定解和一个不稳定解。当轴向速度增加到超过对应于第二临界屈曲的第二临界值2π时，管道有3种平衡状态：不稳定的直线结构、对应于第一和第二屈曲结构的另外两种平衡状态。当流速超过临界流速时，屈曲幅度增大。随着轴向速度增加超过第三临界屈曲值时，管道呈现出与3种屈曲构型相关的3种非平凡平衡。对应于图 3所示的3种屈曲配置。

图 4显示了在初始轴向力的作用下，管道的前3个屈曲配置的静态分叉图。在未达到初始临界值时，管道在其原始静态平衡位置稳定。当初始轴向力达到第一临界值时，未偏转位置变得不稳定。超过该临界载荷后，直线构型变得不稳定，输流管获得其他稳定平衡位置，即已知的屈曲构型。当第二次临界屈曲后轴向载荷增加时，管道具有3个平衡。当轴向荷载继续增加超过第三临界屈曲荷载时，管道表现出同屈曲配置相对应的3种非平凡平衡。由此可见，无量纲挠度随着初始轴向力的增加而增加。

图 5研究了FGM管道在不同幂律指数和孔隙度参数下的分叉图。首次屈曲时的无量纲横向位移随速度和初始轴力的变大而变大。FGM管道的屈曲配置如图 5所示，显示了不同的幂律指数和孔隙度参数的结果。幂律指数和孔隙度参数对临界值的影响是显而易见的。对于幂律指数，当所考虑的控制参数较低时，FGM管道呈现出稳定且未发生分叉的静态平衡。然而，当参数值超过一定值时，管道在经历叉式分叉后又表现出双稳态结构和后屈曲。此外，当幂律指数逐渐增大时，FGM管道的后屈曲位移增加，表明陶瓷组件到金属组件的过渡后屈曲行为。值得注意的是，幂律指数变化的大小影响后屈曲响应的绝对值，这表明仅通过控制幂律指数就可以获得所需要的响应。随着幂律指数增大，同管道不稳定性相关的控制参数显著降低。这意味着随着幂律指数的增加，管道在较低的控制参数下发生屈曲。

由表 1可知，SUS304的杨氏模量远小于Si₃N₄的杨氏模量，这表明当幂律指数逐渐增大时，材料由Si₃N₄组件向SUS304组件过渡，杨氏模量逐渐减小，导致材料的刚度降低，从而使得管道的抗屈曲能力减小。

根据式(45)，将流速和初始轴向力作为物理控制参数。图 6显示了不同孔隙度参数下FGM管道屈曲行为的分叉图。可以观察到，这样的控制参数可以改变临界值。故FGM管道在控制参数较大时表现出稳定和未变形的静态平衡状态；然而，当流速超过第一临界值时，管道在经历叉式分叉后呈现双稳态结构和后屈曲。

在后屈曲状态下，图 7为第一阶屈曲配置周围无量纲自然频率随流体流速和初始轴向力的变化。在后屈曲区域，第一自然频率随着控制参数的改变而改变，而第二、三自然频率基本保持不变，即流体流速和初始轴向力在后屈曲状态下对高阶的频率影响比较小。

图 8(a)显示了流体流速和初始轴向力对第一阶屈曲配置周围的第一阶无量纲频率的影响。初始轴向力从零开始增大时(此时为拉力)，使系统发生失稳的临界值也随着初始轴向力的增加而增加。不同幂律指数和孔隙度下流体流速第一阶无量纲频率的影响。在后屈曲区域中，流速较低时，幂律指数对频率的影响较小，而孔隙度参数对频率的影响在低流速比高流速更明显。图 8(c)显示了温度对频率的影响，当输流管存在内外温差时，就会产生热轴力。由图 8(c)中可知，随着内外温差的加大，使得在较小流速下，系统就会发生失稳。而随着流速的增加，系统频率也会迅速增加。

为了展示孔隙度体积分数和温度变化对屈曲前后固有频率以及临界温度的影响。图 9给出了在不同的孔隙度(α=0, 0.1, 0.2)下，频率和温度的变化关系。在发生屈曲前，无量纲频率随着温度的升高而降低，直至接近零。接近零点的点为该系统发生屈曲的临界点。在临界温度后，频率随着温度的增加而增加。从图 9中还可以看出，临界温度随着空隙度的增加而升高。从式(4)可知，孔隙度的改变会改变材料的杨氏模量、密度和热膨胀系数等，进而影响系统的刚度。

图 10展示了不同的功能梯度体积分数下，温度的变化对频率的影响。在没有达到临界温度时，无量纲频率随着温度的增加而减小，在达到临界温度之后，这种行为则完全相反。由式(3)可知，当梯度指数为零时，输流管材质为陶瓷单一材料。而当梯度指数为∞时，输流管的材质为金属单一材料。从图 9可以看出，临界温度随着梯度指数的增加而减小，这是因为当梯度指数增加时，材料中的金属含量增加，而金属更易受温度的影响，陶瓷的耐热性能更好。

图 11显示了不同的长径比下，温度变化对频率的影响。长径比小会使频率下降得更加缓慢，而且频率变化的曲线向高温方向移动。长径比大会使得输流管的结构刚度降低。故长径比越大，发生屈曲的临界值就越小。系统就越容易发生屈曲，在实际工程中，通过调整长径比来改变系统的稳定性也是较为常见的方法。

6 结论

(1) 对于较低的流速和初始轴向力，FGM管道保持稳定，静态分叉变形率为零。当流速和初始轴向力进一步升高超过各自的临界值时，系统发生叉式分叉并进入不稳定的屈曲状态。

(2) 在一阶屈曲配置附近，横向位移与幂律指数成正比。相反，随着幂律指数的增加，控制参数的临界值降低。而孔隙度体积分数与横向位移成反比。与幂律指数不同的是，随着孔隙度体积分数的增加，控制参数的临界值有明显提高。

(3) 屈曲前后某一阶模态的固有频率是近似相等的。而控制参数对后屈曲的高阶频率几乎没有影响。

(4) 初始轴向力、温度、幂律指数、孔隙度体积分数和长径比对第一阶屈曲位形的动态特性有显著的影响。