蒲 育，周鳳璽，任永忠，劉 君

(1. 蘭州工業(yè)學院土木工程學院，蘭州 730050；2. 蘭州理工大學土木工程學院，蘭州 730050)

功能梯度材料(FGM)梁作為當前現(xiàn)代工程領域中廣泛使用的結構形式之一，其穩(wěn)定性和振動特性關乎結構的安全設計和功能設計，一直備受人們重視。隨著近些年FGM 梁廣泛服役于高溫環(huán)境或特定復雜工況，F(xiàn)GM 梁的熱-力耦合振動特性及穩(wěn)定性研究顯得尤為必要。

例如，Pradhan 和Murmu[1]基于Euler 梁理論(CBT)，應用微分求積法(DQM)數(shù)值研究了變彈性地基FGM 夾層梁的振動特性。趙鳳群和王忠民[2]采用小波微分求積法(WDQ)分析了熱載荷和切向隨從力耦合下FGM CBT 梁的穩(wěn)定性及振動特性。考慮沿梁厚度方向非線性升溫，Mahi 等[3]獲得了材料物性關于幾何中面對稱FGM 梁自由振動響應的精確解。文獻[4]基于Reddy 三階剪切梁理論(TBT)，應用Ritz 法探討了均勻升溫下FGM 梁的振動與屈曲特性。之后，文獻[5]采用微分變換法(DTM)數(shù)值分析了兩端彈性約束下多孔FGM CBT 梁的線性及非線性振動特性。Trinh 等[6]基于TBT 并采用狀態(tài)空間法(SSM)，分析了熱-力耦合下無孔FGM 梁的振動及屈曲特性。文獻[7 ? 10]則從材料加工缺陷或特殊功能要求角度出發(fā)，僅考慮簡支邊界梁，應用Navier 法研究了含孔隙FGM梁的耦合振動及屈曲特性。Ebrahimi 和Jafari[7 ? 8]先后采用TBT 及雙曲型(HBT)兩種高階剪切梁理論，研究了多孔FGM 梁的熱-力耦合振動特性。蘇盛開和黃懷緯[9]采用CBT 及正弦型剪切梁理論(SBT)，分析了多孔FGM 梁的熱-力耦合屈曲行為。最近，筆者研究團隊[10 ? 11]擴展并提出一種廣義梁理論(GBT)，分別選用兩組不同的位移量用來描述位移場，探討了濕-熱-機耦合作用下多孔FGM 梁的振動及屈曲特性。此外，文獻[12 ? 19]則考慮了幾何非線性，采用不同的梁理論或應用不同的分析方法研究了FGM 梁的熱后屈曲及非線性振動特性。

綜上所述，目前針對熱-力耦合多孔FGM 梁振動及屈曲特性的研究十分有限，且兩者都需要求解微分方程組的特征值問題，解耦不易，獲得解析解十分困難，故大多只限于簡支邊界梁。其次，臨界載荷與頻率作為刻畫FGM 梁靜動態(tài)力學行為的重要指標，它們對FGM 梁的安全設計和功能設計具有重要的指導意義，但目前針對這兩種力學行為之間耦聯(lián)性的研究也十分少見。最后，考慮多場耦合及多因素影響下，其耦合影響作用機理尚不十分清楚，仍需進一步揭示。

本文采用GBT，基于兩類問題的二元耦聯(lián)性，在Hamilton 體系下統(tǒng)一建立描述多孔FGM梁熱-力耦合振動與耦合屈曲問題力學模型的控制方程。在前期研究基礎之上，采用MGDQ 法[20]獲得數(shù)值解。通過算例，著重探討多參數(shù)對多孔FGM 梁耦合振動及耦合屈曲特性的影響。揭示多孔FGM 梁屈曲與振動這兩種靜動態(tài)力學行為之間的二元耦聯(lián)性，重點分析多參數(shù)對頻率、臨界載荷、臨界升溫值影響的耦合作用機理。

1 控制微分方程

1.1 材料屬性描述

