Mindlin 矩形微板的熱彈性阻尼解析解1)

馬航空 周晨陽 李世榮

(揚州大學建筑科學與工程學院，江蘇揚州 225127)

引言

微/納諧振器作為微/納電機系統(MEMS/NEMS)中的基本器件被廣泛用作頻率控制器件、邏輯開關、傳感器、驅動器、陀螺儀、能量收集器等，在生物醫學、自動控制、航空電子設備、國防科技等眾多領域有著十分廣泛的應用.為了設計和制造高性能的微/納電機系統，需要諧振器在工作時耗能越小越好，即具有更高的品質因子.然而，在微/納電機系統中不可避免地存在各種耗能機制.除了空氣阻尼、支撐阻尼、擠壓模阻尼等外部能量耗散外還存在內部能量耗散，例如熱彈性阻尼(thermoelastic damping,TED)，或內摩擦[1].實踐證明，可以通過完美的設計和制造工藝最大限度地降低或消除外阻尼引起的耗量能散.但是，TED 是系統所固有的能耗機制，是在結構周期振動時材料內部的熱彈性耦合變形引起的內部阻尼，它不能通過改善外部條件而消除.因此，TED 有時會成為微電機系統的主要能耗形式，或影響品質因子的重要因素.要在理論上分析微/納結構在振動中的熱彈性阻尼，就需要基于熱--彈耦合動力學理論建立結構振動的數學模型，通過聯立求解熱--彈耦合的運動方程和熱傳導方程獲得系統的動力響應，從而求得表征熱彈性阻尼的品質因子.TED 與微結構的材料性質、幾何尺寸、支承條件以及環境溫度等諸多因素密切相關.因此，定量地分析和預測TED 對高品質微/納諧振器的研究和設計具有重要意義.

大多數諧振器力學模型可以簡化為彈性梁或彈性板.最早關于諧振器熱彈性阻尼的理論研究當屬Zener 的工作[1].早在1938 年，Zener 首先基于Euler-Bernoulli (E-B)梁理論和單相耦合的準一維熱傳導方程求得了矩形截面微梁諧振器的TED近似解析解[1]，被稱為著名的Zener 公式.2000 年，在Zener工作的基礎上Lifshitz 和Roukes(L-R)[2]首次給出了E-B 微梁的準一維熱傳導方程的解析解，進而求解振動方程的特征值問題，獲得了更為精確的矩形截面微/納梁諧振器逆品質因子和頻移的解析解，被稱為著名的L-R 公式.Zener 和L-R 的杰出工作對后來關于微結構的TED 的進一步研究具有重要的參考和借鑒作用，他們的成果在后來研究者的論文中被高頻次引用.研究者采用不同的結構變形理論、熱傳導理論和分析方法開展了關于微/納梁板諧振器的熱彈性阻尼研究.關于逆品質因子的計算主要采用能量法[1]和復頻率法[2].

基于經典板理論，文獻[3-15]報道了均勻各向同性材料彈性薄板諧振器的TED[3-15].Nayfeh 和Youns[3]采用攝動法求解了引入剛度阻尼的振動方程，獲得了受靜電載荷和殘余面內應力共同作用的矩形微板TED 的近似解析解.Sun 等[3-5]、Ali 和Mohammadi[6]采用復頻率法分別研究了微圓/環板的TED，給出了L-R 形式的解析解.Salajeghe 等[7]在運動方程的建立中考慮了變形的幾何非線性，定量分析了軸對稱振動微圓板諧振器的振幅對熱彈性阻尼的影響規律.Li 等[8]采用準一維熱傳導方程，通過計算一個周期內彈性板損失的能量與結構內存儲的最大彈性勢能之比求得了矩形和圓形薄板諧振器L-R 形式的逆品質因子解析解.在此基礎上Fang 等[9-10]分別采用二維和三維熱傳導方程研究了軸對稱自由振動的圓板和矩形板諧振器的TED,分析了熱傳導方程的簡化對熱彈性阻尼值的預測精度的影響.考慮熱傳導過程中溫度熱流運動的延滯效應(或波動效應),文獻[11-14] 利用廣義熱傳導理論(或非傅里葉熱傳導理論)研究了微/納板諧振器的熱彈性耦合振動響應.Sharma 等[11]基于Lord-Shulman 廣義熱傳導理論和經典板理論研究了非軸對稱自由振動圓板諧振器的熱彈性阻尼，分析了延滯時間參數對熱彈性阻尼的影響.在此基礎上，Sharma 和Grover[12]還進一步研究了由于空隙率變化而引起的延滯效應對板式諧振器的影響.Guo 等[13]采用雙向延滯廣義熱傳導模型研究了圓板諧振器的熱彈性阻尼，采用復頻率法獲得逆品質因子的L-R 形式解析解.數值結果發現，在廣義熱傳導理論下逆品質因子隨著厚度的減小出現極小值.Grover 等[14]基于Kelvin-Voigt 材料模型下的廣義熱黏彈性理論研究了均勻微圓板諧振器的阻尼特性，定量地分析了機械松弛時間和熱松弛時間對逆品質因子的影響.

