宋明鑫 張碧星
(1 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)
(2 中國科學院大學 北京 100049)
近年來隨著半導體行業的快速發展,聲表面波(Surface acoustic wave, SAW)器件在智能傳感、氣體檢測以及移動通訊等領域取得了一系列卓越的進步[1?3],而快速、精確的模擬理論是支撐SAW 器件發展的關鍵因素。作為SAW 器件的一種新型數值模擬方法,級聯有限元技術(Hierarchical cascade technique-Finite element method,HCT-FEM)因其計算速度快、模擬精度高以及具備聲場計算能力等優勢在行業內引起了廣泛關注[4?5]。HCT-FEM 最早由Koskela 等[6?8]提出,其核心是利用SAW 器件中叉指換能器(Interdigital transducer, IDT)的周期性,從整體結構中分離出若干相互獨立的基本單元,通過單元邊界矩陣(Boundary matrix)的級聯實現有限長器件結構的模擬。在實際的SAW器件設計中,通過優化IDT的指條間距,激發出不同波長的聲表面波,調整SAW器件的諧振頻率。HCT-FEM 需要針對不同的波長分別建立相應的FEM 單元模型,頻繁的建模、分析以及數據傳輸過程會消耗大量的計算時間與內存資源。針對上述問題,本文提出一種頻率變化時的HCT-FEM快速模擬方法,用于提升SAW器件的模擬效率。首先分析HCT-FEM中基本單元的系統矩陣構成以及不同波長對應單元系統矩陣間的關系。接著引入參考單元模型與波長變化因子,建立單元系統矩陣與波長間的函數關系。利用參考單元的FEM 模型快速計算出不同波長對應的單元系統矩陣。進一步通過單元邊界矩陣的級聯實現SAW 器件頻域響應的計算與快速變頻效果。
SAW 器件中存在大量周期性排布的結構單元,依據HCT-FEM的基本原理[6],首先從有限長SAW器件結構中分離出若干基本單元。如圖1所示,當IDT 的指條間距、電極指寬、膜厚等參數保持不變時,單根指條單元(基本單元1)的不斷排布就構成了完整的SAW器件結構。圖1中基本單元2表示吸收層單元,以同樣的級聯方式疊加聲波吸收效果,實現有限長器件左右兩側的自由邊界狀態[9]。
根據有限元的基本原理,對上述基本單元分別建模分析。在壓電介質中,考慮模型網格劃分后的任一單元,其節點上的自由度包括節點位移ua與節點電勢?a,單元內任意位置(x,y,z)的狀態通過插值的方式獲得:

式(1)中,a表示單元的節點序號,n為節點個數,Na表示該節點對應的形函數。

圖1 同步單端SAW 諧振器的基本單元構成Fig.1 The basic units of synchronous one-port SAW resonator
對于模型網格中的任一單元,其剛度矩陣與質量矩陣為[10]

其中,積分區域V為網格單元的體積,ρ表示介質密度,C為彈性矩陣,由材料的彈性模量與泊松比決定。上標T代表矩陣的轉置,D為微分算子,在笛卡爾坐標下:

將上述單元剛度矩陣Ke、質量矩陣Me集總成模型的整體剛度矩陣K與整體質量矩陣M,構建不考慮阻尼的系統方程[7]:

式(4)中,ω= 2πf,表示頻域分析的角頻率,U為模型的整體自由度矢量,F表示激勵源與邊界條件。令A=(K ?ω2M),稱為基本單元的系統矩陣。
根據模型中節點坐標的位置,將U劃分成[UL,UI,UR,?] 四部分。其中,[UL,UI,UR]分別表示單元FEM 模型左邊界、內部、右邊界上所包含的節點位移與電勢自由度;?為電極區域的節點電勢。模型的系統矩陣A做相應的重排[7]:

其中,Aij是分塊處理后的各子矩陣,τL、τR是左右邊界上的節點應力與電荷值,q為電極電荷。由模型內部應力與電荷平衡狀態可知,UI對應的線性方程組等式為零。即

將式(6)代入式(5),消去內部自由度UI對應的方程式,獲得新的模型邊界方程組,其系數矩陣稱為邊界矩陣[7],即B矩陣:

其中,邊界矩陣的各元素由系統矩陣的各分塊子矩陣Aij決定,表達式為

