非完全覆蓋雙層微梁諧振器熱彈性阻尼建模

楊龍飛 李普 葉一舟

摘要 熱彈性阻尼是決定微機(jī)械諧振器品質(zhì)因子上限的關(guān)鍵參數(shù)之一。以往熱彈性阻尼解析解只適用于完全覆蓋多層微梁結(jié)構(gòu)。由于制造工藝和實(shí)際功能需求，非完全覆蓋雙層梁為代表的復(fù)雜結(jié)構(gòu)形式更普遍?；诟道锶~傳熱定律，推導(dǎo)出非完全覆蓋雙層微梁諧振器熱彈性阻尼的解析解。同時(shí)利用數(shù)值方法和實(shí)驗(yàn)驗(yàn)證了該解析解的有效性。并分析了鍍層厚度、長(zhǎng)度和位置對(duì)熱彈性阻尼的影響。

關(guān)鍵詞 微機(jī)電系統(tǒng)諧振器; 熱彈性阻尼; 品質(zhì)因子; 熱彈性理論

引 言

微機(jī)械諧振器用途廣泛，常作為加速度傳感器、陀螺儀、濾波器、能量收集器等微機(jī)電系統(tǒng)的核心元器件。品質(zhì)因子是微諧振器性能參數(shù)的核心參數(shù)之一，與其頻率穩(wěn)定性、相位噪音及分辨率密切相關(guān)［1］。決定品質(zhì)因子的因素是微諧振器工作過程出現(xiàn)的各種能量損耗，主要包括空氣阻尼、支撐阻尼和熱彈性阻尼等阻尼機(jī)制。空氣阻尼可通過真空封裝消除，支撐阻尼可通過合理結(jié)構(gòu)設(shè)計(jì)降低到可忽略?？諝庾枘岷椭巫枘釋儆谕獠孔枘?，熱彈性阻尼是諧振器內(nèi)部產(chǎn)生的不可逆熱流導(dǎo)致的能量損失。熱彈性阻尼已被實(shí)驗(yàn)證實(shí)屬于固有阻尼［2］，只能通過合理的結(jié)構(gòu)設(shè)計(jì)優(yōu)化減小，而不能完全消除。因此，熱彈性阻尼決定了撓性微諧振器的品質(zhì)因子上限，成為近年來研究熱點(diǎn)之一［3］。

Zener最先認(rèn)識(shí)到熱彈性阻尼的重要性，自1937陸續(xù)發(fā)表了一系列成果，為熱彈性阻尼研究奠定了理論基礎(chǔ)。Zener利用三角函數(shù)級(jí)數(shù)疊加法得到矩形梁厚度方向的溫度場(chǎng)函數(shù)，最終推導(dǎo)出熱彈性阻尼的精確解。Lifshitz和Roukes（L?R）［4］利用矩形梁的復(fù)數(shù)形式溫度場(chǎng)函數(shù)，通過復(fù)頻率法得到熱彈性阻尼解析解。L?R解析解還可計(jì)算由于熱彈性阻尼造成的頻率偏移。L?R解析解與Zener解析解吻合很好，是后續(xù)研究熱彈性阻尼機(jī)理應(yīng)用最廣泛的兩個(gè)理論框架。矩形梁是微諧振器件常采用的結(jié)構(gòu)形式之一，Zener和L?R均以矩形梁為對(duì)象展開研究。Kumar和Haque［5］考慮靜態(tài)軸向拉應(yīng)力對(duì)微梁諧振器的影響，推導(dǎo)出熱彈性阻尼解析解。研究發(fā)現(xiàn)：軸向拉應(yīng)力可提高諧振頻率，同時(shí)降低熱彈性阻尼。Prabhakar和Vengallatore［6］開發(fā)了中空和開槽的單層微梁的熱彈性阻尼解析解，研究發(fā)現(xiàn)開槽可阻斷不可逆熱流擴(kuò)散，有效降低熱彈性阻尼。馬航空等［7］利用復(fù)頻率法推導(dǎo)出Mindlin矩形微板的熱彈性阻尼解析解。Kumar和Mukhopadhyay［8］提出了一種基于Moore?Gibson?Thompson廣義熱彈性理論和修正偶應(yīng)力理論分析微梁諧振器的熱彈性阻尼和動(dòng)力學(xué)特性。Gu等［9］研究了應(yīng)力非局部和高階應(yīng)變梯度效應(yīng)對(duì)微梁諧振器熱彈性阻尼的影響。