如圖1 所示，xoy 面為梁的幾何中面，現(xiàn)考慮一長、寬、高分別為L×b×h 且含均勻孔隙的FGM梁，上表面為陶瓷，下表面為金屬。梁兩端受初始軸向機械力為N(假設壓力為正)，溫度T 沿梁厚度方向穩(wěn)態(tài)分布，材料的物性沿梁厚呈梯度連續(xù)分布且與溫度相關，采用含孔隙率β 修正后的Voigt 混合冪律模型，物性參數(shù)(如彈性模量E、泊松比ν、熱膨脹系數(shù) α、熱傳導率κ)與坐標z 及溫度T 可采用統(tǒng)一式表述為[10]：

圖1 多孔FGM 梁的幾何尺寸Fig.1 Geometry of a porous FGM beam

金屬與陶瓷這兩種材料的某一物性參數(shù) X隨溫度 T的變化可統(tǒng)一表述為[11]：

式中， Pi(i=?1,0,1,2,3)表示隨溫度變化的材料系數(shù)。表1 給出了金屬(SUS 304)和陶瓷(Si3N4)與溫度相關的材料系數(shù)[11]。

表1 金屬(SUS 304)和陶瓷(Si3N4)兩種材料隨溫度變化的物性系數(shù)[11]Table1 Temperature-dependent coefficients for metal (SUS 304) and ceramic (Si3N4)[11]

1.2 升溫類型

本文考慮以下三種升溫類型：

類型1 均勻升溫(UTR)：

式中， T0表示無應力狀態(tài)時的參考溫度，本文取T0=305 K， ?T表示升溫值。

類型2 線性升溫(LTR)：

式中，?T =Tc?Tm

類型3 非線性升溫(NLTR)：

1.3 多孔FGM 梁振動及屈曲的微分方程

位移場可描述為：

由小變形的幾何方程可得非零應變分量：

g(z)= f′(z)

式中，為橫向切應變形函數(shù)。

本構方程為：

G(z,T)=E(z,T)/2[1+ν(z,T)]

式中，為剪切模量。

動能的變分為：

應變能的變分為：

內(nèi)力由式(11)和式(12)定義并可由位移分量表示為：

熱-力耦合載荷做功的變分為：

對系統(tǒng)應用Hamilton 原理：

將式(6)～式(13)代入式(14)并作一些變分和積分運算化簡后可得由位移分量表示的多孔FGM梁自由振動微分方程：

顯然，彈性系數(shù)與慣性系數(shù)項中存在拉/壓-剪切-彎曲的耦合。當式(15)中與時間t無關時，則退化為多孔FGM 梁熱-力耦合屈曲的控制方程。

由式(14)導出的梁自然邊界條件由內(nèi)力式(11)化簡后可得由位移分量表示的以下3 種邊界：

兩端固支(C-C)：

兩端簡支(S-S)：

左端固支-右端簡支(C-S)：

1.4 多孔FGM 梁振動及屈曲的控制方程

對于諧振動，位移分量可設為：

式中： ω為固有頻率；U(x)、Ψ(x)、W(x)為振型位移；i為虛數(shù)單位；e 為自然底數(shù)。

將式(19)代入式(15)可得多孔FGM 梁熱-力耦合振動及屈曲的統(tǒng)一控制微分方程：

2 MGDQ 法的實施及特征值問題

不失一般性，各參數(shù)無量綱化如下：

2.1 邊界控制參數(shù)

通過引入邊界控制參數(shù)以實施C-C、C-S、SS 這3 種邊界梁動態(tài)響應求解的MATLAB 統(tǒng)一化編程，進而可優(yōu)化GDQ 法這一數(shù)值方法。首先，基于GDQ 法[20]，位移振型函數(shù)在離散節(jié)點ξ=ξi處的k階導數(shù)可表示為：

其次，這里定義ξ=0 及ξ=1處任意未知函數(shù)0 階導數(shù)的權系數(shù)：

據(jù)此定義，3 種梁邊界對應的式(16)～式(18)可由GDQ 法統(tǒng)一離散為：

