基于改進型GDQ法FGM納米梁的熱-機耦合振動及屈曲特性分析

周鳳璽， 蒲 育，2

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

近些年，隨著微納米測試與微機械加工技術的迅猛發(fā)展，許多材料力學試驗已表明:微/納米結構的靜動態(tài)力學特性均表現(xiàn)出顯著的尺度效應。Eringen[1]非局部彈性理論含有1個內稟尺度效應參數，很好地解釋了尺度效應，從而基于該理論對微/納結構力學特性分析及應用日益廣泛[2]。功能梯度材料（functionally graded material，F(xiàn)GM）微/納米梁模型以其簡單而高效的結構形式，F(xiàn)GM優(yōu)異的熱穩(wěn)定性、可設計性與智能可控性等，目前常見于微機電系統(tǒng)（micro-electromechanical-systems，MEMS）中的核心組件，諸如微執(zhí)行器、微傳感器、微諧振器等。分析此類微/納結構的力學特性關乎MEMS核心組件的安全與設計、功能與優(yōu)化、智能與控制等，這也是復合材料納米力學優(yōu)先發(fā)展的前沿課題，因而備受學者們關注。

基于Eringen非局部彈性理論，Thai[3]采用一種拋物線型剪切梁理論，應用Navier法分析了簡支納米梁的彎曲、屈曲及自由振動。Ke等[4]基于Timoshenke梁理論（first-order shear deformation beam theory，F(xiàn)SBT），采用微分求積法（differential quadrature method，DQM）研究了壓電納米梁的熱-機-電耦合振動特性。文獻[5-6]考慮了各向同性材料，均采用Euler-Bernoulli梁理論（classical beam theory，CBT），劉燦昌等研究了納米梁的非線性振動問題，張大鵬等利用傳遞函數法分析了黏彈性地基上納米梁的振動特性。Rahmani等[7]基于FSBT并采用Navier法，首次研究了FGM納米簡支梁的自由振動，其結果表明:尺度效應對FGM納米梁的固有頻率影響顯著，納米尺度因子越大，則頻率越低。Ebrahimi等[8-11]考慮了不同溫度分布類型下材料物性隨溫度變化的相關性，采用CBT，F(xiàn)SBT并應用Navier法及微分變換法（differential transform method，DTM），先后研究了FGM納米梁的撓曲振動、熱-機耦合振動、熱屈曲問題。最近，文獻[12-14]均采用CBT，Ebrahimi等[12]應用DTM數值研究了FGM納米梁的濕-熱-機耦合振動特性；Hosseini等[13]應用Navier法獲得了FGM納米簡支曲梁自由振動響應的精確解；Nejad等[14]應用廣義微分求積法（generalized differential quadrature，GDQ）數值研究了雙向FGM納米梁的屈曲特性，但未考慮熱環(huán)境工況。Ansari等[15]基于FSBT考慮了幾何非線性，將Newton-Raphson迭代法和GDQ法相結合，數值研究了壓電FGM納米梁在熱-機-電耦合作用下的非線性自由振動問題。Ebrahimi等[16]基于Reddy三階剪切梁理論（third-order shear deformation beam theory，TBT），考慮了微孔隙的影響，應用Navier法研究了彈性地基上壓磁/電FGM納米簡支梁的磁-電耦合振動特性。

綜上所述，目前針對熱-機耦合作用下FGM納米梁振動與屈曲特性的研究仍十分有限，且這兩類問題都需要求解復雜耦合微分方程組的特征值問題，獲得解析解十分困難，故大多只限于簡支梁而采用Navier法獲得精確解。其次，為了將研究問題簡化分析，大多學者采用了CBT模型，由于該梁理論忽略了剪切作用的影響，故而由此預測的FGM納米梁靜動態(tài)響應其精確性仍有待商榷。此外，已有文獻對熱-機耦合作用下FGM納米梁屈曲與振動這兩類靜動態(tài)力學行為耦聯(lián)性的揭示目前十分少見，仍需進一步研究。

