基于离散几何法的多物理场耦合问题研究

高 展，徐小宇，闫 帅，吕鹏飞，任卓翔

（1.中国科学院微电子研究所 北京100029；2.中国科学院大学 北京100049；3.三维及纳米集成电路设计自动化技术北京市重点实验室 北京100029；4.巴黎六大L2E实验室 巴黎75005）

基于离散几何法的多物理场耦合问题研究

高 展1，2，徐小宇1，3，闫 帅1，3，吕鹏飞1，3，任卓翔1，4

（1.中国科学院微电子研究所 北京100029；2.中国科学院大学 北京100049；3.三维及纳米集成电路设计自动化技术北京市重点实验室 北京100029；4.巴黎六大L2E实验室 巴黎75005）

Tonti图揭示并高度概括了多种物理场所共同具备的数学结构，为求解不同物理场提供相同的拓扑算子以构建刚度矩阵。离散几何法是基于离散空间中互为正交的原始和对偶网格，利用其几何尺度直接导出本构矩阵，从而快速建立代数方程的数值解法，它的形式结构直观简单，便于计算。本文根据Tonti图将离散几何法应用于三维集成电路中TSV阵列的电-热-力场的耦合问题，并将计算结果与基于有限元法的商业软件COMSOL对比。实验结果表明：在多物理耦合问题上，离散几何法精度可靠，从而有望应用推广至电子设计自动化工具之中。

离散几何法；电-热-力耦合问题；多物理；Tonti图；电子设计自动化

随着集成电路工艺尺寸达到18nm，甚至更小的技术代，电路中多种物理场的耦合作用越来越明显，尤其是电-热-力场的耦合。具体来说，就是电路密度的增大，使得功耗密度随之增加，电路内部温度升高，进而产生热应力，导致器件与互连线产生形变，进而影响电路的可靠性，甚至使其失效。因此开发一种能高效处理多物理场耦合问题的算法来分析热失效、热应力可靠性问题变得尤为重要[1]。从技术基础及实现角度来看，就是为电子设计自动化（Electronic Design Automation，EDA）工具提供其所需要的耦合建模技术和场求解技术。

求解场的本质就是求解不同物理场的偏微分方程（Partial Differential Equations，PDE）。因实际电路的复杂结构使得解析求解难以实现，人们只能寻求数值解法，在误差允许的范围内，用近似解代替精确解。当前应用于多物理场耦合问题的数值算法主要为有限元法（Finite Element Method，FEM），它具有精确度高、适应性强、能处理复杂几何模型、材料及非线性问题等优点[2]。但是，有限元方法需要计算插值函数的微分与积分，当采用矢量位作为未知变量时，受限于多连通域内的多值问题，必须在区域间建立联接（link），求解形式就比较复杂，自由度（未知量个数）也比较大[3-4]。离散几何法（DiscreteGeometricMethod，DGM）则采用相互正交的原始网格与对偶网格之间的几何尺度（例如长度、面积）比值来构造单元材料矩阵，构造形式非常简单统一，同时，该材料矩阵具有对角阵的特点，对于采用Yee格式的高频电磁场时域分析十分有利[5-6]，可能为集成电路内部场的时域仿真提供高效的手段。

1 离散几何法算法及实施

1.1 Tonti图

E.Tonti等人根据对多种物理定律以及物理量的分析，提出了形式统一、适用于互补的离散几何空间、表征多种物理定律的Tonti图。其中，稳恒电流场，热传导场和弹性力学场的Tonti图如图1所示[12]。根据图中所标的原始循环路径（primal cycle），表征稳恒电流场的整体线性方程为：

其中，未知量v为原始网格节点上的标量电势，G为离散空间的关联矩阵，Mσ为电导率相关的本构矩阵。

表征稳态热传导的方程为：

其中，未知量T为原始网格节点上的温度值，G与（2）式一致，是相同的关联矩阵，Mλ为热传导系数相关的本构矩阵。

表征弹性热应力场的线性方程为：