式中：n0=0 或1，n1=0 或1，其為引入的邊界控制參數(shù)。n0=0，n1=0 表示C-C 梁；n1=1，n1=1表示S-S 梁；n0=0，n1=1 表示C-S 梁。

2.2 推導MGDQ 法的權系數(shù)

應用傳統(tǒng)GDQ 離散控制方程與邊界方程后，方程總數(shù)多于振型位移未知量的數(shù)目，不能直接轉化為特征值問題求解。因此本文在前期研究基礎之上，采用一種改進型GDQ 法(MGDQ)獲得數(shù)值解[20]。由式(22)及式(24)推導易得：

式中：

2.3 振動與屈曲問題的MGDQ 法離散化

將式(21)代入控制方程式(20)并應用MGDQ法可統(tǒng)一獲得C-C、C-S、S-S 這3 種邊界多孔FGM梁熱-力耦合振動及屈曲問題的無量綱離散方程：

2.4 振動與屈曲問題的特征值

通過應用MGDQ 法離散化后，方程組式(29)可化為標準特征值：

式中：[K]為彈性剛度矩陣；[M]為質量矩陣；[KT]與[KN]分別為熱載荷及軸向機械力作用對應的剛度矩陣；{ X}為節(jié)點位移列向量：

求解式(30)的特征值問題可獲得頻率和相應的振型；通過振動與屈曲的二元耦聯(lián)性，編寫MATLAB 循環(huán)子程序可求解得出臨界載荷、臨界升溫值及臨界跨厚比等重要參數(shù)指標。

3 算例與討論

3.1 GBT 的有效性及階數(shù)n 取值的探討

首先，將問題退化，令β=0, 作為本文算例數(shù)值結果的有效性驗證：考慮了LTR 及NLTR 這2 種升溫類型，表2 給出了C-C、C-S 及S-S 這3 種邊界下無孔FGM 梁(p=1, λ=20, ΔT=20 K, n=3)的無量綱基頻 ?1。本文采用MGDQ 法得出的結果與文獻[6]采用SSM、文獻[21]采用DTM 分析獲得的結果基本吻合。

表2 FGM 梁無量綱基頻的結果比較Table2 Comparison of dimensionless fundamental frequencies of FGM beams.

其次，圖2 反映了NLTR 下FGM 固支長梁(p=1, λ=20, ΔT=20 K)的基頻 ?1隨GBT 階數(shù) n變化的關系曲線：總的來看， ?1隨階數(shù)n的變化波動極小，換言之，不同梁理論對FGM 長梁基頻的預測值相差無幾；取n=2,4,6···偶數(shù)時略微高于取n=3,5,7···奇數(shù)時預測的基頻；n =2 與n=3(TBT)預測的基頻值分別為最大和最小，而CBT(n =1)預測值介于其間，極差僅為0.06666，CBT 相對TBT 的預測誤差為0.3363%；隨著 n 的增加，基頻 ?1趨于FSBT 的預測值。

圖3 給出了NLTR 下FGM 固支短梁(p=1, λ=5,ΔT=20 K)的無量綱臨界載荷Ncr與GBT 階數(shù) n之間的關系曲線：CBT 相對于TBT 明顯高估臨界載荷，相對誤差可達7.05967%，其不容小覷。這是由于剪切作用對短梁影響明顯，而CBT 忽略了剪切變形作用的影響，故而明顯高估了臨界載荷。綜上所述，為提高梁理論預測值的精確性并便于在工程中使用，GBT 階數(shù)宜取n≥3的奇數(shù)最為理想。

圖2 NLTR 升溫下FGM 固支梁的無量綱基頻 ?1與GBT 階數(shù)n 的關系曲線Fig.2 Curve of non-dimensional fundamental frequency ?1 of FGM C-C beam under NLTR versus order n of GBT

圖3 NLTR 升溫下FGM 固支梁的無量綱臨界載荷Ncr 與GBT 階數(shù)n 的關系曲線Fig.3 Curve of non-dimensional critical load Ncrof FGM CC beam under NLTR versus order n of GBT