本文采用n階廣義梁理論（generalized beam theory，GBT）[17]，基于非局部線彈性理論，在Hamilton體系下統(tǒng)一建立熱-機耦合FGM納米梁的振動與屈曲控制方程。應用前期研究基礎之上提出的一種改進型廣義微分求積法（modified generalized differential quadrature，MGDQ）[18-19]求解動態(tài)響應，通過屈曲與振動這兩類力學行為之間的二元耦聯(lián)性，編寫相應的循環(huán)子程序用來求解屈曲靜態(tài)響應。通過算例刻畫了多參數（因素）共同對FGM納米梁振動及屈曲特性的影響，分析了其影響作用機理。

1 控制微分方程

1.1 材料屬性

如圖1所示，一長寬高分別為L×b×h的FGM納米梁，上表面為純陶瓷，下表面為純金屬。梁受初始軸向機械載荷為N（假設壓力為正），考慮溫度T沿梁的厚度方向溫態(tài)分布，材料的物性沿梁厚按Voigt混合冪律呈梯度分布且依賴于溫度變化。彈性模E，熱膨脹系數α，熱傳導率κ，泊松比υ，密度ρ等物性參數均是坐標z和溫度T的函數，可用統(tǒng)一式表示為

圖1 FGM納米梁的幾何尺寸Fig.1 Geometry of a FGM nanobeam

式中:Pm（T）和Pc（T）分別為金屬與陶瓷材料與溫度相關的某一物性參數；p為FGM納米梁材料組分的梯度指標。顯然，當p=0時，F(xiàn)GM退化為純陶瓷材料；當p=∞時，退化為純金屬材料。Pm（T）和Pc（T）隨溫度T的變化可統(tǒng)一表述為

式中，Pi（i=-1，0，1，2，3）為與溫度相關的材料系數。表1給出了金屬（SUS 304）和陶瓷（Al2O3）這兩種材料與溫度相關的材料系數[20]。

表1 金屬（SUS 304）和陶瓷（Al2 O3）兩種材料隨溫度變化的物性系數Tab.1 Temperature-dependent coefficients for metal（SUS304）and ceramic（Al2 O3）

1.2 升溫類型

本文考慮以下3種升溫類型:

類型1——均勻升溫（uniform temperature rise，UTR）

式中:T0為無應力狀態(tài)時的參考溫度，本文取T0=300 K；ΔT為升溫值。

類型2——線性升溫（linear temperature rise，LTR）

式中，ΔT=Tc-Tm。

類型3——非線性升溫（non-linear temperature rise，NLTR）

1.3 運動方程

梁內任意一點在t時刻沿x，y，z方向的位移分量可描述為

式中:u（x，t），wb（x，t），ws（x，t）分別為梁軸線上一點的軸向位移、彎曲變形項撓度及剪切變形項撓度；f（z）=h（2z/h）n/（2n）為n階廣義梁理論對應的橫向切應力形函數。顯然，GBT可退化為3種著名的梁理論:當n=1時，退化為CBT；當n=3時，退化為TBT；當n=∞時，退化為FSBT。

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

式中，g（z）=1-f′（z）為橫向切應變形函數。

由Eringen非局部線彈性理論，本構方程為

式中:e0為材料的非局部常數；a為內部特征長度；μ=（e0a）2定義為尺度效應非局部參數；L為外部特征長度；G（z，T）為剪切模量；?2為一維Laplace算子。

結構動能的變分

結構應變能的變分

初始軸向機械載和熱載荷做功的變分

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

將式（7）～式（11）代入式（12）并作一些變分和積分運算，化簡后可得運動方程及邊界條件

顯然，當式（13）中與時間t無關時，則問題退化為FGM納米梁的熱-力耦合屈曲方程。此外，需要指出，嚴格意義上來講，Hamilton原理導出的邊界條件為尺度效應相關的非局部邊界條件，但基于Saint-Venant原理，可略去梁左右兩端小邊界處尺度效應相關的高階微量，進而經簡化并化簡后，考慮以下3種梁邊界:

（i）左端固支-右端固支（C-C）

（ii）左端簡支-右端簡支（S-S）

（iii）左端固支-右端簡支（C-S）

1.4 振動控制微分方程

自由振動為簡諧運動，位移分量可設為

式中:U（x），Wb（x），Ws（x）為振型函數；ω為振動頻率；i為虛數單位。

將式（17）代入運動方程式（13）可得FGM納米梁熱-機耦合振動的控制微分方程

2 MGDQ法與特征值問題

首先，GDQ實施離散化前，不失一般性，各參數采用如下的無量綱化

式中:λ為梁的跨厚比；I為慣性矩；Np和Ncr分別為無量綱初始軸向機械載荷和臨界屈曲載荷；Ω為無量綱頻率。

2.1 引入邊界控制參數

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

其次，可定義ξ=0及ξ=1處任意未知函數0階導數的權系數

據此定義，納米梁的3種邊界條件方程式（14）～式（16）可由GDQ統(tǒng)一離散化為

式中:n0=0或1；n1=0或1；n0和n1為邊界條件控制參數；n0=0，n1=0為C-C梁；n0=0，n1=1為CS梁；n0=1，n1=1為S-S梁。

2.2 MGDQ法的無量綱離散方程

考慮到應用傳統(tǒng)的GDQ離散控制方程與邊界方程后，方程總數和節(jié)點振型位移未知量的數目不相等，不能直接將振動問題轉化為特征值問題求解其動態(tài)響應。因此本文在前期研究基礎之上，采用一種改進型GDQ法進行求解。即從離散邊界條件方程式（22）出發(fā)，消去邊界點處的6個節(jié)點位移，同時解出內點處的這4個節(jié)點位移，使控制方程離散數目與振型位移獨立未知量的數目相等。其改進型權系數矩陣具體推導可參考蒲育等的研究。

將無量綱式（19）代入控制方程式（18）并應用MGDQ法離散后可得3種邊界下FGM納米梁熱-機耦合振動的統(tǒng)一化無量綱離散線性方程組

2.3 振動與屈曲兩類問題的特征值

線性代數方程式（23）可用分塊矩陣表示為

式中:[K]與[M]分別為彈性剛度矩陣與質量矩陣，且均為（3N-10）階方陣；[KN]與[KT]分別為軸向機械載荷及熱載荷作用對應引起的子剛度矩陣；{X}為獨立節(jié)點振型位移列向量，其表達式為

顯然，求解線性方程式（24）的特征值問題可獲得FGM納米梁熱-機耦合振動的頻率和相應的振型。通過振動與屈曲的二元耦聯(lián)性，編寫MATLAB循環(huán)子程序可求解得出臨界屈曲載荷、熱屈曲及熱-力耦合屈曲臨界升溫值、臨界跨厚比等重要參數指標。

3 算例分析與討論

3.1 GBT階數n值的探討及結果的有效性

本文取離散節(jié)點個數N=17，納米梁的長寬尺寸分別為L=104nm，b=103nm，非局部參數采用 μ×10-6（nm）2的形式表示量級和量綱。本節(jié)主要探討不同梁理論對FGM納米梁預測頻率的影響并對數值結果的有效性進行驗證。

表2給出了均勻升溫（UTR）熱載荷作用下C-C，C-S及S-S這3種邊界FGM納米梁（λ=20，μ=1，ΔT=20 K）的前三階無量綱頻率Ω并與文獻[9]的結果進行了比較。本文采用TBT（n=3），而文獻[9]為了簡化分析采用了CBT（n=1）。不難看出:兩者預測的基頻結果比較吻合。注意到算例中跨厚比λ=20為細長梁，故而兩種梁理論預測的基頻結果比較接近。但隨著振型階數的增加，對三階無量綱頻率Ω3值而言，文獻[9]的結果明顯偏高。這是由于CBT忽略了剪切變形作用的影響而高估了梁的整體剛度。換言之，即使對于納米尺度細長梁而言，CBT比TBT預測的高階頻率值也明顯偏高。此外，由表2可見:約束越強，相應的頻率值越大；梯度指標p越大，頻率值越小。