考慮材料性質沿單軸方向階梯變化，Bishop和Kinra[15-16]最早采用能量法開展了關于復合材料層合微/納諧振器的熱彈性阻尼研究.在表面邊絕熱條件和內部界面處的非完善協調條件下，求得了準一維熱傳導方程的解析解，并通過計算一個振動周期內由不可逆傳熱產生的總熱量與彈性總彈性勢能之比給出了諧振器的逆品質因子[15].并定量地分析了對稱鋪設的三層矩形微板的彈性熱動阻尼特性[16].Sun 等[17]基于準一維熱傳導方程，采用復頻率法研究了對稱鋪設的三層微圓板在軸對稱自由振動下的熱彈性阻尼.采用與文獻[16]中相同的能量方法，忽略溫度梯度在面內的變化，Zuo 等[18]研究了雙層微板的熱彈性阻尼，給出了周邊夾緊矩形板和圓板的逆品質因子解析解.考慮溫度梯度在面內的變化，Liu 等[19-20]采用格林函數法分別求解了二維和三維熱傳導方程，利用能量法獲得了雙層圓板和非對稱鋪設三層矩形板在給定邊界條件下的逆品質因子的級數形式解析解.這里需要說明的是，文獻[18-20]通過引入物理中面的概念消去了幾何中面的面內位移，簡化了幾何方程.但是，當材料性質變化關于幾何中面非對稱時板內將會產生熱薄膜力，而熱薄膜力對熱彈性阻尼具有貢獻[29].顯然，在上述文獻中忽略了熱薄膜力對熱彈性阻尼的影響.

近年來，亦有少量的文獻研究了材料性質沿厚度連續變化的功能梯度材料(functionally graded material,FGM)微/納梁板結構的熱彈性阻尼.Azizi 等[21]研究了由硅和壓電材料為組分材料復合而成的兩端夾緊的矩形截面FGM微梁橫向自由振動時的TED.采用Galerkin 法求得了FGM 微梁的復頻率進而獲得反映熱彈性阻尼的逆品質因子.數據結果表明可以通過壓電材料組分的合理調節來降低系統的TED.基于一階剪切變形理論和單向耦合的準一維熱傳導方程，Emami 和Alibeigloo[22]研究了材料性質沿橫向冪函數變化的FGM 微梁的TED.將熱傳導方程的系數和溫度場同時展開成關于橫向坐標的泰勒級數形式，獲得了變系數熱傳導方程的級數形式的解析解，進而獲得了兩端簡支FGM 微梁的TED.在此基礎上，文獻[22]的作者中又基于一階剪切變形板理論和修正的應變梯度理論，研究了四邊簡支功能梯度中厚度矩形板的熱彈性阻尼[23].但是，通過研究發現其中給出的關于均勻材料方板的熱彈性阻尼的計算結果是錯誤的.由此無法肯定其中關于功能梯度材料微板熱彈性阻尼數值結果的可靠性.考慮材料性質在厚度方向按指數函數分布特殊情形，Zhong 等[24]采用解析方法獲得了功能梯度E-B 微梁的熱彈性阻尼解析解.但是，作者在應變的計算中卻直接忽略了由于材料性質在橫向的非均勻性而導致的拉--彎耦合變形.許新等[25-26]基于單向耦合的準一維熱傳導方程研究了材料性質在橫向任意連續變化的FGM E-B 微梁的熱彈性阻尼.發展了求解橫向非均勻微梁的復雜變系數熱傳導方程的分層均勻化方法，從而獲得FGM 微梁熱彈性阻尼的半解析解,其中精確地考慮了熱軸力對TED 的影響.最近又將分層均勻化方法成功地應用于功能梯度微板的熱彈性阻尼的求解[27-28].研究了材料梯度指數、環境溫度、邊厚比、振動模態等對熱彈性阻尼的影響規律，首次定量地分析了熱薄膜力對熱彈性阻尼的影響程度.

