多環(huán)諧振微陀螺性能參數(shù)的有限元計(jì)算*

盧浩琳，張衛(wèi)平，谷留濤，馮 軍

（1.上海交通大學(xué) 微米納米加工技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室·上?！?00240；2.上海交通大學(xué) 電子信息與電氣工程學(xué)院 微納電子學(xué)系·上海·200240）

0 引 言

陀螺儀作為一種慣性傳感器，用于檢測(cè)物體在慣性坐標(biāo)系中轉(zhuǎn)動(dòng)的角度或角速度，是導(dǎo)航或制導(dǎo)系統(tǒng)的關(guān)鍵組成部分，在軍事和民用領(lǐng)域均發(fā)揮著重要作用[1]。相較于三維陀螺儀[2-4]，基于微機(jī)電系統(tǒng)（Micro-Electro-Mechanical System，MEMS）技術(shù)制造的多環(huán)諧振微陀螺儀具有體積小、功耗低、成本低、可批量生產(chǎn)等優(yōu)勢(shì)[5-6]，因此受到了廣泛關(guān)注。但是對(duì)于多環(huán)諧振微陀螺而言，由于結(jié)構(gòu)較復(fù)雜，不易求解，其解析需要借助有限元方法（Finite Element Method，F(xiàn)EM）。

有限元方法的核心思想是離散化，即將一個(gè)連續(xù)體離散為有限個(gè)單元，通過求解各單元上的解，從而得到整個(gè)連續(xù)體的近似解[7]。有限元方法自提出以來，經(jīng)過幾十年的發(fā)展，已經(jīng)成為工程應(yīng)用中常用的一種數(shù)值計(jì)算方法。并且結(jié)合計(jì)算機(jī)技術(shù)開發(fā)的有限元軟件也不斷出現(xiàn)，例如 ANSYS、ABAQUS、COMSOL Multiphysics（簡(jiǎn)稱 COMSOL）等[8-11],使用這些有限元軟件可以很方便地求解一些復(fù)雜的大規(guī)模的工程問題。COMSOL 提供了許多預(yù)定義的物理場(chǎng)接口，具有強(qiáng)大的多物理場(chǎng)耦合計(jì)算功能。該軟件最早源于 MATLAB 的Toolbox工具箱，目前仍提供與 MATLAB 協(xié)同計(jì)算的接口（COMSOL with MATLAB）。此外，它還提供 Microsoft Excel、SolidWorks、AutoCAD 等軟件的接口，為了更方便地實(shí)現(xiàn)與這些軟件之間的數(shù)據(jù)交流，本文主要基于 COMSOL 進(jìn)行仿真模擬。

1 機(jī)械靈敏度與機(jī)械熱噪聲

1.1 機(jī)械靈敏度

機(jī)械靈敏度定義為檢測(cè)軸方向的振動(dòng)幅度與輸入角速度之比，它表示微陀螺對(duì)于輸入角速度的靈敏性，單位為m/（rad/s），其表達(dá)式為

（1）

其中，qA0=（QAF0）/kA，F(xiàn)0表示驅(qū)動(dòng)力的幅度，kA表示驅(qū)動(dòng)模態(tài)的振動(dòng)位移；Ag為陀螺儀的角度增益；ωA、ωB分別為驅(qū)動(dòng)模態(tài)和檢測(cè)模態(tài)的諧振頻率；QA、QB分別為驅(qū)動(dòng)模態(tài)和檢測(cè)模態(tài)的品質(zhì)因子。當(dāng)不存在頻率裂解（ωA=ωB）時(shí)，機(jī)械靈敏度可以達(dá)到最大[12]

（2）

其中，ω0表示不存在頻率裂解時(shí)的諧振頻率，滿足ω0=ωA=ωB；Q表示不存在頻率裂解時(shí)的品質(zhì)因子，滿足Q=QA=QB。

1.2 機(jī)械熱噪聲

微陀螺的噪聲包括由結(jié)構(gòu)本身帶來的機(jī)械熱噪聲和由測(cè)控電路造成的電學(xué)噪聲[12]，如式（3）所示

（3）

其中，Ωoverall表示微陀螺的整體噪聲；Ωmech表示機(jī)械熱噪聲；Ωelect表示電學(xué)噪聲。陀螺儀的結(jié)構(gòu)本身主要受機(jī)械熱噪聲的影響，其表達(dá)式為[13]

