多環諧振微陀螺性能參數的有限元計算*

盧浩琳，張衛平，谷留濤，馮 軍

（1.上海交通大學 微米納米加工技術國家級重點實驗室·上海·200240；2.上海交通大學 電子信息與電氣工程學院 微納電子學系·上海·200240）

0 引 言

陀螺儀作為一種慣性傳感器，用于檢測物體在慣性坐標系中轉動的角度或角速度，是導航或制導系統的關鍵組成部分，在軍事和民用領域均發揮著重要作用[1]。相較于三維陀螺儀[2-4]，基于微機電系統（Micro-Electro-Mechanical System，MEMS）技術制造的多環諧振微陀螺儀具有體積小、功耗低、成本低、可批量生產等優勢[5-6]，因此受到了廣泛關注。但是對于多環諧振微陀螺而言，由于結構較復雜，不易求解，其解析需要借助有限元方法（Finite Element Method，FEM）。

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

1 機械靈敏度與機械熱噪聲

1.1 機械靈敏度

機械靈敏度定義為檢測軸方向的振動幅度與輸入角速度之比，它表示微陀螺對于輸入角速度的靈敏性，單位為m/（rad/s），其表達式為

（1）

其中，qA0=（QAF0）/kA，F0表示驅動力的幅度，kA表示驅動模態的振動位移；Ag為陀螺儀的角度增益；ωA、ωB分別為驅動模態和檢測模態的諧振頻率；QA、QB分別為驅動模態和檢測模態的品質因子。當不存在頻率裂解（ωA=ωB）時，機械靈敏度可以達到最大[12]

（2）

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

1.2 機械熱噪聲

微陀螺的噪聲包括由結構本身帶來的機械熱噪聲和由測控電路造成的電學噪聲[12]，如式（3）所示

（3）

其中，Ωoverall表示微陀螺的整體噪聲；Ωmech表示機械熱噪聲；Ωelect表示電學噪聲。陀螺儀的結構本身主要受機械熱噪聲的影響，其表達式為[13]

（4）

（5）

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

2 模態分析

微陀螺模型由SolidWorks 建立，結構參數標注如圖1所示，包括諧振子半徑R，中央錨點半徑r，諧振環寬度Wr，輻條寬度Ws，輻條長度Ls，諧振子厚度H，電容間隙d0，諧振環數目N。

（a）俯視圖

（b）側視圖圖1 多環諧振微陀螺的結構參數Fig.1 The structural parameters of disk resonator gyroscope

本文研究了有限元仿真平臺的搭建，對模型進行了一定簡化，基于如下兩點假設：

1）微陀螺為全對稱結構，并且不存在加工誤差，各個結構參數的取值如表1所示；

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

本文在沒有特殊說明的情況下，都是基于以上兩點假設進行計算的。

表1 用于有限元計算的結構參數取值

表2 用于有限元計算的材料參數取值

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

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

3 有效質量與角度增益的求解

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

（6）

其中，u1、v1、w1分別表示驅動模態每個單元在x、y、z三個方向上的位移；q1表示驅動模態的波腹點位移。同理，對于檢測模態而言，振型函數可表示為

（7）

其中，u2、v2、w2分別表示檢測模態每個單元在x、y、z三個方向上的位移；q2表示檢測模態的波腹點位移。

文獻[6]給出了有效質量和角度增益的計算公式，分別為式（8）和式（9）。式（8）事實上是驅動模態的有效質量， 理想情況下檢測模態的有效質量可由式（9）表示。當考慮微陀螺的不對稱性時，求解檢測模態的有效質量則需要將式（8）中的振型函數用式（7）替代

（8）

（9）

其中，V表示微陀螺的體積；n表示模態階數。定義科氏質量γcor為

（10）

則角度增益與有效質量以及科氏質量的關系為

（11）

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

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

（12）

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

（13）

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