文獻調研表明，截至目前關于均勻微板熱彈性阻尼的解析解是基于Kirchhoff經典彈性板理論推導出的，其中忽略了橫向剪切變形對振動響應的影響.本文基于考慮一階剪切變形效應的Mindlin 中厚板理論和單向耦合熱傳導理論建立四邊簡支均勻材料微板的熱--彈耦合自由振動的控制微分方程.消去兩個獨立的轉角變量將振動方程轉化為只包含撓度函數和熱彎矩的四階偏微分方程.忽略溫度梯度在面內的變化，在上下表面絕熱邊界條件下求得了用撓度函數表示的熱傳導方程的解析解.從而將包含熱彎曲內力的結構振動方程轉化為只包含撓度振幅的偏微分方程.利用特征值問題的數學相似性,在四邊簡支條件下獲得了用均勻無阻尼Kirchhoff微板的固有頻率表示的Mindlin 矩形微板的復頻率，從而由復頻率法給出了反映熱彈性阻尼水平的逆品質因子解析解.通過數值結果定量地分析了剪切變形對熱彈性阻尼值的影響程度.

1 問題的控制方程

1.1 運動方程

考慮長度a、寬度為b、厚度為h的等厚度矩形微/納板.選取直角坐標系，板在空間的區域定義為：0 ＜x＜a，0 ＜y＜b，-h/2 ＜z＜h/2 (見圖1).基于Mindlin 板理論[29-33]，板內任意一點的位移分量可表示為

其中，t為時間；u，v和w分別為板內任意一點沿x，y和z軸方向的位移分量；φx和φy為中面法線分別繞y軸和x軸的轉角.

圖1 微板的幾何尺寸和坐標示意圖Fig.1 Schematic illustration of a micro plate and the coordinate system

將式(1)代入彈性力學的幾何方程可得板的應變場

根據胡克定律得到應力場

其中，θ(x,y,z,t)=T(x,y,z,t)-T0為變溫場，T為瞬態溫度，T0為初始溫度.E，v和α 分別為彈性模量、泊松比和熱膨脹系數.這里的溫度場是由熱--彈耦合振動引起的.

將彈性力學空間問題應力形式的運動方程在厚度進行靜力等效，并利用上下表面的應力邊界條件可得用等效內力和內力矩表示的Mindlin 板[29-33]自由振動微分方程

其中Ii為等效慣性參數，具體表示為

ρ 為板的質量密度.式(4)～式(6)中的等效內力和內力矩可分別表示為

將式(8)代入式(4)和式(5)，分別關于坐標x和y求導后相加，利用式(6)可得

將式(8d)和式(8e)代入式(6)可得

最后，將式(11)代入式(9),可得只用橫向位移w0表示的運動方程

上式中熱彎矩由熱--彈耦合產生的變溫場確定.

1.2 熱傳導方程

根據單向耦合的熱彈性動力學理論，可給出微/納板的準一維熱傳導方程[3-8]

其中，κ，C分別是熱傳導系數和比熱，e是體積應變.Mindlin 板的體積應變可表示為