隨著MEMS制造工藝的精進(jìn)，多層結(jié)構(gòu)因其功能多樣性應(yīng)用越來越廣，比如：金屬膜經(jīng)常用于電極、質(zhì)量檢測(cè)器、光學(xué)反射、磁性單元和熱導(dǎo)體等［10］。SiO2層可用于提高微諧振器的溫度頻率性能等。Nourmohammadi等［11］提出了雙層微梁諧振器考慮厚度熱傳導(dǎo)的一維熱彈性阻尼解析解，首次發(fā)現(xiàn)SiO2/Si梁雙德拜峰現(xiàn)象并作了相應(yīng)分析。左萬里等［12］建立了雙層矩形板微諧振器熱彈性阻尼解析解。Yang等［13］推導(dǎo)得到考慮長(zhǎng)度和厚度熱傳導(dǎo)的熱彈性阻尼解析解，并總結(jié)出導(dǎo)致二維熱彈性阻尼模型與一維模型差異的兩個(gè)因素。

以上雙層及多層微梁的熱彈性阻尼解析解均以上下層完全覆蓋為出發(fā)點(diǎn)進(jìn)行研究。實(shí)際上，由于機(jī)械夾緊或電絕緣，鍍層很難完全覆蓋于基底層。Sandberg等［14］通過實(shí)驗(yàn)證實(shí)，即使很薄的金屬鍍層也將導(dǎo)致品質(zhì)因子劇烈下降。文中建議可在基底選擇性鍍膜，而非完全覆蓋基底層，以實(shí)現(xiàn)對(duì)熱彈性阻尼的有效控制。本文從理論層面推導(dǎo)出非完全覆蓋雙層微梁諧振器熱彈性阻尼解析解，并考慮了沿長(zhǎng)度和厚度兩個(gè)方向熱傳導(dǎo)產(chǎn)生的能量損失。該解析解可退化到Y(jié)ang等［13］提出的完全覆蓋雙層微梁熱彈性阻尼解析解。利用數(shù)值方法和實(shí)驗(yàn)驗(yàn)證了當(dāng)前解析解的有效性。通過分析鍍層厚度、長(zhǎng)度和位置對(duì)熱彈性阻尼的影響，為降低微諧振器熱彈性阻尼提供了一種思路。

1 理論模型的建立

1.1 動(dòng)力學(xué)方程及求解過程

非完全覆蓋雙層微梁諧振器結(jié)構(gòu)示意如圖1所示，基底層的長(zhǎng)度和厚度分別為L(zhǎng)1和h1，鍍層的長(zhǎng)度和厚度分別為L(zhǎng)2和h2。鍍層左端與基底層左端距離為a，且在寬度方向重合，均為矩形截面，組成材料均質(zhì)且各向同性。

鍍層改變了覆蓋區(qū)域的剛度，因此可將該結(jié)構(gòu)整體分為三個(gè)子區(qū)域：左側(cè)單層區(qū)域{1}、雙層區(qū)域{2}及右側(cè)單層區(qū)域{3}。圖1的直角坐標(biāo)系中，x，y和z分別代表梁的長(zhǎng)度、厚度和寬度方向。

簡(jiǎn)諧激勵(lì)下諧振器橫向垂直位移可表示為：

式中 w為諧振頻率，上標(biāo) p = 1，2，3，與左側(cè)單層區(qū)域、中間雙層區(qū)域和右側(cè)單層區(qū)域一一對(duì)應(yīng)。

振動(dòng)時(shí)歐拉?伯努利梁的曲率半徑遠(yuǎn)大于振幅，可忽略單層與雙層區(qū)域過渡位置處應(yīng)力集中的影響［6］。因此，動(dòng)力學(xué)控制方程可根據(jù)線彈性振動(dòng)理論方程設(shè)為［15］：

式中 （EI）{p}和（ρA）{p}為等效參數(shù)。

式（3）和（4）中，E1和E2分別表示基底層和鍍層材料的彈性模量。ρ表示材料密度，I，A分別表示截面慣性矩、截面面積。后續(xù)出現(xiàn)的下標(biāo)1和2均分別代表基底層和鍍層。

