雙層微梁結構諧振器中波動溫度場的研究?

左萬里劉 璇楊孫法子

(1.寧波大學機械工程與力學學院，浙江 寧波 315211；2.浙大寧波理工學院計算機與數(shù)據(jù)工程學院，浙江 寧波 315100)

微電子機械系統(tǒng)(MEMS，Micro electro mechanical systems)是在微機械加工與微電子技術基礎上發(fā)展起來的新興交叉領域，涉及熱力學、材料學、光學和電磁學等多種學科和工程技術[1]。微機械諧振器是MEMS中的一類典型器件，其通常工作在一階固有頻率，實現(xiàn)電信號與機械信號之間的相互轉換[2－4]。微小尺度是MEMS的重要特征，當尺度縮小到微米級時，會產生微尺度效應[5－7]，從而使許多物理現(xiàn)象與宏觀世界有很大差別。在宏觀器件中，重力等與尺寸三次方(L3)有關的力起主導作用，熱應力等與尺寸一次方(L)有關的力可忽略不記。然而，在微諧振器中，由于尺寸的減小，熱應力對器件的影響遠遠大于重力。微諧振器中振動發(fā)熱帶來的熵增，會產生熱彈性阻尼，從而降低了器件的品質因數(shù)與靈敏度[8－9]。因此，為了精確計算熱應力的大小與對器件的作用，需要建立諧振器件中的波動溫度場[10－11]。

本文以雙層Euler-Bernoulli微梁結構機械諧振器為研究對像，微梁兩端固定，如圖1所示。微諧振器鍍層常采用真空電子束蒸發(fā)淀積加工，層與層之間具有良好的熱邊界條件[12－14]?；诟窳趾瘮?shù)法，建立一維熱傳導條件下，內部具有熱源項時雙層微梁中的波動溫度場分布函數(shù)。由溫度場函數(shù)，計算了在諧振器中波動溫度場的變化幅值。

圖1 雙層微梁結構示意圖

1 雙層微梁內控制方程

1.1 振動控制方程

當雙層Euler-Bernoulli微梁在角頻率為ω的激振力作用下進行簡諧振動時，位移函數(shù)可以表示為:

式中:w0(x)為梁在靜載荷作用下的位移函數(shù)?？紤]彈塑性材料的熱耦合效應，簡諧變化的應力/應變會產生溫度變化，而溫度場的波動又會導致熱應變的變化，則波動溫度場產生的熱應變可表示為:

式中:?m為第m層梁中相比與平衡溫度T0的波動溫度場，αm為第m層中材料的熱膨脹系數(shù)。

由于微梁是純彎曲振動，每一層中y向和z向正應力可以忽略，即σyy，m(z，t)＝σzz，m(z，t)＝0。則第m層中材料的體應變可記為[15]:

式中:νm為第m層中材料的泊松比。

每一層中軸向正應力可記為[11]:

式中:Em為第m層中材料的楊氏模量。

當梁作純彎曲振動時，橫截面中總的正應力之和為零，即:

將式(4)代入式(5)，可得中性面的位置為:

1.2 熱傳導方程及其熱邊界條件

對于Euler-Bernoulli微梁，垂直方向上的溫度梯度遠大于水平方向上的溫度梯度，因此一般只考慮厚度方向上的熱傳導。則雙層Euler-Bernoulli微梁中含有內部熱源項時，一維熱傅立葉傳導下的溫度場控制方程可記為[16－17]:

式中:km為第m層梁中材料的熱傳遞系數(shù)，Cm為第m層中材料單位體積的比熱，內部熱源項為:

雙層微梁中溫度場θm(x，y，z，t)＝?m(x，y，z)eiωt的邊界條件可分為如下情況:梁的下表面為第二類齊次熱邊界條件，即沒有熱量的傳遞，因此:

此外，第一層梁與第二層梁具有良好的熱接觸，在界面處沒有熵增，即:

同時，第二層梁上表面熱邊界條件可以表示為:

2 格林函數(shù)法求解波動溫度場方程

由格林函數(shù)法可得每一層梁中溫度場方程為[18]:

式中:F1(z′)，F(xiàn)2(z′)分別表示第一層梁和第二層梁中初始時刻的溫度，由于初始時刻梁未振動，每一層中的溫度都等于平衡溫度T0，即有F1(z′)＝F2(z′)＝T0。g1(z′，τ)和g2(z′，τ)分別為第一層梁和第二層梁內部熱源項，由式(8)和(4)可得，即:

平面結構中格林函數(shù)的表達式可記為::

將式(13)和式(14)代入式(12)，可得:

式中:

?m，n和βn分別為雙層微梁溫度場控制方程式(7)對應特性問題的特征函數(shù)和特征向量。上述溫度場所對應的特征方程記為:

微分方程(19)的一般解為:

由邊界條件(9)，在z＝0處，有:

由式(21)可知，A1，n等于0。為了計算方便，可設任何一個不為零的系數(shù)，例如B1，n等于1，即相當于系數(shù)矩陣提取非零項Bm，n，這一點不會影響結果的一般性。則雙層梁中，每一層的特征函數(shù)可分別記為:

式中:βn、An和Bn為待定系數(shù)(z)和?2，n(z)滿足正交關系。

將式(22)和(23)代入到邊界條件式(10)和(11)中，可得:

由式(15)可知，格林函數(shù)法所求得溫度場方程分為三部分，其中只有第二項是周期變化的穩(wěn)態(tài)項，第一項和第三項是瞬時衰減項。當諧振器處于工作狀態(tài)時，衰減項將會消失，即諧振器中溫度場分布函數(shù)為:

將式(28)展開后，每一層中波動溫度場分別可記為:

3 結果與分析

SiC－Si在微機械諧振器中是常見的材料組合[14]，300 K(27℃)時，材料的力學特性如表1所示[15]。對于兩端固定Euler-Bernoulli微梁，當小變形時，其位移函數(shù)w0可用下列多項式來表示[19]

表1 300 K時材料的力學參數(shù)[16]

式中:l為梁的總長度，A0為最大振幅。為了計算與比較的方便，可設定最大振幅A0＝1μm。在諧振器實際工作過程中，振幅一般遠小于1μm。

為了驗證當前所得解析模型的有效性，本文利用ANSYS有限元數(shù)值模型計算了雙層微梁中的溫度場分布情況。在ANSYS中采用了solid226(3維20節(jié)點)單元進行了諧響應計算。由式(34)與式(35)可知，波動溫度場為一個復變量，其包含實部和虛部兩個部分。本文分別從實部和虛部兩個方面比較當前解析模型與有限元數(shù)值模型計算結果。

3.1 溫度場實部

圖2為有限元數(shù)值模型計算所得溫度場實部云圖。圖中所用雙層微梁，Si層厚度為5μm，SiC層厚度為5μm，總長度為200μm，諧振頻率為1×106Hz。顯然易見，溫度場云圖左右對稱。從圖中還可知，微梁最大振幅為1.005 35μm，最低溫度為26.769 9℃，波動溫度為－0.230 1℃(K)，最高溫度為27.187 8℃，波動溫度為0.187 8℃(K)，波動的溫度遠小于初始平衡溫度300 K。

圖2 溫度場實部有限元計算結果

運用MATLAB可繪制當前解析模型所得波動溫度場云圖，如圖3所示。比較圖2和圖3可得，當前解析模型所得結果與有限元數(shù)值模型所得結果具有很高的擬合度。固定端下側(z＝0μm處)為波動溫度場較高部分，上側(z＝10μm處)為較低部分。而中間部分與兩端正好相反，上側為較高部分，下側為較低部分。

圖3 解析模型溫度場實部計算結果

圖4給出了圖2，圖3中三個不同截面處，沿厚度方向上，當前解析模型與有限元數(shù)值模型計算結果的比較。從圖中可以知，上下表面處溫度場波動比較大，中性面處(z0＝6.07μm)溫度場為零。根據(jù)式(8)分析其原因，可知上下表面處的內部熱源項數(shù)值較大，而中性面處內部熱源項為零。

圖4 有限元模型與解析模型結果實部比較

由于上下表面處的溫度場波動較大，圖5給出了最下側(z＝0μm處)沿長度方向上，解析模型與數(shù)值模型之間的誤差值與誤差率，誤差率由式(37)可得

圖5 溫度場實部誤差分析

從圖5中可知，固定端(x＝0μm處)兩種模型間的誤差值為最大，但兩者之間的誤差率在x＝45 μm處達到最大－36.63%。分析其原因可知，由于ANSYS數(shù)值模型計算時只保留了6位有效數(shù)字(如圖2所示)，用prnsol命令提取節(jié)點上溫度時只保留了5位有效數(shù)字，即保留到小數(shù)點后3位。而在x＝45μm處波動溫度場解析模型計算結果為－0.009 468 K，有限元數(shù)值模型計算結果為－0.006，兩種模型計算結果間的誤差值為0.003 468 K，可知數(shù)值模型計算過程中的斷截誤差會對最后結果有很大的影響。

3.2 溫度場虛部

由于溫度場虛部與振動產生的應力/應變之間存在相位差，從而會帶來材料的熱松弛，進而導致能量的損耗[20－22]。因此，對于微機械諧振器，波動溫度場的虛部也是研究的重點。

圖6給出了有限元數(shù)值模型所得波動溫度場虛部的云圖，圖7為當前解析模型所得云圖。兩種模型所得結果具有較高的擬合度。進一步比較實部與虛部云圖，二者數(shù)值相近，云圖基本相似，即實部最大波動處也為虛部最大波動處。

圖6 溫度場虛部有限元計算結果

圖7 解析模型溫度場虛部計算結果