其中，未知量u为原始网格节点的位移；关联矩阵G3的下标3表示，每个节点上有3个未知量，即分别沿x、y、z3个方向的位移量；ME为本构矩阵，表示为ME=ADP，其中A表示对应力进行积分的面矩阵，D是弹性矩阵，只取决于弹性体材料的弹性模量E与泊松比v，P是位移梯度矩阵的对称部分[18]。

图1 Tonti图

1.2 关联矩阵

对偏微分方程实施数值分析时，通常需要对空间进行一定形式的离散，同时即产生相互关联的原始（primal）网格及对偶（dual）网格。对三维空间进行网格剖分时，通常采取Delaunay四面体剖分格式，因为它可以适应任何形状的几何体，相对于长方体和棱柱单元具有更广的适用性。为了更加便于理解，首先给出二维空间的三角形剖分及其对偶网格系统示意图，如图2（a）所示，由此再推及三维非结构化四面体单元网络（图2（b））。在互为对偶关系的两套网格之间，原始网格的外心是对偶网格的顶点，原始网格的线或者面的中垂线为对偶网格的棱边。从原始网格到对偶网格的映射是点到体、线到面、面到线、体到点的映射变换[7]。基于微分几何理论的DGM方法，采用Hodge变换（或即映射），建立了物理量原始序列（施加在原始网格上的物理量）与对偶序列（施加在对偶网格上的物理量）各自内在的以及互相之间的相关关系。

根据定义在对偶网格上的微分形，物理场方程可以在原始网格和对偶网格上以代数形式表达出来。连续偏微分方程中的梯度、旋度和散度算子对映离散空间中的关联矩阵G（点到棱，grad）、C（棱到面，curl）、D（面到体，div）。同理，在对偶网格上，定义关联矩阵根据其拓扑结构的对偶关系[7-8]有：

基于同一套网格信息的不同物理场方程完全可以使用相同的关联矩阵，不同的地方仅在于本构矩阵。

图2 原始网格与对偶网格相互关联示意图，对偶网格的节点采用原始网格的外心构成

1.3 本构矩阵

物理变量根据各自的性质，分别定义在原始和对偶网格的点（如标量电位）、线（如矢量点位，如梯度）、面（如电通量，如旋度）上[9-10]。联接定义在原始和对偶网格上的物理量的方程称为本构方程，本构方程离散形式的系数矩阵称为本构矩阵。基于DGM方法，本构矩阵通过互相正交的原始网格与对偶网格系统中的长度、面积等几何尺寸以及材料的物理特性即可计算得出，它们均是对角阵。在三维空间，两者的维度分别等于原始网格上的棱边数量及面片数量[11]。

以稳恒电流场的本构矩阵为例，其元素是定义在原始网格每条棱上的电导，如图3所示，相当于该棱相邻的所有单元在该棱上的电导的并联。四面体ijks在棱ij上的电导可表达成对偶网格中的面片的几何尺度（即面积）与原始网格上的棱的几何尺度（即长度）之间比值，它对该棱边上总电导的贡献等于：

图3 单元网格棱上电导求解示意图

故而可以得到单元本构矩阵：

方程单元矩阵为KPσ=GTMPeσG，将每个单元矩阵中的元素组装起来就构成整体矩阵，从而计算出未知量v。

稳态热传导可以类比稳恒电流场，电导网络对应热阻网络，即电导率对应热导率。同理求得热场的本构矩阵。

对于应力场来说，应力和应变均为张量，为了获得本构矩阵，需要做一些中间处理，相对位移矢量h与机械应变量通过如下定律联接：

其中，P是位移梯度矩阵的对称部分。单元内应力-应变关系为：

其中，D为弹性矩阵。

其中，为温变的杨氏模量，为泊松比。最后，应力和面积力的柯西关系为。综上所示，力学场的本构矩阵。

1.4 基于DGM的电-热-力耦合分析