將式(14)代入式(13)，可得單向耦合的準一維熱傳導方程

2 熱彈性阻尼的求解

2.1 自由振動響應

假設位移場合溫度場的動態響應同步，則系統的熱--彈耦合自由振動響應可假設為調和形式

可得無量綱控制方程

其中

其中結構振動方程通過無量綱的熱彎矩振幅mT與熱傳導方程耦合.

2.2 熱傳導方程的求解

方程(20)的通解為

微諧振器一般在恒溫真空環境內工作,可認為沒有外界的熱流輸入.而由于熱--力耦合振動在諧振器內部產生的變溫場幅度很小,因而熱流強度也很低.于是可假設在微板的上下表面沒有熱流通過.這樣，在微板的上下表面通常給出絕熱邊界條件[2-28]

利用邊界條件(24)可確定系數A=qK/[pcos(p/2)]和B=0.于是得到特解

2.3 熱彈性阻尼

將式(25)代入熱彎矩的計算公式可得

將式(26)代入式(19),并利用式(11)的無量綱形式

可得只用無量綱撓度表示的結構振動微分方程

微分方程(29)可記為標準形式

其中

另外，由彈性薄板理論可得不考慮熱彈性阻尼的Kirchhoff板的振動微分方程[29]

對于四邊簡支矩形板，Kirchhoff薄板理論下方程(32)的邊界條件可表示為

可以證明(具體證明見附錄)，對于四邊簡支Mindlin 矩形板，振動微分方程(30)的邊界條件亦可表示為

由微分方程邊值問題(30)，(34)與問題(32)，(33)之間的相似性，可得它們的特征值之間的關系

于是，將式(31)代入式(35)，可得特征方程

其中R=μ00-(μ2+csf).由式(27)知f(Ω)是關于頻率Ω 的超越函數，因此方程(36)具有很強非線性.為了簡化求解,可近似地令f(Ω)=f(Ω0)[2-14]，其中Ω0為等溫(不考慮熱彈性阻尼)Mindlin 微板的固有頻率.于是，由式(36)可容易求得用無阻尼Kirchhoff板的頻率表示的具有熱彈性阻尼的Mindlin 微板的復頻率

其中，四邊簡支等溫Kirchhoff矩形微板的固有頻率的解析解為

在式(37)中令f=0 可得

式(39)給出了等溫Mindlin 微板的固有頻率與相應的Kirchhoff板的固有頻率之間的解析關系.將文獻[29]中的式(12.2.16)進行無量綱運算后所得結果與式(39)完全相同.

先由式(39)求得Ω0，然后代入式(27)可得到f=f(Ω0),最后將其代入式(37)即可得到熱--彈耦合自由振動Mindlin 微板的復頻率.根據復頻率法[2-8]即可得到用逆品質因子表示的Mindlin 微板的熱彈性阻尼解析解

其中Re(Ω)和Im(Ω)分別是復頻率的實部和虛部.

可以證明，對于簡支多邊形微板式(33)和式(34)仍然成立[29].因此，式(37)和式(40)也可以計算直邊多邊形中厚度微板的復頻率和熱彈性阻尼.

如果在式(36)中令μ2=cs=0，f=f()，則可得到Kirchhoff微板的復頻率[3-8]

利用式(27)和式(41)可得文獻[8]中 Kirchhoff微板熱彈性阻尼L-R[2]形式的解析解

3 數值結果與討論

數值計算時微板的材料分別選取碳化硅(SiC)、氮化硅(Si3N4)、金屬鋁(Al)和金屬鎳(Ni)[23,29],具體材料參數見表1.環境溫度為T0=300 K.

表1 微板的材料性質(T0=300 K)[23,25]Table 1 Material properties of the micro plate(T0=300 K)[23,25]