3.2 影響因素分析及二元耦聯(lián)性探究

本節(jié)重點探討熱-力耦合作用下多因素(參數(shù))對多孔FGM 梁振動與屈曲特性的影響以及這兩種力學行為之間的二元耦聯(lián)性(取n=3)。

圖4 無量綱基頻Ω1 與無量綱軸向載荷 關系曲線Fig.4 Curves of dimensionless fundamental frequency Ω1 versus non-dimensional axial load

為了便于下文分析，圖5 反映了FGM 梁(p=1,β=0.0)在3 種不同溫度分布下彈性模量系數(shù)S1隨升溫ΔT 的變化：3 種溫度分布對應的S1值均隨升溫ΔT 的增加而明顯減?。簧郎卅 相同時，UTR分布下S1值減小最為顯著，LTR 分布次之，NLTR分布最為緩和；升溫值ΔT 較小時(如ΔT ≤ 150 K)，LTR 和NLTR 對應的S1值十分接近；其他各彈性系數(shù)Si-ΔT 與S1-ΔT 具有相類似的變化關系，即各彈性系數(shù)Si隨ΔT 的增加而減小，進而梁的彈性剛度也就降低了。另一方面，ΔT 增加，增大的熱軸力弱化了梁的整體剛度。綜上分析可知：FGM 梁的整體剛度隨升溫ΔT 的增加而減小了。同時注意到各慣性系數(shù)Ii并不隨升溫而改變，進而梁的整體等效質量也不隨升溫而改變，因此，無孔FGM 梁的 ?1和Ncr都隨ΔT 的增加而減??；升溫相同時，NLTR 分布使梁的頻率、屈曲臨界載荷及臨界升溫值減小最為緩和，LTR 分布次之，UTR分布使其減小則最為明顯；ΔT 較小時，如ΔT=20 K，LTR 與NLTR 分布下同一邊界FGM 梁的基頻值相差很小，表2 的數(shù)值結果則正好驗證了該點。

圖5 不同溫度分布下彈性系數(shù)S1 與升溫ΔT 關系曲線Fig.5 Curves of elastic coefficient S1 versus temperature rise ΔT under different temperature distributions

圖6 無量綱基頻 ?1及臨界載荷Ncr 與升溫ΔT 關系曲線Fig.6 Curves of dimensionless fundamental frequency ?1 and critical load Ncr versus temperature rise ΔT

圖6 主要用于FGM 梁在熱環(huán)境中兩種力學行為之間的二元耦聯(lián)性探討。圖6(a)刻畫了熱載荷作用下無孔FGM 梁(p=1, λ=20, β=0)自由振動的?1-ΔT 曲線，圖6(b)則給出了該梁熱-力耦合屈曲的Ncr-ΔT 曲線。首先，本算例中，邊界條件、溫度分布、升溫對 ?1及Ncr的影響上文已有詳述，這里不再贅述。其次，由 ?1-ΔT 及Ncr-ΔT 兩類曲線不難看出： ?1及Ncr都隨升溫ΔT 的增加而單調減小，直至達到臨界升溫值ΔTcr而發(fā)生熱屈曲，此時， ?1=0，Ncr=0。計算結果表明：相同參數(shù)下兩類曲線與橫軸ΔT 交于同一點，即該點的橫坐標值為熱屈曲ΔTcr。這是因為對梁的熱-力耦合屈曲而言，隨著ΔT 的增加，熱軸力增大，當ΔT=ΔTcr時，熱軸力值剛好達到使梁發(fā)生屈曲的臨界壓力，臨界失穩(wěn)載荷此時完全由熱軸力提供，進而Ncr=0。另一方面，熱屈曲為臨界平衡狀態(tài)，故而此時 ?1=0。換言之，熱屈曲關聯(lián)了兩種力學行為，進而也表明了耦合條件下兩種力學行為之間的耦聯(lián)性。最后，作為ΔTcr值的有效性驗證，文獻[6]考慮了溫度LTR 分布：ΔTcr=1062.5 K (C-C)、612.5 K (C-S)、328.1 K (S-S)；本文結果：ΔTcr=1065.9 K (C-C)、617.6 K (C-S)、332.9 K (S-S)，兩者最大相對誤差不超過1.5%，吻合較好。特別注意到：C-S (LTR)比C-C (UTR)，S-S (LTR)比C-S(UTR)對應的ΔTcr值偏高，這是約束邊界和溫度分布共同作用影響的結果。