將式（1）代入式（2），化簡(jiǎn)后得到：

方程（5）的通解可設(shè)為：

懸臂和雙端固支是微梁諧振器最常使用的兩種支撐方式，其結(jié)構(gòu)邊界條件分別為：

此外，單層區(qū)域過渡到雙層區(qū)域需要滿足以下邊界連續(xù)性條件：

① 位移連續(xù)性：

② 轉(zhuǎn)角連續(xù)性：

③ 彎矩連續(xù)性：

④ 剪力連續(xù)性：

將幾何邊界條件（8）或（9）和連續(xù)條件（10）～（12）代入通解（7）中，得到12個(gè)線性方程。可將該12個(gè)方程組裝成矩陣形式：

為求非平凡解，M（β{p}）行列式值應(yīng)為0：

式（15）僅包含一個(gè)未知變量，即諧振頻率ω。式（15）是一個(gè)復(fù)雜超越方程，可利用二分法等尋根方法求解。求得ω后，利用式（6）和（14）可分別求解特征值β{p}和振型參數(shù)u。

1.2 溫度場(chǎng)求解

依據(jù)熱彈性耦合理論［16］，彎曲振動(dòng)時(shí)體積變化會(huì)產(chǎn)生不均勻溫度場(chǎng)。根據(jù)熵增原理，該溫度場(chǎng)存在一個(gè)平衡溫度T0，溫度場(chǎng)的相對(duì)增量函數(shù)可表示為：

式中 T{p}1和T{2}2是瞬時(shí)溫度場(chǎng)函數(shù)。

其中，h{p}0表示各子區(qū)域中性面位置：

根據(jù)胡克定律，子區(qū)域的應(yīng)力可表示為：

在高品質(zhì)因子MEMS諧振器中，與外界施加的應(yīng)力相比，由溫度變化產(chǎn)生的熱應(yīng)力小到可以忽略不計(jì)。因此，方程中軸向應(yīng)力可近似為：

通常，耦合溫度場(chǎng)是復(fù)數(shù)形式的，表示與施加的應(yīng)力存在相位差。溫度場(chǎng)沿長(zhǎng)度和厚度方向的二維熱傳導(dǎo)控制方程為：

1.3 熱彈性阻尼解析解

根據(jù)熱彈性耦合理論［16］，彈性應(yīng)變與應(yīng)力場(chǎng)同相，因此彈性應(yīng)變沒有能量損失。但熱應(yīng)變場(chǎng)與應(yīng)力場(chǎng)不同相，導(dǎo)致在每個(gè)振動(dòng)周期內(nèi)引起能量損失。根據(jù)Bishop和Kinra提出的理論框架［17］，對(duì)于不規(guī)則結(jié)構(gòu)，可將結(jié)構(gòu)劃分為若干個(gè)規(guī)則子區(qū)域，計(jì)算每個(gè)子區(qū)域的能量耗散，然后疊加得到總的耗散能。因此，部分覆蓋雙層梁的熱彈性阻尼解析解可通過以下方法計(jì)算：

2 理論分析與討論

2.1 與已有實(shí)驗(yàn)對(duì)比及理論驗(yàn)證

Enderling等［18］利用聚焦離子束工藝將金屬Pt分別沉積在13 μm×5 μm的SiC和Si懸臂梁自由端，Pt的厚度為0.5～2.6 μm，以實(shí)現(xiàn)微諧振器較大幅度的頻率調(diào)諧，圖2為掃描電鏡圖。Enderling等［18］發(fā)現(xiàn)自由端淀積Pt后諧振頻率與品質(zhì)因子均出現(xiàn)下降，且下降幅度隨Pt層的厚度增加而變大。

為驗(yàn)證本解有效性，以圖2對(duì)象為例，同時(shí)用有限元法做對(duì)比。微梁上表面受到簡(jiǎn)諧力F0sin（ωt）激勵(lì)，以彎曲振動(dòng)模式工作。

在ANSYS中，使用8節(jié)點(diǎn)Plane223單元，進(jìn)行熱?固耦合的諧響應(yīng)分析，利用能量法即可得到有限元法的熱彈性阻尼值。圖3給出了有限元模型及振型和溫度場(chǎng)云圖。表1列出了本文用到的材料參數(shù)。

