张 程，徐小宇，任卓翔

（1.中国科学院微电子研究所，北京 100029；2.中国科学院大学，北京 100049；3.索邦大学，巴黎75006）

集成电路加工工艺已经进入跨越3 nm 的时代，新结构、新材料的应用以及工艺涨落对于电路性能的影响也更加显著，也对互联寄生参数的提取提出了更高的要求。参数提取主要依赖麦克斯韦方程组的电磁场求解器，其中，广泛采用的方法是有限元法（FEM）[1-6]。总的来说，有限元法将节点标量电势、棱边矢量电位等作为未知量，对于网格，采用插值规则建立代数方程组，其中网格单元可以是不同几何形状[1-7]。所有有限元求解器的作用都针对的是性能提升，并在效率与精度之间进行权衡。近些年来，有限元法的能量互补和对偶性被广泛研究，该特性被用来加速求解收敛[3-4]。所谓的对偶有限元法，即节点有限元和棱边有限元，分别将原始网格上节点标量电位和棱边上的矢量电位作为未知量。电容参数提取是一类静电场问题，采用矢量电位时会出现多值问题，可通过在带电荷区域和多连通域之间建立“连接”来解决[1]。根据欧拉特性，棱边有限元总是会导致出现更多的自由度，因此在对偶网格上采用标量电势作为未知量，可以提供FEM 的另一种解决方案如下。

对偶网格由多边形（二维）及多面体（三维）单元组成，广义重心坐标系（GBCs）可被用来为这些单元提供近似插值函数[7-8]。基于GBCs的伽辽金格式的有限元可称为多边形或多面体有限元。近些年来，基于GBCs的有限元法获得了较多的关注并显示出新的前景[9-13]，且在计算机图形学、拓扑优化[14-15]、力学[9]、热学和电磁学[1,8]等领域已经有许多研究，然而该算法的高复杂度和耗时阻碍了应用推广。为了在三维对偶网格上更加有效地建立变量有限元，文中提出了一种在多棱柱单元上分片插值的方法，并通过算例给出三维电容参数提取结果的收敛性和精确性。

1 原始网格及对偶网格

给定一个Delaunay三角化网格，很容易获得它的关联对偶网格，分别称为原始网格MP及对偶网格MD，如图1 所示。通常情况下，采用基于外心的方式从原始网格获得对偶网格，此时两套网格互相正交。

图1 网格的二维及三维视图

虽然已有PolyMesher 等工具可用来生成多边形网格[16]，但为了从已建立的原始网格中产生对偶网格，应当注意不同区域的交界面。如图2 所示，由于界面的存在，多边形c3、c2被分割为若干个部分。尤其是MP中形态不佳的三角形会导致多边形c2和c5也处于形态不佳的情况，这时可以局部采用重心来构建对偶网格。

图2 网格以及区域界面

2 广义重心坐标系及插值

2.1 广义重心坐标系定义

有限元法中经常采用的重心坐标系仅适用于三角形、四面体等单纯形，对于更加广泛的多边形与多面体，广义重心坐标系（GBCs）可视为重心坐标在多边形和多面体上的延伸。GBCs 函数λi在多边形或多面体Ω内的某个点vi满足非负性、线性完备性、单位性、克罗内克函数特性等性质[1,8]。

GBCs 满足作为代数基的伽辽金逼近函数的充分必要条件，也有大量GBCs 可供选择。在FEM 框架下，需要计算λ与gradλ，并且对gradλ及gradλi、gradλj在单元上进行积分。对于GBCs 而言，该计算非常耗时，可以考虑采用近似插值技术。

2.2 分片插值广义重心坐标系

考虑到在三棱柱上的能量积分相对容易，可以考虑选取一个点将具有n条边的凸多棱柱分为多个三棱柱，其中电场和能量可以分开计算。这种分片插值算法可用来简化GBCs 和其梯度计算，并可以用于FEM 计算。

2.2.1 选取插值点

为了准确地表达多棱柱各节点之间的关联关系，首先需要定位插值点，根据经验可知，多胞体的重心或最大内切球的球心是较好的选择。考虑到可能会存在多个最大内切球，所以选择其中最靠近重心的球心。为简便起见，文中采用重心作为插值点，如图3（a）中的do，n边多棱柱可分为n个三棱柱，其中，do(1)与do(2)分别是do在顶面与底面上的投影，包含了do的横截面如图3（b）所示。设定凸多边形d1…dn的第i个顶点坐标为(xi,yi,zi)。