（4）

（5）

由式（1）、式（2）、式（4）和式（5）可知，要計(jì)算機(jī)械靈敏度和機(jī)械熱噪聲，需要求解陀螺儀的品質(zhì)因子、頻率裂解、角度增益和有效質(zhì)量。

2 模態(tài)分析

微陀螺模型由SolidWorks 建立，結(jié)構(gòu)參數(shù)標(biāo)注如圖1所示，包括諧振子半徑R，中央錨點(diǎn)半徑r，諧振環(huán)寬度Wr，輻條寬度Ws，輻條長(zhǎng)度Ls，諧振子厚度H，電容間隙d0，諧振環(huán)數(shù)目N。

（a）俯視圖

（b）側(cè)視圖圖1 多環(huán)諧振微陀螺的結(jié)構(gòu)參數(shù)Fig.1 The structural parameters of disk resonator gyroscope

本文研究了有限元仿真平臺(tái)的搭建，對(duì)模型進(jìn)行了一定簡(jiǎn)化，基于如下兩點(diǎn)假設(shè)：

1）微陀螺為全對(duì)稱結(jié)構(gòu)，并且不存在加工誤差，各個(gè)結(jié)構(gòu)參數(shù)的取值如表1所示；

2）微陀螺的材料為各向同性的硅材料，其材料參數(shù)的取值如表2所示。

本文在沒有特殊說明的情況下，都是基于以上兩點(diǎn)假設(shè)進(jìn)行計(jì)算的。

表1 用于有限元計(jì)算的結(jié)構(gòu)參數(shù)取值

表2 用于有限元計(jì)算的材料參數(shù)取值

通過模態(tài)分析，可找出微陀螺的四波腹工作模態(tài)及其諧振頻率。從數(shù)值求解的角度來看，模態(tài)分析本質(zhì)上是求解矩陣的特征值。根據(jù)實(shí)際有限元模型的求解需求，選擇合適的特征值求解方法可以有效地提高求解過程的收斂性和運(yùn)算效率。對(duì)微陀螺進(jìn)行模態(tài)分析需要添加特征頻率研究，四波腹模態(tài)的諧振頻率約為22kHz，因此可將特征頻率的搜索基準(zhǔn)值設(shè)置為附近的18kHz，所需特征頻率數(shù)設(shè)置為4個(gè)，并從基準(zhǔn)值開始往更大的頻率進(jìn)行搜索。由于多環(huán)微陀螺具有高度對(duì)稱性，使用三角形與掃掠相結(jié)合的網(wǎng)格劃分方法可以實(shí)現(xiàn)高精度和高速度的計(jì)算。此外，相較于Z軸方向而言，XY平面內(nèi)的網(wǎng)格細(xì)化程度對(duì)結(jié)果的影響更明顯。本文將優(yōu)先采用自由三角形與掃掠相結(jié)合的網(wǎng)格劃分方法，并且重點(diǎn)關(guān)注XY平面內(nèi)的網(wǎng)格細(xì)化程度。圖2所示為使用COMSOL進(jìn)行模態(tài)分析求解得到的n=2模態(tài)振型。

圖2 COSMOL 求解n=2模態(tài)的振型Fig.2 Vibration shapes of n=2 mode in COMSOL

3 有效質(zhì)量與角度增益的求解

完成模態(tài)分析之后,即可求解有效質(zhì)量和角度增益。COMSOL中沒有預(yù)設(shè)的求解器對(duì)有效質(zhì)量和角度增益進(jìn)行求解，因此需要找出有效質(zhì)量和角度增益的有限元計(jì)算公式。引入一個(gè)表征模態(tài)形狀的參量——振型函數(shù)。振型函數(shù)定義為該振型中某個(gè)單元的位移與所有單元的最大位移（即波腹點(diǎn)位移）之比[14]。對(duì)于驅(qū)動(dòng)模態(tài)而言，某一單元x、y、z三個(gè)方向的振型函數(shù)可分別表示為

（6）

其中，u1、v1、w1分別表示驅(qū)動(dòng)模態(tài)每個(gè)單元在x、y、z三個(gè)方向上的位移；q1表示驅(qū)動(dòng)模態(tài)的波腹點(diǎn)位移。同理，對(duì)于檢測(cè)模態(tài)而言，振型函數(shù)可表示為