圖7 與圖8 考慮了LTR 升溫及4 種孔隙率β=0.00, 0.10, 0.15, 0.20, 分別刻畫了FGM C-S 梁(p=2,λ=20)的 ?1-ΔT 曲線及Ncr-ΔT 曲線：兩類曲線的?1與Ncr值都隨ΔT 的增加而單調減小，直至達到ΔTcr時發(fā)生熱屈曲使 ?1=Ncr=0，其ΔTcr=575.1 K(β=0.00)、628.8 K (β=0.10)、658.9 K (β=0.15)、693.6 K(β=0.20)，其極差可達118.5 K。

圖7 無量綱基頻 ?1與線性升溫ΔT 關系曲線Fig.7 Curves of dimensionless fundamental frequency ?1 versus linear temperature rise ΔT

當升溫值ΔT 較小時(如0≤ ΔT≤ 200 K)，從兩類曲線的疏密程度來看， ?1曲線比較密集，Ncr曲線明顯疏離，取相同的ΔT 值，β 越大，Ncr越小，梁越易達到熱-力耦合臨界失穩(wěn)，但 ?1減小甚微。從Ncr線的斜率來看，β 越小，Ncr值減小越快，進而隨著ΔT 的增加，兩種孔隙率β 取值的Ncr線將交于一點。

圖8 無量綱臨界載荷Ncr 與線性升溫ΔT 關系曲線Fig.8 Curves of dimensionless critical load Ncrversus linear temperature rise ΔT

當升溫ΔT 值較大時(如ΔT ≥ 500 K)，兩類曲線明顯都疏離，取相同的ΔT 值，β 越大， ?1與Ncr值均越大，孔隙率越小的梁反而更易達到熱-力耦合臨界失穩(wěn)。隨著升溫ΔT 繼續(xù)增加，熱軸力將繼續(xù)增大，直至達到各自的熱屈曲臨界平衡狀態(tài)。不難看出，β 越大，ΔTcr也越大，即孔隙越小，F(xiàn)GM 梁更易達到熱屈曲而失穩(wěn)。

綜上所述，在升溫的初始階段(ΔT 較小)，熱軸力本身較小，增大孔隙率，熱軸力進一步被削弱，它對梁的整體剛度的削弱影響更小了，另一方面，孔隙率增大同時弱化了梁的整體剛度和有效質量，兩者弱化程度相當時， ?1則隨β 的增大其影響十分有限。由于梁整體剛度的減少，這一階段Ncr隨β 的增大而減小。在升溫增加的后期(ΔT 較大)，熱軸力增大將明顯削弱梁的整體剛度，另一方面，孔隙率增大削弱了熱軸力，則熱軸力對梁整體剛度弱化的程度減小了，換言之，在該階段升溫ΔT 相同條件下，對較大與較小的兩種孔隙率梁而言，前者的整體剛度高于后者，進而使前者達到耦合屈曲所需額外的機械載荷值就更大一些，達到熱屈曲時的熱軸力也更大一些，故而β 越大，Ncr和ΔTcr都越大。同時注意到前者的等效質量低于后者，故而β 越大， ?1也就越大。