图3 多棱柱单元及其横截面图

2.2.2 确定插值点上的标量电势

一旦给出了插值点，施加在该点上的电势值即可被确定。一般采用能量最小原则来确定插值点上的标量电势。

如图3（a）所示，设多棱柱上的节点上的电势为，其中r=1,…,n；s=1,2。假定重心do上的待定电势为vo，每个子区域中的能量可根据公式计算得出，累加所有子区域的能量，就可以得到多棱柱单元的总能量We。

显然，We是vo的函数，可通过局部极值点的计算来查找vo以使局部能量最小化，也就是寻求对该点电势的偏导数为零的能量：

根据计算得到电势为：

其中，权重系数γi可表示为：

可以证明权重系数在每个单元内的和为1，即满足单元性。并且也可以得到：

因此，基于式（2）获得的值是最小值。

2.2.3 形函数

如图3（a）所示，考察落入子区域dodidi+1的d*点上的电势，可以推出单元dodidi+1上的形函数N(s)为：

其中，s取值1、2，分别表示上、下表面，且有：

关于棱柱单元d1...dn之中电势的连续性，根据式（5）～（7）可以看出，在插值点、顶点及连线上是0阶连续，在所有棱、面上是1 阶连续，在其他地方是无穷阶连续。

2.2.4 单元及总体刚度矩阵

根据局部区域能量最小的原理以及式（5）～（7），对于计算区域Ω上的偏微分方程问题，可以很容易得到n边多棱柱单元的系数矩阵元素。

这样就建立起了单元和总体刚度矩阵。最后，求解所得到稀疏矩阵，即可得到每一个节点上的标量电势。

3 数值分析算例

寄生参数提取是集成电路工业中的重要议题，它是典型的静电场问题。文中利用所提出的有限元法在对偶网格MD中提取8 导体结构的电容值。如图4 所示，C8被设为主导体并施加1 V 电压，其他导体作为环境导体并施加0 V 电压。

图4 三维多导体系统的电容参数提取算例结构示意图

图4 中描述的是8 导体结构，导体被相对介电常数为3.70 的材料（图中所示的尺寸非等比例）包围。

如图5 所示，计算区域经Delaunay 三角化离散形成原始网格MP，然后通过外心形成对偶网格MD。文中采用了标量的原始有限元（primal FEM，标记为pFEM）、矢量对偶有限元（dual FEM，标记为dFEM），以及所提出的多胞体（如多棱柱）有限元（polytope FEM，标记为tFEM）。

图5 网格单元

电势分布及能量的计算结果如图6 所示。由于对偶有限元的两种方法具有能量互补特征，其均值也被用来作为参考基准。针对原始网格上的算例#1 与算例#2，算例#1 具有9 126 个节点、16 150 个单元，算例#2 具有41 834 个节点、77 550 个单元，计算相关概要信息，如表1、表2 所示。对应情况下，指定的主导体C8的自电容与互电容结果如表3、表4 所示。

图6 不同有限元法计算获得的电势分布计算结果

表1 不同算例网格的概要信息

表2 不同算例有限元法计算的概要信息

表3 不同有限元法的算例#1电容参数提取结果

表4 不同有限元法的算例#2电容参数提取结果

随着网格加密的电容参数提取结果如图7 所示，结果表明，多胞体有限元法具有良好的收敛性和精度，与原始有限元法相比，虽然耗时略多，但明显具有更高的精度。

图7 导体C4与C8之间的电容

4 结论

与基于原始网格的传统标量有限元法及矢量有限元法不同，文中给出了一种基于对偶网格的标量有限元法，该方法采用广义重心坐标系进行插值函数的构建。这类方法近些年才重新兴起，并在包括电磁场分析的领域有所应用，但是广义重心坐标及相应的有限元法计算复杂度过高，文中针对三维多棱柱单元，提出一种基于分片插值的方法，以快速获得伽辽金格式的有限元法的形函数，并将其用于集成电路的参数提取。对数值的分析结果表明，该方法具有良好的精度，且随着网格的加密也呈现出稳定的收敛性。该方法容易推广到其他单元类型，并且对偶网格也很容易由原始网格导出。

从经验上来看，多胞体有限元法比原始网格上的标量有限元法的精度高，该方法实际上隐含了在多个细密的网格上进行形函数导出的过程。为了将多胞体有限元应用于实际工程，有必要对网格离散误差和有限元计算误差的估计进行理论研究。