（7）

其中，u2、v2、w2分別表示檢測(cè)模態(tài)每個(gè)單元在x、y、z三個(gè)方向上的位移；q2表示檢測(cè)模態(tài)的波腹點(diǎn)位移。

文獻(xiàn)[6]給出了有效質(zhì)量和角度增益的計(jì)算公式，分別為式（8）和式（9）。式（8）事實(shí)上是驅(qū)動(dòng)模態(tài)的有效質(zhì)量， 理想情況下檢測(cè)模態(tài)的有效質(zhì)量可由式（9）表示。當(dāng)考慮微陀螺的不對(duì)稱性時(shí)，求解檢測(cè)模態(tài)的有效質(zhì)量則需要將式（8）中的振型函數(shù)用式（7）替代

（8）

（9）

其中，V表示微陀螺的體積；n表示模態(tài)階數(shù)。定義科氏質(zhì)量γcor為

（10）

則角度增益與有效質(zhì)量以及科氏質(zhì)量的關(guān)系為

（11）

對(duì)于有限元模型，式（8）和式（10）可寫成如下離散形式

[φz1（i）]2}δV（i）

（12）

φx2（i）φy1（i）}δV（i）

（13）

其中，δV（i）表示第i個(gè)單元的體積；Ne表示單元的總數(shù)目；i=1,2，…，Ne表示第i個(gè)單元。要求出有效質(zhì)量和角度增益，應(yīng)該先提取每個(gè)單元的位移，然后求出其振型函數(shù)，再對(duì)振型函數(shù)進(jìn)行體積分求得有效質(zhì)量和角度增益。

根據(jù)上述分析，求解有效質(zhì)量和角度增益的關(guān)鍵在于提取各網(wǎng)格單元的位移、體積等數(shù)據(jù)。MATLAB 與COMSOL可以實(shí)現(xiàn)很好的協(xié)同計(jì)算，因此可以采用MATLAB提取COMSOL有限元模型的節(jié)點(diǎn)數(shù)據(jù)，并在MATLAB中對(duì)這些節(jié)點(diǎn)數(shù)據(jù)進(jìn)行處理和計(jì)算。本文采用COMSOL with MATLAB進(jìn)行計(jì)算。計(jì)算過程如下：

第一步，在COMSOL中進(jìn)行模態(tài)分析。由于前面已經(jīng)得到特征頻率較精確的值，因此在該值附近提取特征頻率，可以設(shè)置所需特征頻率數(shù)為2，得到的二階模態(tài)分別對(duì)應(yīng)驅(qū)動(dòng)模態(tài)和檢測(cè)模態(tài)。

第二步，通過MATLAB的mpheval命令提取節(jié)點(diǎn)數(shù)據(jù)，需要提取的節(jié)點(diǎn)數(shù)據(jù)包括x、y、z方向的位移以及體積。使用mpheval命令在 MATLAB 中存儲(chǔ)一個(gè)結(jié)構(gòu)體，該結(jié)構(gòu)體包括：d1為所計(jì)算表達(dá)式在每個(gè)節(jié)點(diǎn)處的結(jié)果（如位移、體積） ；p為節(jié)點(diǎn)的坐標(biāo)；t為與節(jié)點(diǎn)坐標(biāo)對(duì)應(yīng)的節(jié)點(diǎn)列索引，其每列對(duì)應(yīng)于這些節(jié)點(diǎn)所構(gòu)成的一個(gè)網(wǎng)格單元，因此t的列數(shù)就是網(wǎng)格單元的個(gè)數(shù)，t的行數(shù)是網(wǎng)格單元對(duì)應(yīng)的節(jié)點(diǎn)數(shù)；ve為每個(gè)節(jié)點(diǎn)的網(wǎng)格單元索引；unit為所計(jì)算表達(dá)式的單位。