在同一坐標系下，圖9 刻畫了NLTR 升溫ΔT 時梯度指標p 對多孔FGM C-S 梁(β=0.1, λ=20,Nˉ=10) ?1與Ncr的影響：首先，3 種梯度指標值對應的兩類曲線 ?1和Ncr值都隨NLTR 升溫ΔT 的增加而單調減小，直至達到各自的臨界升溫值ΔTcr而耦合失穩(wěn)，且兩類對應曲線與橫軸交于同一點，該ΔTcr值分別為508.3 K (p=0)、285.4 K(p=1)、136.8 K (p=10)，本算例也表明了熱-力耦合屈曲與耦合振動這兩類力學行為之間的二元耦聯(lián)性。其次，取相同的升溫值(ΔT≤ 136.8 K)，p 值越大， ?1和Ncr值越小，即 ?1和Ncr值都隨p 的增大而減小。這是由于隨著p 的增大，陶瓷(Si3N4)成份減少而降低了梁的整體剛度，進而ΔTcr值也隨p 的增加而減小了。

圖9 熱-力耦合作用下多孔FGM C-S 梁無量綱基頻 ?1及臨界載荷Ncr 與非線性升溫ΔT 的關系曲線Fig.9 Curves of dimensionless fundamental frequency ?1 and buckling load Ncr of porous FGM C-S beam subjected to thermal-mechanical loads versus nonlinear temperature rise ΔT

圖10 考慮了NLTR 升溫ΔT=300 K 及3 種梯度指標p=0, 1, 10，在同一坐標系下，刻畫了熱-力耦合下多孔FGM C-S 梁(β=0.10, Nˉ=10) ?1與Ncr隨跨厚比λ 的影響關系曲線：隨著λ 的增加，兩類曲線 ?1和Ncr值均先增加，隨后單調減小，直至達到熱-力耦合屈曲相應的臨界跨厚比 λcr而失穩(wěn)，該值分別為：27.0 (p=0)、19.7 (p=1)、13.1 (p=10)。且計算結果表明：3 種梯度指標p 所對應曲線的駐點值 λ雖有所不同，但均 λ≤10，這通常定義為短梁的界限范圍，換言之，對長梁來說， ?1與Ncr都隨λ 的增大而單調減小，直至屈曲時 ?1=Ncr=0。這是由跨厚比、梯度指標、熱-力耦合效應共同影響導致的。此外，不難看出，p 值越大， λcr值越小，即 λcr也隨著陶瓷(Si3N4)成份的減少而迅速降低了。

圖10 熱-力耦合作用下多孔FGM C-S 梁無量綱基頻 ?1及臨界載荷Ncr 與跨厚比的關系曲線Fig.10 Curves of dimensionless fundamental frequency ?1 and buckling load Ncr of porous FGM C-S beam subjected to thermal-mechanical loads versus slenderness ratios λ

4 結論

首先，在Hamilton 體系下統(tǒng)一建立了描述多孔FGM 梁熱-力耦合振動與屈曲問題力學模型的控制方程。為了便于解耦并提高計算效率，通過引入邊界控制參數(shù)，實施了3 種典型邊界梁動態(tài)響應MGDQ 法求解的MATLAB 統(tǒng)一化編程?；趦深悊栴}的二元耦聯(lián)性解耦，通過循環(huán)子程序求解屈曲響應，避免了再次解耦，二次提高了計算效率。問題退化后的數(shù)值結果對比表明：該分析方法行之有效，具有精度高、高效且易實施、計算工作量小、通用性好等優(yōu)點。

其次，算例結果表明：熱-力耦合作用下多孔FGM 梁的頻率、臨界載荷、臨界升溫值、臨界跨厚比等重要指標受多因素共同的影響，其影響作用機理也表現(xiàn)出復雜的耦合特性。

(1) 只有機械載荷作用時，含孔比無孔更易使FGM 梁發(fā)生屈曲；考慮熱-力耦合作用時，在升溫值增加的前期，孔隙率越大，臨界載荷明顯越小，頻率則略有減少；在升溫值增加的后期，孔隙率越大，臨界載荷及頻率反而越大，達到熱屈曲時，臨界升溫值也越大。

(2) 熱-力耦合作用下多孔FGM 梁的頻率和臨界載荷隨跨厚比的增加先增大，隨后單調減小，直至達到耦合屈曲臨界跨厚比時而失穩(wěn)。

最后，GBT 具有重要的理論和工程應用價值，建議推廣并在工程中使用。