表2 不同邊界FGM納米梁均勻升溫下的前三階無量綱頻率（λ=20，μ=1，ΔT=20 K）Tab.2 First three dimensionless frequencies of FGM nanobeams under the UTR for different boundary conditions（λ=20，μ=1，ΔT=20 K）

圖2主要用于驗證并探討GBT階數n的理想取值?？紤]了UTR熱載荷作用，圖2（a）～圖2（c）分別給出了FGM納米固支梁（p=1，μ=1，λ=20，ΔT=40 K）的前三階無量綱頻率Ω隨GBT階數n變化的關系曲線:總的來看，各階頻率Ω都隨n的變化波動很小，n取偶數略微高于n取奇數（除n=1外）時預測的頻率，隨著n的增加，各階Ω值趨于穩(wěn)定，將趨于FSBT預測的頻率值。n=2時各Ω值為最大，n=3時（TBT）各 Ω值為最小，n=1（CBT）預測的頻率值介于期間。需要指出，取n=1實施計算過程中，結果顯示剛度矩陣奇異，此時可取n≈1代替（如取n=1.001），將此作為退化的CBT。與文獻[9]CBT預測的相應頻率值Ω1=14.513 1，Ω2=35.406 4，Ω3=59.492 1進行比較:由圖2可見，本文退化的CBT預測各階頻率值均偏低，且更為接近高階剪切梁理論TBT的預測值。

圖2 UTR升溫下FGM納米固支梁的前三階無量綱頻率Ω與GBT階數n的關系曲線Fig.2 First three dimensionless frequencies of a C-C FGM nanobeam under the UTR versus orders of n for the GBT

綜上所述，本文所采用的MGDQ法切實可行，行之有效。另一方面，即使對于納米尺度細長梁而言，CBT比TBT預測的高階頻率值也明顯偏高。因此，為使在高頻振動中獲得精確的輸出響應、同時也便于工程中使用，建議GBT的階數取n≥3的奇數為宜。

3.2 FGM納米梁的耦合振動及耦合屈曲特性

本節(jié)采用TBT，主要探討FGM納米梁的熱-機耦合振動及屈曲特性。圖3（a）～圖3（c）分別刻畫了C-C，C-S，S-S這3種邊界FGM納米梁（p=1，λ=20）的Ω1-Np關系曲線:無量綱基頻Ω1都隨無量綱軸向機械壓力Np的增加而單調減小，直至達到各自臨界屈曲載荷Np=Ncr時而發(fā)生失穩(wěn)，此時Ω1=0。換言之，由Ω1-Np曲線與橫軸交點的坐標值可確定Ncr值。這是由于軸向機械壓力削弱了梁結構的整體剛度，當Np=Ncr時，F(xiàn)GM納米梁達到屈曲狀態(tài)而發(fā)生臨界失穩(wěn)，由于臨界失穩(wěn)為靜平衡狀態(tài)，故而Ω1=0。這也反映了FGM納米梁屈曲和振動這兩類靜動態(tài)力學行為之間的轉化和二元耦聯(lián)性。本文基于該二元耦聯(lián)性，通過編寫相應的MATLAB循環(huán)子程序用以求解FGM納米梁屈曲靜態(tài)響應，這樣可避免屈曲問題求解的二次解耦，二次提高計算效率。此外，由圖3可見:對μ=0，μ=1，μ=2，μ=3，μ=4這5種取值而言，非局部參數μ值越大，基頻Ω1和臨界載荷Ncr越小。這是由于μ值越大，則尺度效應越顯著，納米梁的整體剛度將明顯弱化的緣故。比較圖3（a）～圖3（c）可見:取相同的參數值，約束越強，頻率和臨界載荷值越大，則納米梁越不易失穩(wěn)。這與宏觀尺度（μ=0）FGM梁類似，邊界約束強，梁抵抗變形的能力就強，進而抗失穩(wěn)能力也就提高了。