根據上述分析，求解有效質量和角度增益的關鍵在于提取各網格單元的位移、體積等數據。MATLAB 與COMSOL可以實現很好的協同計算，因此可以采用MATLAB提取COMSOL有限元模型的節點數據，并在MATLAB中對這些節點數據進行處理和計算。本文采用COMSOL with MATLAB進行計算。計算過程如下：

第一步，在COMSOL中進行模態分析。由于前面已經得到特征頻率較精確的值，因此在該值附近提取特征頻率，可以設置所需特征頻率數為2，得到的二階模態分別對應驅動模態和檢測模態。

第二步，通過MATLAB的mpheval命令提取節點數據，需要提取的節點數據包括x、y、z方向的位移以及體積。使用mpheval命令在 MATLAB 中存儲一個結構體，該結構體包括：d1為所計算表達式在每個節點處的結果（如位移、體積） ；p為節點的坐標；t為與節點坐標對應的節點列索引，其每列對應于這些節點所構成的一個網格單元，因此t的列數就是網格單元的個數，t的行數是網格單元對應的節點數；ve為每個節點的網格單元索引；unit為所計算表達式的單位。

第三步，采用節點平均法對節點數據進行處理。由于MATLAB提取出來的是節點數據，并不是網格單元的數據，因此需要將節點數據轉化為單元數據。每個網格單元由多個節點構成，通過求解這些節點數據的平均值可以得到對應網格單元的數據。該步驟的一個關鍵點在于如何將節點與單元對應起來。由于上述結構體中的列索引t存儲了每個單元對應的節點，因此只需要依次求解t中每一列所有節點數據的平均值，即可得到每個單元數據的值，其中節點數據的值可以由d1得到。

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

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

3.1 求解meff

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

方法2：將驅動模態和檢測模態所有單元的總位移分別構成列向量d1和d2，所有單元在x、y、z方向的位移分別構成列向量u1、v1、w1和u2、v2、w2， 所有單元的體積構成列向量δV。根據式（8）可得

（14）

其中，diag運算表示以向量的形式返回矩陣上的對角線元素。

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

meff=λmMreal

（15）

其中，Mreal定義為微陀螺的實際質量。有效質量的大小為實際質量乘以一個系數λm，其表達式為

（16）

將λm定義為有效質量系數，表示有效質量占實際質量的比例。可以看出,它是一個與模態形狀有關的量，與波腹點位移相近的單元數越多則λm越大。通過該方法求解有效質量,需要先求解微陀螺的實際質量Mreal和總體積Vtotal。

3.2 求解Ag

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

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

（17）

采用非結構化的自由四面體單元進行網格劃分，單元大小設置為極細化，計算結果匯總于表3。 求解Ag的兩種方法在結果上有些許差異， 但這種差異是由于計算機精度造成的誤差，在合理范圍內。

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

通過前文分析可知，在相同的精度下，采用非結構化網格的計算量和計算時間遠遠超過部分結構化網格。COMSOL with MATLAB協同計算只限于非結構化網格劃分，將造成計算緩慢且內存需求很大。若利用COMSOL單獨計算，可不受網格劃分方式的限制。因此，有必要研究單獨利用COMSOL計算有效質量和角度增益的方法。

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

相較于COMSOL with MATLAB協同計算，采用COMSOL單獨計算不受網格劃分方式的限制，可以提高計算精度和計算速度。

4 品質因子的求解

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

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

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

5 結 論

本文主要介紹了多環諧振微陀螺性能參數的有限元計算方法。 在陀螺儀基本研究的基礎上，搭建了基于COMSOL 和 MATLAB 等軟件的有限元計算平臺，實現了對諧振頻率、有效質量、角度增益和品質因子等性能參數的有限元求解。并且這些性能參數既可由 COMSOL with MATLAB 協同計算，也可由 COMSOL 單獨計算。本文搭建的有限元計算平臺為研究微陀螺結構的性能參數和優化設計提供了有力的工具。