第三步，采用節(jié)點(diǎn)平均法對(duì)節(jié)點(diǎn)數(shù)據(jù)進(jìn)行處理。由于MATLAB提取出來的是節(jié)點(diǎn)數(shù)據(jù)，并不是網(wǎng)格單元的數(shù)據(jù)，因此需要將節(jié)點(diǎn)數(shù)據(jù)轉(zhuǎn)化為單元數(shù)據(jù)。每個(gè)網(wǎng)格單元由多個(gè)節(jié)點(diǎn)構(gòu)成，通過求解這些節(jié)點(diǎn)數(shù)據(jù)的平均值可以得到對(duì)應(yīng)網(wǎng)格單元的數(shù)據(jù)。該步驟的一個(gè)關(guān)鍵點(diǎn)在于如何將節(jié)點(diǎn)與單元對(duì)應(yīng)起來。由于上述結(jié)構(gòu)體中的列索引t存儲(chǔ)了每個(gè)單元對(duì)應(yīng)的節(jié)點(diǎn)，因此只需要依次求解t中每一列所有節(jié)點(diǎn)數(shù)據(jù)的平均值，即可得到每個(gè)單元數(shù)據(jù)的值，其中節(jié)點(diǎn)數(shù)據(jù)的值可以由d1得到。

第四步，通過前面三個(gè)步驟已經(jīng)得到每個(gè)單元的數(shù)據(jù)（位移與體積），接下來只需要根據(jù)式（11）、式（12）和式（13）即可求解有效質(zhì)量和角度增益。

這里給出了三種求解meff的方法和兩種求解Ag的方法。

3.1 求解meff

方法1：提取每個(gè)單元在x、y、z方向的位移以及波腹點(diǎn)位移，求出每個(gè)單元的振型函數(shù)，然后根據(jù)式（12）即可求出meff。

方法2：將驅(qū)動(dòng)模態(tài)和檢測(cè)模態(tài)所有單元的總位移分別構(gòu)成列向量d1和d2，所有單元在x、y、z方向的位移分別構(gòu)成列向量u1、v1、w1和u2、v2、w2， 所有單元的體積構(gòu)成列向量δV。根據(jù)式（8）可得

（14）

其中，diag運(yùn)算表示以向量的形式返回矩陣上的對(duì)角線元素。

方法3：式（8）可以改寫為

meff=λmMreal

（15）

其中，Mreal定義為微陀螺的實(shí)際質(zhì)量。有效質(zhì)量的大小為實(shí)際質(zhì)量乘以一個(gè)系數(shù)λm，其表達(dá)式為

（16）

將λm定義為有效質(zhì)量系數(shù)，表示有效質(zhì)量占實(shí)際質(zhì)量的比例?？梢钥闯?它是一個(gè)與模態(tài)形狀有關(guān)的量，與波腹點(diǎn)位移相近的單元數(shù)越多則λm越大。通過該方法求解有效質(zhì)量,需要先求解微陀螺的實(shí)際質(zhì)量Mreal和總體積Vtotal。

3.2 求解Ag

方法1：提取每個(gè)單元在x、y、z方向的位移以及波腹點(diǎn)位移， 求出每個(gè)單元的振型函數(shù)， 然后根據(jù)式（6）、式（7）和式（8）即可求出Ag。

方法2：與求解meff的方法2類似，角度增益可以由式（17）求出

（17）

采用非結(jié)構(gòu)化的自由四面體單元進(jìn)行網(wǎng)格劃分，單元大小設(shè)置為極細(xì)化，計(jì)算結(jié)果匯總于表3。 求解Ag的兩種方法在結(jié)果上有些許差異， 但這種差異是由于計(jì)算機(jī)精度造成的誤差，在合理范圍內(nèi)。

當(dāng)采用自由三角形與掃掠相結(jié)合的網(wǎng)格劃分方法時(shí)，利用上述方法進(jìn)行計(jì)算，得到的結(jié)果為：meff=7.3469×10-7kg，Mreal=2.3232×10-6kg，Vtotal=9.9753×10-10m3。這與表3的結(jié)果有很大差異， 并且Vtotal的值與實(shí)際值不一致（COMSOL 與SolidWorks 中測(cè)量得到的實(shí)際模型體積為 0.3324mm3），這說明該計(jì)算結(jié)果不正確。通過分析發(fā)現(xiàn)，此時(shí)結(jié)構(gòu)體中的t仍然是4行的矩陣，t的列數(shù)與實(shí)際的網(wǎng)格單元數(shù)也不相等。這說明此時(shí)在 MATLAB 中，每個(gè)網(wǎng)格單元仍被視為四面體單元，這與實(shí)際情況是不相符的。事實(shí)上，當(dāng)采用部分結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分時(shí)，每個(gè)單元為三棱柱結(jié)構(gòu)，而并非四面體結(jié)構(gòu)。 因此COMSOL with MATLAB協(xié)同計(jì)算僅僅適用于自由四面體單元的非結(jié)構(gòu)化網(wǎng)格劃分。