圖8為圖6，圖7中波動溫度場虛部，三個不同截面處沿厚度方向上兩種模型計算結果的比較。從圖中可知，在x＝50μm和100μm處，兩種模型計算結果基本重合。在固定端x＝0μm處，兩種模型計算結果存在一定的誤差，且越遠離中性面，誤差值越大。

圖8 有限元模型與解析模型結果虛部比較

與圖5相似，圖9給出了最下側(z＝0μm處)沿長度方向上，波動溫度場虛部解析模型與數(shù)值模型之間的誤差值與誤差率。從圖中可知，誤差值與誤差率最大值都出現(xiàn)在x＝0μm處，最大誤差率為18.5%。在x＝30μm至60μm之間，誤差率在－5%至5%之間波動。分析其原因，x＝30μm至60μm之間的誤差率波動主要是由于數(shù)值模型的截斷誤差來的。而在固定端(x＝0μm)，誤差值為0.026 K，遠大于數(shù)值模型計算時保留的最后一位有效值(0.000 01 K，5位有效數(shù)值)。因此，固定端的誤差不再是由數(shù)值模型截斷所帶來的。在解析模型中，由于只考慮了厚度方向上的一維熱傳導，而在固定端附近，長度方向上的溫度梯度只略小于厚度方向上的梯度，(如圖2，圖6所示)，一維熱傳導假設不再成立，從而會帶來誤差。

圖9 溫度場虛部誤差分析

3.3 波動溫度場的影響因數(shù)

基于本文所建立的解析模型，進一步分析激振頻率，微梁的長度，材料比等對波動溫度場的影響。

圖10為x＝100μm截面處，不同激振頻率時，沿厚度方向上的溫度場虛部變化圖。圖中雙層微梁的長度為200μm，Si和SiC厚度都為5μm。由圖可知，當激振頻率為1×104Hz和1×108Hz時，波動溫度場變化極??；在1×106Hz時，波動溫度場幅值最大。

圖10 激振頻率對溫度場的影響

根據(jù)式(8)和式(7)可知，當激振頻率較小時，內部熱源項發(fā)熱較小，且每個周期的時間較長，熱量有足夠的時間從高溫部分傳遞到低溫部分，每半個振動周期內材料中的溫度場相當于等溫。且由于振動的交替性，對于某一局部，變?yōu)楦邷貢r輸出的熱量與變?yōu)榈蜏貢r輸入的熱量基本當?shù)?，每一個周期內變化的熱量總合基本為0，所以波動溫度場變化不大。

當激振頻率較大時，雖然內部熱源項較大，但頻率過大時，每個周期的時間就較短，由于振動的交替性，高溫部分的熱量還沒來得及向低溫部分傳遞，高溫部分就逆轉為低溫部分，熱量就要回傳，此時材料相當于絕熱，單位周期內總的熱量變化也基本為0，則整體的波動溫度場變化也不大。而在1×105Hz和1×106Hz頻段，內部熱源項和單位周期時間適中，溫度場波動相對較大。

圖11和12分別為在中間位置和l/4長度兩個截面處，波動溫度場虛部沿厚度方向上的變化曲線。圖中微梁的激振頻率為1×106Hz，Si和SiC厚度都為5μm，梁的長度不同。由式(36)可知，在中間位置x＝l/2處，不同長度微梁位移函數(shù)w0的振幅都相同，其值為A0。同理，在中間位置x＝l/4處，不同長度微梁位移函數(shù)w0的振幅都為9A0/16。而由圖中可知，波動溫度場隨梁的總長度變化而不同，且隨著總長度的增加，波動溫度場幅值減小。

圖11 在x＝l/2處總長度對溫度場的影響

圖12 在x＝l/4處總長度對溫度場的影響

圖13為不同材料比情況下，中間截面處，波動溫度場虛部沿厚度方向上的變化曲線。圖中微梁的激振頻率為1×106Hz，總長l＝200μm，總厚為z2＝10μm，設定Si層厚度z1變化。從圖中可知，隨著Si層材料的體積比變化，波動溫度場也變化，且為非線性變化。

圖13 材料比對波動溫度場的影響

4 結論

本文針對兩端固定雙層Euler-Bernoulli微梁結構諧振器，利用格林函數(shù)求解了一維傅立葉熱傳遞條件下的溫度場分布函數(shù)，得出了波動溫度場的解析模型，通過與有限元數(shù)值模型計算的比較，驗證了本文所得解析模型的正確性。主要結論如下:①由于振動發(fā)熱而產生的波動溫度場相比于平衡溫度非常小，僅為平衡溫度的0.06%，則由于波動溫度場而產生的熱應變也非常??；②數(shù)值模型在計算溫度場實部時，由于計算時有效數(shù)字的原因，會帶來截斷誤差，此為數(shù)值模型的局限性；③解析模型在固定端，由于水平方向上的溫度梯度也較大，會產生18.5%的誤差，此為解析模型的局限性；④激振頻率對波動溫度的影響比較大，激振頻率過低或過高時，微梁中的波動溫度場基本為零；⑤微梁的幾何尺寸、材料比對波動溫度場也有影響。