圖3 不同邊界下FGM納米梁的無量綱基頻Ω1與無量綱初始軸向壓力N p的關系曲線Fig.3 Curves of dimensionless fundamental frequencyΩ1 of FGM nanobeams versus dimensionless initial axial compressive force N p for different boundary conditions

為了便于下文分析不同升溫類型對FGM納米梁頻率和臨界載荷的影響及作用機理，考慮了UTR，LTR，NLTR這3種升溫類型，圖4給出了FGM納米梁（p=1）的拉/壓彈性模量系數S1隨升溫ΔT的變化關系:各S1值均隨ΔT的增加而單調減小。取相同的ΔT值，UTR使S1減小最為顯著，LTR次之，NLTR最為不明顯。其余彈性模量系數Si與之類似，這里不再一一驗證。換言之，由于本文考慮了材料物性與溫度的相關性，對這3種溫度分布升溫類型而言，梁結構的彈性剛度均隨升溫ΔT的增加而降低了，且UTR使其減小最為明顯，NLTR使其減小最為緩和，進而UTR使頻率和臨界載荷減小最為顯著。

圖4 不同升溫類型下彈性系數S1與升溫ΔT關系曲線Fig.4 Curves of elastic coefficient S1 versus different types of temperature riseΔT

圖5（a）、圖5（b）分別給出了UTR和LTR升溫下3種邊界FGM納米梁（p=1，μ=1，λ=20）的Ω1-ΔT和Ncr-ΔT關系曲線:各基頻Ω1和臨界載荷Ncr曲線值均隨升溫ΔT的增加而單調減小，直至達到各自的臨界升溫值ΔTcr時發(fā)生熱屈曲而失穩(wěn)，此時Ω1=0，Ncr=0，兩類曲線與橫軸交點的坐標值幾乎相同，該值即為相應的ΔTcr。這是由于隨著ΔT的增加，各彈性系數減小了，進而結構的彈性剛度降低了；另一方面，增大的熱軸力為壓軸力，這也愈加削弱了梁結構的整體剛度，當ΔT=ΔTcr達到熱屈曲狀態(tài)時，臨界失穩(wěn)載荷完全由熱軸力提供，不需要額外的軸向機械壓力，故而Ncr=0，且熱屈曲為臨界平衡狀態(tài)，從而Ω1=0。這也再次表明了FGM納米梁兩類靜動態(tài)力學行為之間的二元耦聯(lián)性。此外，取相同的ΔT值時，同一邊界下UTR比LTR使Ω1和Ncr值減小更為明顯，進而UTR更易使FGM納米梁達到熱屈曲而失穩(wěn)；同一升溫類型下，邊界約束較強，則頻率、臨界載荷及臨界升溫值均越大。其作用影響機理前文已有分析，這里不再贅述。

圖5 FGM納米梁的無量綱基頻Ω1及無量綱臨界屈曲載荷N cr隨升溫ΔT變化關系曲線Fig.5 Curves of dimensionless fundamental frequencyΩ1 and dimensionless critical buckling load N cr of FGM nanobeams versus temperature riseΔT

圖6考慮了Np=-5，Np=-2，Np=0，Np=2，Np=5這5種取值，刻畫了FGM納米梁（p=1，μ=1，λ=20）在熱-機耦合作用下FGM C-S納米梁的Ω1-ΔT關系曲線:各Ω1曲線值都隨線性升溫ΔT的增加而單調減小，當ΔT=ΔTcr時，納米梁發(fā)生熱-機耦合屈曲而失穩(wěn)。此外，取相同的ΔT值，由Np=2和Np=5對應的這2條基頻曲線來看:初始軸向機械壓力越大，則Ω1和ΔTcr值越?。挥蒒p=-2和Np=-5對應的2條基頻Ω1曲線來看:初始軸向機械拉力越大，則Ω1和ΔTcr值也越大。這是由于初始軸向機械拉力提高了納米梁抗彎曲變形的能力。