圖4給出了利用有限元法和解析解與實(shí)驗(yàn)值對(duì)比的熱彈性阻尼的計(jì)算結(jié)果。QPt和QnoPt分別表示有Pt鍍膜和無Pt鍍膜的品質(zhì)因子。隨著Pt層厚度增加，品質(zhì)因子的降低幅度逐漸增加，且?guī)缀醭示€性趨勢(shì)。同時(shí)，解析解與有限元數(shù)值解吻合很好，與實(shí)驗(yàn)結(jié)果趨勢(shì)保持一致。

此外，完全覆蓋雙層微梁諧振器是部分覆蓋雙層梁的特例。鍍層與基底層長(zhǎng)度相同時(shí)，即L1=L2，當(dāng)前熱彈性阻尼解析解可以退化為完全覆蓋雙層梁的解。由于推導(dǎo)過程簡(jiǎn)單，此處不再贅述。綜上所述，當(dāng)前解析解是可靠有效的。

2.2 鍍層厚度、長(zhǎng)度和位置對(duì)熱彈性阻尼的影響

為方便分析，對(duì)尺寸參數(shù)進(jìn)行無參化表示，hr表示鍍層與基底層的厚度比，hr=h2/h1；Lr表示長(zhǎng)度比，Lr=L2/L1；Pr表示位置比，Pr=a/L1；Tr表示基底層的長(zhǎng)度與總厚度之比，Tr=L1/h。本文以下分析均基于Tr=40，h1=10 μm的假設(shè)。

微機(jī)械諧振器通常工作在第一階固有頻率（基頻）附近。圖5繪制了基頻下長(zhǎng)度比Lr從0～1的Al/Si懸臂梁中熱彈性阻尼曲線。鍍層常用兩種布置方案：（a） Pr=0，即從固定端延長(zhǎng)到自由端；（b） Pr=1-Lr，從自由端延伸到固定端。厚度比為hr=0.01，0.05和0.1。如圖5（a）所示，隨著Lr從0增加到0.6，熱彈性阻尼值急劇增加。然而，當(dāng)0≤Lr≤0.6時(shí)，圖5（b）中的熱彈性阻尼值幾乎保持不變。當(dāng)Lr≥0.6時(shí)，其值迅速增大。Pr=0.6時(shí)Lr=1的熱彈性阻尼值比Lr=0.6增大206%。此外，隨著hr增大，其值隨之增大。

圖6給出了基頻下位置比Pr從0增加到1-Lr的Al/Si懸臂梁中熱彈性阻尼值，長(zhǎng)度比為L(zhǎng)r=0.6。可看出，隨著Pr增大，即鍍層從夾緊端（Pr=0）向自由端移動(dòng)（Pr=0.4），阻尼值幾乎呈線性下降。Al鍍層在夾緊端處阻尼值比在自由端處高253%。

此外，圖5和6均給出有限元結(jié)果，可看出本解析解與有限元數(shù)值解吻合很好，證明了其有效性。更重要的是，圖5和圖6提示我們：布置在懸臂梁自由端（Pr=1-Lr）的金屬鍍層長(zhǎng)度在Lr≤0.6范圍內(nèi)不會(huì)引起熱彈性阻尼的顯著增大。

2.3 鍍層長(zhǎng)度和位置對(duì)SiO2/Si雙德拜峰的影響

頻率譜是實(shí)驗(yàn)中確定微諧振器品質(zhì)因子的關(guān)鍵依據(jù)。Nourmohammadi等［11］發(fā)現(xiàn)在h1/h2=1的SiO2/Si完全覆蓋雙層梁的熱彈性頻率譜中有兩個(gè)顯著的阻尼峰。圖7給出了長(zhǎng)度比Lr從0增長(zhǎng)到1的SiO2/Si懸臂梁熱彈性阻尼頻率譜：（a） Pr=0，（b） Pr=1-Lr。

從圖7（a）可看出，當(dāng)Lr≥0.2時(shí)，頻率譜中便出現(xiàn)了明顯的雙峰。但隨著Lr增大，高頻峰值沒有增大，反而降低。在Lr=1時(shí)，高頻峰值最小，說明SiO2對(duì)熱彈性阻尼有明顯的擬制作用，但在低頻處引入了一個(gè)額外峰。圖7（b）顯示，當(dāng)SiO2層覆蓋基底夾緊區(qū)域時(shí)（Lr≥0.8），才出現(xiàn)雙峰現(xiàn)象。