通過前文分析可知，在相同的精度下，采用非結(jié)構(gòu)化網(wǎng)格的計(jì)算量和計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)超過部分結(jié)構(gòu)化網(wǎng)格。COMSOL with MATLAB協(xié)同計(jì)算只限于非結(jié)構(gòu)化網(wǎng)格劃分，將造成計(jì)算緩慢且內(nèi)存需求很大。若利用COMSOL單獨(dú)計(jì)算，可不受網(wǎng)格劃分方式的限制。因此，有必要研究單獨(dú)利用COMSOL計(jì)算有效質(zhì)量和角度增益的方法。

COMSOL提供了一個(gè)可提取模型中任意解的算子withsol，該算子既可用在計(jì)算過程中，也可以用在數(shù)據(jù)的后處理當(dāng)中，使用起來特別方便。例如，使用“withsol（′sol1′,expr,setind（lambda,1））” 命令,即可提取第一階模態(tài)的表達(dá)式，其中setind用于控制所提取數(shù)據(jù)的模態(tài)階數(shù)，expr表示所要提取的表達(dá)式（x、y、z方向的位移、 波腹點(diǎn)位移、 體積等）。 采用withsol算子設(shè)置如下幾個(gè)變量：波腹點(diǎn)位移、振型函數(shù)、有效質(zhì)量、角度增益。其中還需要定義兩個(gè)算子：最大值算子maxop1，用于提取波腹點(diǎn)位移；積分算子intop1， 用于進(jìn)行積分運(yùn)算。 同樣采用部分結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分， 計(jì)算結(jié)果為：meff=2.4489×10-7kg，Ag=0.3958。 該結(jié)果與前文基于COMSOL with MATLAB的計(jì)算結(jié)果基本一致，其差異是由于網(wǎng)格劃分方式不同造成的。

相較于COMSOL with MATLAB協(xié)同計(jì)算，采用COMSOL單獨(dú)計(jì)算不受網(wǎng)格劃分方式的限制，可以提高計(jì)算精度和計(jì)算速度。

4 品質(zhì)因子的求解

由于品質(zhì)因子可以由熱彈性品質(zhì)因子來表示，因此求解品質(zhì)因子等同于求解熱彈性阻尼。研究熱彈性阻尼需要同時(shí)添加力學(xué)場(chǎng)和熱學(xué)場(chǎng)，COMSOL中的熱彈性力學(xué)物理場(chǎng)耦合了固體力學(xué)與固體傳熱，可以用于熱彈性阻尼的分析。同時(shí)添加熱擾動(dòng)和特征頻率研究，初始溫度T0設(shè)置為室溫293.15K，熱量來源于熱彈性阻尼[15]。

COMSOL內(nèi)置耦合熱彈性方程，并將熱擾動(dòng)下的特征頻率用復(fù)數(shù)表示，利用前文所述的復(fù)頻率方法可以得出熱彈性品質(zhì)因子QTED，求解得到QTED的值為72418。同時(shí)，COMSOL還可以求得微陀螺儀中的溫度分布，如圖3所示，由于各處溫度分布不均勻，導(dǎo)致了振動(dòng)能量以熱流的形式耗散。

圖3 多環(huán)諧振微陀螺中的溫度分布Fig.3 Temperature distribution in the disk resonator gyroscope

5 結(jié) 論

本文主要介紹了多環(huán)諧振微陀螺性能參數(shù)的有限元計(jì)算方法。 在陀螺儀基本研究的基礎(chǔ)上，搭建了基于COMSOL 和 MATLAB 等軟件的有限元計(jì)算平臺(tái)，實(shí)現(xiàn)了對(duì)諧振頻率、有效質(zhì)量、角度增益和品質(zhì)因子等性能參數(shù)的有限元求解。并且這些性能參數(shù)既可由 COMSOL with MATLAB 協(xié)同計(jì)算，也可由 COMSOL 單獨(dú)計(jì)算。本文搭建的有限元計(jì)算平臺(tái)為研究微陀螺結(jié)構(gòu)的性能參數(shù)和優(yōu)化設(shè)計(jì)提供了有力的工具。