圖6 C-S邊界FGM納米梁的無量綱基頻Ω1隨線性升溫ΔT變化關系曲線Fig.6 Curves of dimensionless fundamental frequencyΩ1 of C-SFGM nanobeam versus linear temperature riseΔT

圖7反映了NLTR升溫ΔT=200 K及5種非局部參數μ=0，μ=1，μ=2，μ=3，μ=4取值時FGM納米簡支梁（p=1，Np=5）熱-機耦合振動的Ω1-λ關系曲線:隨著跨厚比λ的增大，各基頻曲線Ω1值均先增大，隨后單調減小，達到各自臨界跨厚比λcr時，發(fā)生熱-機耦合屈曲而失穩(wěn)，此時Ω1=0。這是由于熱載荷、機械載荷及跨厚比多因素共同作用引起的，此時FGM納米梁的力學行為表現(xiàn)出復雜的耦合特性。此外，取相同的跨厚比λ，非局部參數μ值越大，基頻Ω1越小，納米梁達到熱-機耦合屈曲時，λcr值也越小。即尺度效應越顯著，頻率和臨界跨厚比越小，F(xiàn)GM納米梁也更易失穩(wěn)。

圖7 FGM納米簡支梁熱-機耦合振動的無量綱基頻Ω1隨跨厚比λ變化關系曲線Fig.7 Curves of dimensionless fundamental frequency Ω1 of S-SFGM nanobeam subjected to thermal-mechanical loads versus slenderness ratiosλ

圖8反映了3種邊界FGM納米梁（μ=1，λ=20，Np=5，ΔT=100 K）熱-機耦合振動的Ω1-p關系曲線:各基頻曲線Ω1值先隨材料組分梯度指標p的增加而顯著減小，隨后減小趨于緩慢。這是由于p增加的初始階段，陶瓷成份Al2O3驟減，F(xiàn)GM納米梁的整體剛度迅速降低而導致的。

圖8 熱-機耦合作用下FGM納米梁的無量綱基頻Ω1隨梯度指標p變化關系曲線Fig.8 Curves of dimensionless fundamentalΩ1 of FGM nanobeams subjected to thermal-mechanical loads versus material graded index p

4 結 論

通過引入邊界條件控制參數，采用MGDQ法實施了3種邊界FGM納米梁動態(tài)響應求解的MATLAB統(tǒng)一化編程，初次優(yōu)化了數值分析方法?；陟o動態(tài)響應之間的二元耦聯(lián)性實施解耦，通過循環(huán)子程序獲得屈曲靜態(tài)響應，二次優(yōu)化了MGDQ法并提高了計算效率，為FGM微/納結構的力學行為研究提供了一種切實可行、行之有效的分析方法。且研究結果表明:

（1）即使對于納米尺度細長梁而言，隨著階次的增加，CBT仍明顯高估納米梁的高階頻率。而GBT具有重要的理論和工程應用價值，建議推廣使用。

（2）FGM納米梁的熱-機耦合振動頻率、臨界載荷、臨界升溫值均隨尺度效應非局部參數、軸向機械壓力、升溫、梯度指標的增加而減小。溫度分布不同，其影響有顯著差異。即UTR使FGM納米梁頻率減小最快，更易使梁失穩(wěn)、LTR次之、NLTR則最為緩和。

（3）隨著跨厚比的增加，F(xiàn)GM納米梁熱-機耦合振動的基頻先增加，隨后減小，達到臨界跨厚比時發(fā)生熱-機耦合屈曲而失穩(wěn)，此時基頻減小為零；微/納尺度效應越明顯，臨界跨厚比則越小，F(xiàn)GM納米梁此時更易失穩(wěn)。