电流场产生焦耳热，焦耳热可以成为温度场场源；同时，温度场会引起热应变（力），作为应力场的场源。另一方面，温度场通过模型材料的电阻率温变特性而影响电场分布。换句话说，三个物理场通过彼此提供物理场场源或改变材料特性的方式进行强耦合。对于随温度变化不大的材料参数，如材料密度、热容、泊松比等，通常取其室温时的值；需考虑温变的材料参数包括电导率σ（T）、热导率λ（T）、杨氏模量E（T）和热膨胀系数α（T）等。

首先，对于电-热耦合，焦耳热公式为

其次，对于热-力耦合，热量载荷作为外力源表示为：

fTH=其中，T为达到稳态平衡时的温度，T0为物体初始温度，αTH为材料的热膨胀系数。整理后，得到电-热-力三种物理场的耦合方程组为：

联立求解，即可得到各个节点上的标量电势、温度，以及形变，进一步则可以推导出其他诸如焦耳损失、应力分布等物理特性[13]。

2 算法验证

鉴于三维集成电路中的热管理问题非常突出，利用典型的三维集成电路结构——硅通孔（Though silicon via，TSV）阵列作为算例，研究其电-热-力耦合问题，以便对DGM做初步验证。热源来自通孔内部电流产生的焦耳热。模型被剖分为100 787个四面体（节点个数为 18 990），总的未知量个数为93174，COMSOL软件和DGM程序使用同一套网格信息。通孔上铜衬垫加0.05V电压，底部的铜衬垫固定无位移。计算的模型和结果如图4所示。从图中可知，阵列的最大应力值出现在通孔与二氧化硅交界面处。由于不同TSV发生相互作用，使得在相邻通孔之间的热应力也较密集。对相同算例，采用多物理耦合仿真软件COMSOL（基于有限元法）进行对比验证，对应的结果非常吻合，计算得到的结论与文献[16]也相一致，DGM的可靠性得到了验证。表征三个物理场的参数的计算结果如表1所示。从表1的对比结果可知，文中所述的DGM结果与先进的多物理场仿真软件COMSOL具有相同的精度。

图4 （a）单个TSV结构示意图1至5的结构材料分别是铜（Cu）、二氧化硅（SiO2）、硅（Si）、二氧化硅（SiO2）、铜（Cu）。（b）3×3TSV阵列应力示意图，硅衬底和上下的铜衬垫已省略（c）3×3TSV阵列沿Z轴形变量的俯视图（d）沿（c）中黑线上应力大小曲线图

表1 3×3TSV阵列多物理场参数COMSOL和DGM计算结果对比

3 结束语

文中探讨了求解物理场的离散几何法，将该方法成功地用于三维模型的稳态电-热-力耦合问题之中。相对于有限元法使用Whitney棱单元构建本构矩阵[17-18]，离散几何法直接使用网格的几何尺度来获得本构矩阵，计算形式简单直观，且计算结果与仿真结果的精度一致。

文中所述的方法在考虑时变项后可以扩展到瞬态的多物理场耦合问题中去。利用Tonti图的对偶循环路径，及对偶空间的能量互补特性，求解基于原始和对偶路径的两套方程，能在粗网格条件下得出精确解[7]。

[1]International technology roadmap for semiconductors（ITRS）reports 2012[EB/OL].[2016-04-13] http://www.itrs.net.

[2]谢德鑫，杨仕友.工程电磁场数值分析与综合[M].北京：机械工业出版社，2009.

[3]REN Zhuo-xiang.2-D dual finite-element formulations for the fast extraction of circuit parameters in VLSI[J].IEEE Transactions on Magnetics，2003，39（3）:1590-1593.

[4]REN Zhuo-xiang.A 3D vector potential formulation using edge element for electrostatic-field computation[J].IEEE Transactions on Magnetics，1995，31（3）:1520-1523.

[5]REN Zhuo-xiang，XU Xiao-yu.Dual discrete geometric methods in terms of scalar potential on unstructured mesh in electrostatics[J].IEEE Transaction on Magnetics，2014，50（2）:37-40.

[6]BossavitA，Kettunen L.Yee-like schemes on staggered cellular grids:a synthesis between FIT and FEM approaches[J].IE EE Transactions on Magnetics，2000，36（4）:861-867.

