雙層微梁結(jié)構(gòu)諧振器中波動(dòng)溫度場的研究?

左萬里劉 璇楊孫法子

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

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

本文以雙層Euler-Bernoulli微梁結(jié)構(gòu)機(jī)械諧振器為研究對(duì)像，微梁兩端固定，如圖1所示。微諧振器鍍層常采用真空電子束蒸發(fā)淀積加工，層與層之間具有良好的熱邊界條件[12－14]。基于格林函數(shù)法，建立一維熱傳導(dǎo)條件下，內(nèi)部具有熱源項(xiàng)時(shí)雙層微梁中的波動(dòng)溫度場分布函數(shù)。由溫度場函數(shù)，計(jì)算了在諧振器中波動(dòng)溫度場的變化幅值。

圖1 雙層微梁結(jié)構(gòu)示意圖

1 雙層微梁內(nèi)控制方程

1.1 振動(dòng)控制方程

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

式中:w0(x)為梁在靜載荷作用下的位移函數(shù)。考慮彈塑性材料的熱耦合效應(yīng)，簡諧變化的應(yīng)力/應(yīng)變會(huì)產(chǎn)生溫度變化，而溫度場的波動(dòng)又會(huì)導(dǎo)致熱應(yīng)變的變化，則波動(dòng)溫度場產(chǎn)生的熱應(yīng)變可表示為:

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

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

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

每一層中軸向正應(yīng)力可記為[11]:

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

當(dāng)梁作純彎曲振動(dòng)時(shí)，橫截面中總的正應(yīng)力之和為零，即:

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

1.2 熱傳導(dǎo)方程及其熱邊界條件

對(duì)于Euler-Bernoulli微梁，垂直方向上的溫度梯度遠(yuǎn)大于水平方向上的溫度梯度，因此一般只考慮厚度方向上的熱傳導(dǎo)。則雙層Euler-Bernoulli微梁中含有內(nèi)部熱源項(xiàng)時(shí)，一維熱傅立葉傳導(dǎo)下的溫度場控制方程可記為[16－17]:

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

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

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

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

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

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

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

平面結(jié)構(gòu)中格林函數(shù)的表達(dá)式可記為::

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

式中:

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

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

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

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

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

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

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

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

3 結(jié)果與分析

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

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

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

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

3.1 溫度場實(shí)部

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

圖2 溫度場實(shí)部有限元計(jì)算結(jié)果

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

圖3 解析模型溫度場實(shí)部計(jì)算結(jié)果

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

圖4 有限元模型與解析模型結(jié)果實(shí)部比較

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

圖5 溫度場實(shí)部誤差分析

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

3.2 溫度場虛部

由于溫度場虛部與振動(dòng)產(chǎn)生的應(yīng)力/應(yīng)變之間存在相位差，從而會(huì)帶來材料的熱松弛，進(jìn)而導(dǎo)致能量的損耗[20－22]。因此，對(duì)于微機(jī)械諧振器，波動(dòng)溫度場的虛部也是研究的重點(diǎn)。

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

圖6 溫度場虛部有限元計(jì)算結(jié)果

圖7 解析模型溫度場虛部計(jì)算結(jié)果

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

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

與圖5相似，圖9給出了最下側(cè)(z＝0μm處)沿長度方向上，波動(dòng)溫度場虛部解析模型與數(shù)值模型之間的誤差值與誤差率。從圖中可知，誤差值與誤差率最大值都出現(xiàn)在x＝0μm處，最大誤差率為18.5%。在x＝30μm至60μm之間，誤差率在－5%至5%之間波動(dòng)。分析其原因，x＝30μm至60μm之間的誤差率波動(dòng)主要是由于數(shù)值模型的截?cái)嗾`差來的。而在固定端(x＝0μm)，誤差值為0.026 K，遠(yuǎn)大于數(shù)值模型計(jì)算時(shí)保留的最后一位有效值(0.000 01 K，5位有效數(shù)值)。因此，固定端的誤差不再是由數(shù)值模型截?cái)嗨鶐淼摹Ｔ诮馕瞿Ｐ椭校捎谥豢紤]了厚度方向上的一維熱傳導(dǎo)，而在固定端附近，長度方向上的溫度梯度只略小于厚度方向上的梯度，(如圖2，圖6所示)，一維熱傳導(dǎo)假設(shè)不再成立，從而會(huì)帶來誤差。

圖9 溫度場虛部誤差分析

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

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

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

圖10 激振頻率對(duì)溫度場的影響

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

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

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

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

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

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

圖13 材料比對(duì)波動(dòng)溫度場的影響

4 結(jié)論

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