首先，針對碳化硅(SiC)方板，在表2 中給出了在不同振動模態下，基于Mindlin 板理論得到的熱彈性阻尼(40)與文獻[8]中采用經典理論所得到的解析解(42)的比較.由表中的數據可見，隨著微板的邊/厚比的減小(或邊長的減小)以及模態階數的增大Mindlin 微板的逆品質因子與Kirchhoff微板的逆品質因子之間的相對誤差Error=100%(-逐漸增大.在a/h=10 和基頻振動下，兩種理論預測結果之間的相對誤差僅為2.66%,而當微板以高階模態(r,s)=(3,3)振動時，相對誤差卻達到了8.94%.對應同樣的幾何參數、材料性質參數和振動模態，在表3 中給出了文獻[23]中的數值結果.與表2 對比后可見，文獻[23]中的結果與本文結果相差甚遠.而且在表3 中利用文獻[8]中的Kirchhoff微板的解析公式(42)計算得到的數值結果與本文利用同樣公式計算結果亦有很大差別.特別是在邊厚比a/h=10、振動模態為(3,3)的條件下，本文所得Mindlin 板和Kirchhoff板的TED 分別為8.04×10-5和8.83×10-5，相對誤差為8.94%；而在文獻[23]中的相應結果則分別為6.4×10-5和1.30×10-4，兩種理論的預測值之間的相對誤差竟然達到了50.7%.而對于長厚比a/h=5 的正方形微板,在振動模態為(3,3)的條件下，本文預測的兩種理論下的熱彈性阻尼之間的相對誤差也僅為28.4%.

表2 兩種板理論下正方形形陶瓷(SiC)微板的熱彈性阻尼(Q-1×104)比較(a=b,h=1 μm)Table 2 Comparison of TED(Q-1×104)of the square micro plate of ceramic(SiC)based on both of the Kirchhoffand Mindlin plate theories(a=b,h=1 μm)

表3 具有不同邊/厚比的正方形微板在前六階模態下的熱彈性阻尼比較TED(Q-1×103)(h=1 μm)[23]Table 3 Comparison of TED(Q-1×103)of square microplate with different side-to-thickness ratios for the first six modes(h=1 μm)[23]

為了進一步驗證本文解析解的可靠性，選擇a=b=100 μm 的陶瓷(SiC)微板，在圖2 中分別給出了由式(34)和式(40)預測的TED 和采用有限元法得到數值結果的比較.有限元法是基于全耦合熱彈性動力學理論、采用空間六面體單元獲得的.由圖中結果可見，在板厚較小時二者吻合很好.隨著厚度的增大，有限元解答與本文解析解稍有偏離.這是因為在板理論中，將側面上的應力邊界條件沿著厚度靜力等效為幾何中面周邊上的合力矩邊界條件，根據圣維南原理隨著厚度的增大板的自然邊界條件與空間彈性力學的應力邊界條件的差別將會增大.實際上，位移邊界條件的差別也在增大.

圖2 分別由本文方法和有限元法預測的正方形陶瓷(SiC)微板的熱彈性阻尼值的比較Fig.2 Comparison between the values of TED of square ceramic micro plate by FEM and present approach,respectively

為了定量地分析剪切變形對熱彈性阻尼預測值的影響，給定不同邊/厚比a/h，在圖2 中繪出了正方形金屬(Ni)板諧振器在基頻振動時的熱彈性阻尼隨板厚的連續變化曲線.并在表4 中列出了圖2 中熱彈性阻尼的最大值Q-1max和對應的臨界厚度hcr.數值結果表明，Kirchhoff板理論的熱彈性阻尼預測值始終大于Mindlin 板理論的預測值.兩種理論預測的熱彈性阻尼之間的差值在臨界厚度附近十分顯著.另外，隨著微板的邊/厚比的增大，Mindlin 微板的熱彈性阻尼最大值單調增加，而Kirchhoff微板的熱彈性阻尼最大值卻保持不變.從表4 中數據可見，Mindlin 微板的臨界厚度大于Kirchhoff微板的臨界厚度.為了更加清晰地反映兩種理論預測值的差值的變化，在圖3 中繪出了差值-隨厚度的連續變化曲線.圖中結果再次清楚地表明在臨界厚度附近差值變化最大.但是，在不考慮熱彈性阻尼的情況下由式(39)可知Mindlin 板的無量綱固有頻率只與邊厚比有關.在給定邊厚比后無量綱頻率為常數.因此，兩種理論預測的等溫板的無量綱故有頻率差值為常數.

表4 具有不同邊/厚比正方形微板(Ni)的最大熱彈性阻尼和相應的臨界厚度(一階模態)Table 4 The maximum TED and related critical thickness of square micro plate(Ni)for different side-to-thickness ratios(the first mode)

圖3 兩種板理論下正方形金屬(Ni)微板的熱彈性阻尼隨厚度的變化曲線(一階模態)Fig.3 Curves of the TED versus the plate thickness of a square micro metal plate(Ni)based on the two plate theories(in the first mode)

圖4 中繪出了表1 中所列4 種材料的正方形微板在一階模態振動時的熱彈性阻尼隨厚度的連續變化曲線.其中給出了兩種板理論預測結果的對比.從圖中可見，金屬微板的熱彈性阻尼明顯大于陶瓷微板的熱彈性阻尼.并且兩種理論下金屬微板的最大熱彈性阻尼之間的差值比陶瓷微板顯著.圖5 給出了a/h=10 時，具有不同長寬比的矩形陶瓷(SiC)微板在基頻振動下的TED 與厚度的關系曲線.結果表明，隨著長寬比的增大(相當于彎曲剛度增大)兩種板理論預測的熱彈性阻尼的最大值之間的差值顯著增大，而經典理論下的最大值卻保持不變.圖6 則為a/h=10 的中等厚度板金屬(Al)微板在前四階振動模態下的熱彈性阻尼隨厚度的變化曲線.其中的變化規律與圖5 相似.隨著振動模態階數的增大，板的固有頻率增大.這相當于彎曲剛度在增加，從而導致剪切變形對熱彈性阻尼的影響加大.

圖4 給定不同邊厚比時，兩種板理論預測的正方形金屬(Ni)微板的熱彈性阻尼差值隨厚度的變化曲線(一階模態)Fig.4 Difference between the values of TED of a square metal(Ni)micro plate evaluated by the two plate theories varying with the thickness for some specified values of a/h(in the first mode)

圖5 一階模態下不同材料正方形微板的熱彈性阻尼與板厚之間的關系曲線Fig.5 Relationship curves of TED versus the plate thickness of the square micro plates of different materials in the first mode

圖6 具有不同長寬比(a/b)的矩形陶瓷(SiC)微板的熱彈性阻尼與板厚之間的關系曲線(一階模態)Fig.6 TED versus the plate thickness in ceramic(SiC)micro rectangular plate with different length to width ratio(in the first mode)

圖7 不同振動模態下金屬(Al)微板的熱彈性阻尼隨厚度的關系曲線(a/h=10)Fig.7 TED versus the plate thickness in a square metal(Al)micro plate in different vibration modes(a/h=10)

4 結論

基于一階剪切變形理論和單向耦合準一維熱傳導理論，采用解析方法研究了四邊簡支中厚度矩形微板諧振器的熱彈性阻尼.利用微分方程邊值問題之間的相似性，推導出了用無阻尼Kirchhoff微板的固有頻率表示的Mindlin 微板熱--彈耦合自由振動復頻率的解析解，從而采用復頻率法給出了表征Mindlin 微板熱彈性阻尼的逆品質因子解析解.分別選定兩種金屬和兩種陶瓷材料微板，通過參數變化的數值結果定量地分析了剪切變形對熱彈性阻尼預測值的影響規律.數值結果表明，經典板理論預測的熱彈性阻尼值大于一階剪切變形板理論的預測值.通過研究還發現，兩種理論的預測的熱彈性阻尼值的差別在微板的臨界厚度附近最為顯著.

致謝根據審稿專家的建議,在修改稿中增加了本文解析解與有限元法所得數值結果的比較.揚州大學建筑科學與工程學院的張志超博士和他指導的碩士研究生曹靜同學在有限元數值計算中給予作者重要的技術幫助.在此表示感謝.