通過邊界矩陣的級聯實現SAW 器件的快速模擬[11]。從以上HCT-FEM的基本算法流程可知,不同的聲表面波波長與系統矩陣A之間存在一一對應的關系。在實際器件的模擬中,需要根據不同的波長,建立相應的基本單元FEM 模型。構建出單元的系統矩陣后通過分塊矩陣求逆獲得最終的單元邊界矩陣。當波長優化調整時,頻繁的單元建模、分析以及數據傳輸過程降低了HCT-FEM 的模擬效率。
由上述有限元模擬方法可知, 系統矩陣A= (K ?ω2M)的每一個矩陣元都與波長有關,因此,當波長(頻率)改變時,需要對每一波長分別進行模擬和計算,在有限元計算中,模型的單元網格數很大,因而系統矩陣的維數也很大,導致聲表面波器件模擬的計算時間很長,在實際中并不方便。為此,本文提出了一種計算不同波長時的聲表面波器件變頻快速模擬方法。
IDT 的指條間距決定了基片中激發出的聲表面波波長,是SAW器件優化設計中的重要參數。在選定基片的壓電材料以后,不同指條間距的IDT 實現不同頻率的聲波諧振[12]。依據HCT-FEM 的基本原理,首先建立如圖2所示SAW 器件的準三維單電極模型,并進行網格劃分。其結構參數如表1所示。其中,λ為波長,其余參數均為λ的倍數。圖2中頂端黃色區域表示金屬電極,藍色區域為壓電基片,基片末端設置完美匹配層(Perfect matched layer,PML)。

圖2 SAW 諧振器單根電極的準三維模型Fig.2 Quasi three dimensional model of a single electrode of SAW resonator

表1 單根電極單元的幾何結構參數Table 1 Geometric parameter setting in reference unit
在聲表面波器件結構中,聲表面波的能量分布區域是與波長λ成正比的,從而,在有限元模型中,模型中每一網格的邊長(長、寬、高)以及集總成總模型的邊長都與波長λ成正比,使得每一網格體積和集總成總體積都與波長的立方成正比。由有限元模擬方法可知,當波長變化時,公式(2)中的形函數與波長無關,只與模型中網格的相對位置有關。
為明確不同波長對應的單元系統矩陣間的關系,首先考慮單根電極模型網格劃分后的任意一個長方體單元R1,其長、寬、高分別用a、b、c表示,如圖3所示。當波長增加α倍(由λ增大到αλ)時,單元R1的長、寬、高在y、x、z三個方向上同比例拉伸α倍,得到單元R2。相比于R1,長方體單元R2的體積V擴大α3倍。根據公式(2)可知,節點形函數Na與介質密度ρ保持不變,積分區域擴大了α3倍,R2的單元質量矩陣Me是R1單元質量矩陣的α3倍。從而,集總成的整體質量矩陣也增為原來的α3倍。

圖3 長方體單元在各方向上同比例的拉伸示意圖Fig.3 Schematic diagram of cuboid unit stretched in the same proportion in all directions
對于式(2)中的剛度矩陣,首先考慮微分算子D的作用。在有限元計算中,微分算子D是由相鄰節點間的函數差(這個差值與波長無關)除以單元長度得到的,由于單元長度與波長成正比,使得微分算子D與波長成反比,剛度矩陣中包含了兩次微分算子。從而,在公式(2)中,被積函數與波長的平方成反比,由于單元網格體積與波長的立方成正比,使得剛度矩陣與波長成正比。
因此,在固定的網格劃分模式下,當波長增加α倍時,對于模型中的每一網格單元,其質量矩陣和剛度矩陣分別增加α3和α倍,進行裝配后,模型的整體質量矩陣M與剛度矩陣K中每個元素分別擴大α3與α倍數。
根據以上分析,選取一個參考波長λ0,對應的頻域分析角頻率為ω0,這時,模型的整體剛度矩陣與質量矩陣分別為K(λ0)和M(λ0),其系統矩陣為

在第i次模擬中,波長為λi,其剛度矩陣與質量矩陣分別為K(λi)和M(λi),從而系統矩陣為

由前面的分析,可知

其中α=λi/λ0,因此,

