宋明鑫 张碧星

(1 中国科学院声学研究所 声场声信息国家重点实验室 北京 100190)

(2 中国科学院大学 北京 100049)

0 引言

近年来随着半导体行业的快速发展，声表面波(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 器件频域响应的计算与快速变频效果。

1 单元系统矩阵

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 的模拟效率。

2 变频时的快速模拟方法

由上述有限元模拟方法可知， 系统矩阵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 器件模拟与变频优化时的计算效率。

3 仿真验证

结合上述快速计算方法，利用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)可获得任意波长的单元系统矩阵，该过程只涉及矩阵的简易乘法，避免了需要对每一波长值进行建模和相关处理过程，从而节省了大量的时间，模拟效率大大提升。

4 结论

本文研究了一种基于级联有限元的声表面波器件变频快速模拟方法。在SAW 器件的级联模拟算法中需要根据不同的波长分别建立FEM 单元模型，本文引入参考波长，建立了单元系统矩阵与波长变化的简易函数关系。利用参考单元的FEM 模型实现了任意波长对应单元系统矩阵的快速计算，提升了HCT-FEM 的快速变频模拟效率，在SAW 的优化设计中具有较强的实用价值。仿真结果验证了该方法模拟的可靠性。在SAW 器件的结构参数优化中，存在波长增加，同时控制膜厚不变的情况，后续研究会围绕这种情况开展进一步研究。