通過對(duì)比發(fā)現(xiàn)，相同長(zhǎng)度比Lr時(shí)，采用圖7（b）鍍層布置方案的熱彈性阻尼峰值要低于圖7（a）。因此，為降低熱彈性阻尼，應(yīng)優(yōu)先把鍍層布置在懸臂梁基底自由端位置。

單層梁的弛豫時(shí)間τ=h2/（π2α），在頻率f=1/τ=（π2α）/h2時(shí)，熱彈性阻尼達(dá)到峰值［4］。故推斷：雙峰現(xiàn)象是由于SiO2和Si的熱擴(kuò)散率α相差過大導(dǎo)致的。從表1可知SiO2和Si的熱擴(kuò)散率之比為αr=α1/α2=117。為解釋該推斷，有必要研究鍍層與基層之間不可逆熱流的耗散機(jī)理。

圖8繪制了單層Si梁、單層SiO2梁和雙層SiO2/Si （hr=1）梁在臨界頻率處沿厚度方向耦合溫度場(chǎng)虛部函數(shù)曲線。圖8（a）顯示了低頻峰對(duì)應(yīng)頻率處的溫度場(chǎng)虛部。從圖8（a）可知，單層Si梁和SiO2梁中均存在顯著的溫度梯度。假設(shè)SiO2/Si雙層梁的界面熱接觸是完美的，界面熱阻為0。由于αr=117，Si層的熱擴(kuò)散率遠(yuǎn)大于SiO2層。在SiO2層的臨界頻率處，SiO2層的能量耗散達(dá)到峰值，溫度梯度驅(qū)動(dòng)SiO2層的熱流擴(kuò)散到Si層，從而引起在低頻范圍內(nèi)出現(xiàn)峰值。

圖8（b）顯示了高頻峰對(duì)應(yīng)頻率處的溫度場(chǎng)虛部。由于高臨界頻率遠(yuǎn)大于低臨界頻率，單層SiO2梁處于絕熱態(tài)，能量耗散極小，溫度增量幾乎為0。但單層Si梁在臨界頻率處耗散最大，達(dá)到峰值。至于雙層SiO2/Si梁，由于SiO2層的熱擴(kuò)散率極低，因此在高頻處抑制了Si層的熱流耗散。因此，熱彈性阻尼在高臨界頻率處的峰值顯著減小。

3 結(jié) 論

本文提出了考慮沿厚度和長(zhǎng)度方向熱傳導(dǎo)的非完全覆蓋雙層微梁諧振器熱彈性阻尼解析解。根據(jù)綜合分析，簡(jiǎn)要概括出以下結(jié)論：

（1）該解析解可退化為完全覆蓋雙層梁熱彈性阻尼解，并且與有限元數(shù)值解吻合很好。

（2）為減少熱彈性阻尼，金屬鍍層應(yīng)遠(yuǎn)離梁夾緊端。對(duì)于懸臂梁，應(yīng)將鍍層放置在自由端而不是夾緊端。鍍層與基層的長(zhǎng)度比不超過0.6為宜。

（3）鍍層SiO2對(duì)熱彈性阻尼頻率譜具有顯著的擬制作用，但在低頻段引入一個(gè)額外低頻峰。雙峰現(xiàn)象是SiO2和Si的熱擴(kuò)散率α相差過大引起的。

參考文獻(xiàn)

1Bhugra H， Piazza P. Piezoelectric MEMS Resonators［M］. New York， NY， USA： Springer International Publishing， 2017.

2Mohanty P， Harrington D A， Ekinci K L， et al. Intrinsic dissipation in high-frequency micromechanical resonators［J］. Physical Review B， 2002， 66（8）： 085416.

3張文明，閆寒，彭志科，等.微納機(jī)械諧振器能量耗散機(jī)理研究進(jìn)展［J］.科學(xué)通報(bào)，2017，62（19）：2077-2093.