通過式(12),建立起了不同波長所對應的系統矩陣之間的簡易變換關系。只需要選取一個參考波長λ0,并得到整體的系統矩陣A(λ0),然后,利用式(12),即可得到任意波長λi時的系統矩陣A(λi)。在特定的頻域區間內,細化分析頻點,獲得完整的SAW器件變頻后的導納響應。當波長改變時,這種計算方法大大減少了有限元的模擬計算時間,提高了HCT-FEM 用于SAW 器件模擬與變頻優化時的計算效率。
結合上述快速計算方法,利用HCT-FEM 實現同步單端SAW 諧振器的模擬與快速變頻過程。諧振器的襯底材料設置為42?YX鉭酸鋰,其材料參數由Kovacs 等[13]提出的方法計算獲得。電極材料為金屬鋁,密度為2700 kg/m3,泊松比為0.35,楊氏模量為70×109Pa。單電極參考單元模型的結構參數如表2所示,其中IDT 根數為121,左右反射柵為40 根。

表2 單電極參考單元的結構參數Table 2 Structure parameters of single electrode reference unit
依據表2中的結構參數,首先建立波長為4 μm的參考單元FEM 模型,提取模型的質量矩陣與剛度矩陣,并構建出其系統矩陣A0。當波長調整為4.01 μm 時,波長的變化因子α為1.0025。利用式(12)直接獲得新的系統矩陣A1,該過程只包含簡單的矩陣乘法,計算速度很快。接著利用單元級聯算法,在950~1080 MHz之間,分別計算兩種波長的SAW 諧振器頻域響應。圖4顯示了該有限長SAW諧振器的導納計算結果。圖中紅色線為HCT-FEM基于參考單元計算得到的諧振器導納幅值與實部;黑色線是利用快速變頻方法將波長為4.01 μm時的計算結果。從圖4中可知,參考模型的諧振頻率與反諧振頻率為987.7 MHz 與1027.2 MHz。當波長調整為4.01 μm時,導納曲線的位置整體左移,諧振頻點與反諧振頻點出現在985.1 MHz 與1024.6 MHz,頻移的大小為2.6 MHz。兩者諧振點處幅值的微小差異來自于電極阻抗的變化。

圖4 SAW 諧振器的導納計算結果Fig.4 Calculation result of SAW resonator admittance
為進一步驗證該快速變頻模擬方法的可靠性,利用有限元方法直接建立波長為4.01 μm時的基本單元,構建其系統矩陣后,將單元邊界矩陣的級聯結果與快速變頻結果對比,如圖5所示。從圖5中可以看出,兩者吻合較好。

圖5 波長為4.01 μm 時的直接建模與快速變頻結果對比Fig.5 Comparison between the direct modeling and fast frequency change at IDT wavelength of 4.01 μm
在上述變頻模擬中,單根電極單元的FEM 模型自由度總數為3174,基本單元系統矩陣A的維度是3174×3174,消去模型內部與前后側面上的自由度后,基本單元B矩陣的維度是512×512。快速變頻方法優化前后的時間對比如表3所示。其中,優化前的傳統方法完成單元模型一次FEM 分析的時間大致為10 s,單次數據傳輸的時間為14.2 s。而在快速變頻方法中,用一次矩陣數乘代替了建模、分析、數據傳輸過程,單次矩陣乘數的時間為0.02 s,模擬效率得到顯著提升。所用的計算機主板配置為B360M,CPU 處理器為6 核i5-9400F,主頻率2.9 GHz。運行內存為16.0 GB,內存工作的主頻率為2666 MHz。

表3 快速變頻方法優化前后的時間對比Table 3 Comparison of time before and after optimization by fast frequency conversion method
在SAW 器件參數設計中,往往存在著幾百次的優化過程,即需要對幾百個波長值分別進行模擬。利用本文快速計算方法,只需建立一個參考模型,該模型對應著一個波長值,然后,利用式(12)可獲得任意波長的單元系統矩陣,該過程只涉及矩陣的簡易乘法,避免了需要對每一波長值進行建模和相關處理過程,從而節省了大量的時間,模擬效率大大提升。
本文研究了一種基于級聯有限元的聲表面波器件變頻快速模擬方法。在SAW 器件的級聯模擬算法中需要根據不同的波長分別建立FEM 單元模型,本文引入參考波長,建立了單元系統矩陣與波長變化的簡易函數關系。利用參考單元的FEM 模型實現了任意波長對應單元系統矩陣的快速計算,提升了HCT-FEM 的快速變頻模擬效率,在SAW 的優化設計中具有較強的實用價值。仿真結果驗證了該方法模擬的可靠性。在SAW 器件的結構參數優化中,存在波長增加,同時控制膜厚不變的情況,后續研究會圍繞這種情況開展進一步研究。