[7]REN Zhuo-xiang.On the complementarity of dual formulations on dual meshes [J].IEEE Transactions on Magnetics，2009，45（3）:1284-1287.

[8]Geoffrey D.Flexible implementation of the finite element method applied to 3D coupled prob-lems considering convective effects [D].Leuven-Heverlee，Belgi:Katholieke Universiteit Leuven，2003.

[9]Gross P W，Kotiuga P R.Electromagnetic the-ory and computation:a topologicalapproach [M]. Cambridge:Cambridge Univ.Press，2004.

[10]Tonti E.The mathematical structure of classical and relativistic physics:a general classification diagram[M]，Switzerland:Birkh user，2013.

[11]Tonti E.Finite formulation of electromagnetic field [J].IEEE Transactions on Magnetics，2002，38（1）: 333-336.

[12]XU Xiao-yu，REN Zhuo-xiang，QU Hui，REN Dan. 3-D IC interconnect capacitance extraction using dual discrete geometric methods with prism elements [J].，IEEE Transaction on Very Large Scale Integration（VLSI）Systems，2015，2:1305-1309.

[13]Alotto P，Freschi F，Repetto M，et al.The cell method for electrical engineering and mul-tiphysics problem [M].Berlin:Springer-Verlag Berlin Heidelberg，2013.

[14]REN Dan，XU Xiao-yu，QU Hui，REN Zhuo-xiang，Two-dimensional parasitic capacitance extraction for integrated circuit with dual discrete geometric methods[J].Journal of Semiconductors，2015，36（4）:045008-1-7.

[15]Alotto P，Freschi F，Repetto M.Multiphy-sics problems via the cell method:the role of Tonti Diagrams[J].IEEE Transaction on Magnetics，2010，46（8）:2959-2962.

[16]MITRA J，JUNG M，RYU S K，et al.A fast simulation framework for full-chip thermo-mechanical stress and reliability analysis of through-siliconvia based 3D ICs[C]//IEEE Confererence on Electronic Components and Technology Conference，2011:746-753.

[17]Specogna R，Trevisan F，Discrete constitu-tive equations in geometric eddy-current formula-tion [J]，IEEE Transaction Magnetics，2005，41（4）: 1259-1263.

[18]王勖成 ，邵敏.有限单元法基本原理和数值方法[M].2版.北京:清华大学出版社，1997.

Study of multiphysics problem based on discrete geometric method

GAO Zhan1，2，XU Xiao-yu1，3，YAN Shuai1，3，LV Peng-fei1，3，REN Zhuo-xiang1，4
（1.Institute of Microelectronics of Chinese Academy of Sciences，Beijing 100029，China；2.University of Chinese Academy of Sciences，Beijing 100049，China；3.Beijing Key Laboratory of 3D&Nano IC Electronic Design Automatic Technologies，Beijing 100029，China；4.Sorbonne Universities，UPMC Univ Paris 06，UR2，L2E，Paris 75005，France）

The Tonti diagrams disclose a general common mathematical structure for several physical fields，and hence provide uniform topological operators to establish global stiffness matrices for various field solvers in differential form.The discrete geometric method implements the Hodge discretization of constitutive laws through the division of geometric dimensions of mutually orthogonal primal and dual meshes.Complying with Tonti diagrams，it serves as an elegant solution for multiphysics problems intuitively and practicably.The stationary electro-thermo-mechanical coupled problem， including building of the constitutive equations and treating of the coupling terms，is studied in this paper.A typical problem with regard to electro-thermo-mechanical stress induced mainly by resistive losses in the TSV structures is investigated.The proposed method offers accurate solutions as the commercial software COMSOL based on the finite element method，and is expected to applied in the utilities of electronic design automation.

discrete geometric method；electro-thermo-mechanical coupled problem；multiphysics；tonti diagram；electronic design automation（EDA）

TN402

：A

：1674－6236（2017）01-0166-05

2016-04-16稿件编号：201604170

国家自然科学基金项目（61574167）

高 展（1991—），女，河南平顶山人，硕士研究生。研究方向：电磁及多物理耦合仿真，IC设计。