Zhang W M， Yan H， Peng Z K， et al. Research progress on energy dissipation mechanisms in micro-and nano-mechanical resonators［J］. Chinese Science Bulletin， 2017， 62（19）： 2077-2093.

4Lifshitz R， Roukes M L. Thermoelastic damping in micro-and nanomechanical systems［J］. Physical Review B， 2000， 61（8）： 5600-5609.

5Kumar S， Haque M A. Reduction of thermo-elastic damping with a secondary elastic field［J］. Journal of Sound and Vibration， 2008， 318（3）： 423-427.

6Prabhakar S， Vengallatore S. Thermoelastic damping in hollow and slotted microresonators［J］. Journal of Microelectromechanical Systems， 2009， 18（3）： 725-735.

7馬航空， 周晨陽， 李世榮. Mindlin矩形微板的熱彈性阻尼解析解［J］. 力學(xué)學(xué)報(bào)， 2020， 52（5）： 1383-1393

Ma H K， Zhou C Y， Li S R. Analytical soltion of thermoelastic damping in rectangular Mindlin micro plates［J］. Chinese Journal of Theoretical and Applied Mechanics， 2020， 52（5）： 1383-1393

8Kumar H， Mukhopadhyay S. Thermoelastic damping analysis in microbeam resonators based on Moore-Gibson-Thompson generalized thermoelasticity theory［J］. Acta Mechanica， 2020，231： 3003-3015.

9Gu B D， He T H， Ma Y B. Thermoelastic damping analysis in micro-beam resonators considering nonlocal strain gradient based on dual-phase-lag model［J］. International Journal of Heat and Mass Transfer， 2021， 180： 121771.

10Gad?el?Hak M. The MEMS Handbook［M］. Boca Raton： CRC Press， 2002

11Nourmohammadi Z， Prabhakar S， Vengallatore S. Thermoelastic damping in layered microresonators： critical frequencies， peak values， and rule of mixture［J］. Journal of Microelectromechanical Systems， 2013， 22（3）： 747-754.

12左萬里，黃家瀚.雙層矩形微板諧振器件中熱彈性阻尼機(jī)理研究［J］.傳感技術(shù)學(xué)報(bào)，2019，32（1）：50-56.

Zuo W L， Huang J H.Research on the thermoelastic damping mechanism in the double-layer rectangular microplate resonator device［J］. Journal of Transducer Technology， 2019， 32（1）： 50-56

13Yang L F， Li P， Fang Y M， et al. Thermoelastic damping in bilayer microbeam resonators with two-dimensional heat conduction［J］. International Journal of Mechanical Sciences， 2020， 167： 105245.

14Sandberg R， M?lhave K， Boisen A， et al. Effect of gold coating on the Q-factor of a resonant cantilever［J］. Journal of Micromechanics and Microengineering， 2005， 15（12）： 2249-2253

15Singh K V， Li G， Pang S S. Free vibration and physical parameter identification of non-uniform composite beams［J］. Composite Structures， 2006， 74（1）： 37-50.

16Nowick A S， Berry B S， Katz J L. Anelastic Relaxation in Crystalline Solids［M］. Elsevier Inc.， 1972.

17Bishop J E， Kinra V K. Elastothermodynamic damping in laminated composites［J］. International Journal of Solids and Structures， 1997， 34（9）： 1075-1092.

18Enderling S， Hedley J， Jiang L， et al. Characterization of frequency tuning using focused ion beam platinum deposition［J］. Journal of Micromechanics and Microengineering， 2006， 17（2）： 213-219.

Abstract Thermoelastic damping （TED） is one of the key coefficients that determines the upper limit of the quality factor of the micromechanical resonator. In the past， the TED model was only suitable for fully covered multi?layer microbeam structures. However， due to the manufacturing process and actual functional requirements， complex structural forms represented by incompletely covered double?layer beams are more common. Based on Fourier's law of heat conduction， this paper derives an analytical model for the TED of the partially covered bilayer microbeam resonators with the heat conduction along the length and thickness directions. Meanwhile， numerical methods and experimental data are used to verify the effectiveness of the model. The effects of the thickness， length and position of the coating on the TED is also analyzed. The results provide a new idea for reducing TED by optimizing the length and position of the coating.

Keywords micro-electro-mechanical-system （MEMS） resonators; thermoelastic damping; quality factor; thermoelasticity theory