高 展，徐小宇，閆 帥，呂鵬飛，任卓翔

（1.中國科學(xué)院微電子研究所 北京100029；2.中國科學(xué)院大學(xué) 北京100049；3.三維及納米集成電路設(shè)計自動化技術(shù)北京市重點實驗室 北京100029；4.巴黎六大L2E實驗室 巴黎75005）

基于離散幾何法的多物理場耦合問題研究

高 展1，2，徐小宇1，3，閆 帥1，3，呂鵬飛1，3，任卓翔1，4

（1.中國科學(xué)院微電子研究所 北京100029；2.中國科學(xué)院大學(xué) 北京100049；3.三維及納米集成電路設(shè)計自動化技術(shù)北京市重點實驗室 北京100029；4.巴黎六大L2E實驗室 巴黎75005）

Tonti圖揭示并高度概括了多種物理場所共同具備的數(shù)學(xué)結(jié)構(gòu)，為求解不同物理場提供相同的拓撲算子以構(gòu)建剛度矩陣。離散幾何法是基于離散空間中互為正交的原始和對偶網(wǎng)格，利用其幾何尺度直接導(dǎo)出本構(gòu)矩陣，從而快速建立代數(shù)方程的數(shù)值解法，它的形式結(jié)構(gòu)直觀簡單，便于計算。本文根據(jù)Tonti圖將離散幾何法應(yīng)用于三維集成電路中TSV陣列的電-熱-力場的耦合問題，并將計算結(jié)果與基于有限元法的商業(yè)軟件COMSOL對比。實驗結(jié)果表明：在多物理耦合問題上，離散幾何法精度可靠，從而有望應(yīng)用推廣至電子設(shè)計自動化工具之中。

離散幾何法；電-熱-力耦合問題；多物理；Tonti圖；電子設(shè)計自動化

隨著集成電路工藝尺寸達到18nm，甚至更小的技術(shù)代，電路中多種物理場的耦合作用越來越明顯，尤其是電-熱-力場的耦合。具體來說，就是電路密度的增大，使得功耗密度隨之增加，電路內(nèi)部溫度升高，進而產(chǎn)生熱應(yīng)力，導(dǎo)致器件與互連線產(chǎn)生形變，進而影響電路的可靠性，甚至使其失效。因此開發(fā)一種能高效處理多物理場耦合問題的算法來分析熱失效、熱應(yīng)力可靠性問題變得尤為重要[1]。從技術(shù)基礎(chǔ)及實現(xiàn)角度來看，就是為電子設(shè)計自動化（Electronic Design Automation，EDA）工具提供其所需要的耦合建模技術(shù)和場求解技術(shù)。

求解場的本質(zhì)就是求解不同物理場的偏微分方程（Partial Differential Equations，PDE）。因?qū)嶋H電路的復(fù)雜結(jié)構(gòu)使得解析求解難以實現(xiàn)，人們只能尋求數(shù)值解法，在誤差允許的范圍內(nèi)，用近似解代替精確解。當(dāng)前應(yīng)用于多物理場耦合問題的數(shù)值算法主要為有限元法（Finite Element Method，F(xiàn)EM），它具有精確度高、適應(yīng)性強、能處理復(fù)雜幾何模型、材料及非線性問題等優(yōu)點[2]。但是，有限元方法需要計算插值函數(shù)的微分與積分，當(dāng)采用矢量位作為未知變量時，受限于多連通域內(nèi)的多值問題，必須在區(qū)域間建立聯(lián)接（link），求解形式就比較復(fù)雜，自由度（未知量個數(shù)）也比較大[3-4]。離散幾何法（DiscreteGeometricMethod，DGM）則采用相互正交的原始網(wǎng)格與對偶網(wǎng)格之間的幾何尺度（例如長度、面積）比值來構(gòu)造單元材料矩陣，構(gòu)造形式非常簡單統(tǒng)一，同時，該材料矩陣具有對角陣的特點，對于采用Yee格式的高頻電磁場時域分析十分有利[5-6]，可能為集成電路內(nèi)部場的時域仿真提供高效的手段。

1 離散幾何法算法及實施

1.1 Tonti圖

E.Tonti等人根據(jù)對多種物理定律以及物理量的分析，提出了形式統(tǒng)一、適用于互補的離散幾何空間、表征多種物理定律的Tonti圖。其中，穩(wěn)恒電流場，熱傳導(dǎo)場和彈性力學(xué)場的Tonti圖如圖1所示[12]。根據(jù)圖中所標的原始循環(huán)路徑（primal cycle），表征穩(wěn)恒電流場的整體線性方程為：

其中，未知量v為原始網(wǎng)格節(jié)點上的標量電勢，G為離散空間的關(guān)聯(lián)矩陣，Mσ為電導(dǎo)率相關(guān)的本構(gòu)矩陣。

表征穩(wěn)態(tài)熱傳導(dǎo)的方程為：

其中，未知量T為原始網(wǎng)格節(jié)點上的溫度值，G與（2）式一致，是相同的關(guān)聯(lián)矩陣，Mλ為熱傳導(dǎo)系數(shù)相關(guān)的本構(gòu)矩陣。

表征彈性熱應(yīng)力場的線性方程為：

其中，未知量u為原始網(wǎng)格節(jié)點的位移；關(guān)聯(lián)矩陣G3的下標3表示，每個節(jié)點上有3個未知量，即分別沿x、y、z3個方向的位移量；ME為本構(gòu)矩陣，表示為ME=ADP，其中A表示對應(yīng)力進行積分的面矩陣，D是彈性矩陣，只取決于彈性體材料的彈性模量E與泊松比v，P是位移梯度矩陣的對稱部分[18]。

圖1 Tonti圖

1.2 關(guān)聯(lián)矩陣

對偏微分方程實施數(shù)值分析時，通常需要對空間進行一定形式的離散，同時即產(chǎn)生相互關(guān)聯(lián)的原始（primal）網(wǎng)格及對偶（dual）網(wǎng)格。對三維空間進行網(wǎng)格剖分時，通常采取Delaunay四面體剖分格式，因為它可以適應(yīng)任何形狀的幾何體，相對于長方體和棱柱單元具有更廣的適用性。為了更加便于理解，首先給出二維空間的三角形剖分及其對偶網(wǎng)格系統(tǒng)示意圖，如圖2（a）所示，由此再推及三維非結(jié)構(gòu)化四面體單元網(wǎng)絡(luò)（圖2（b））。在互為對偶關(guān)系的兩套網(wǎng)格之間，原始網(wǎng)格的外心是對偶網(wǎng)格的頂點，原始網(wǎng)格的線或者面的中垂線為對偶網(wǎng)格的棱邊。從原始網(wǎng)格到對偶網(wǎng)格的映射是點到體、線到面、面到線、體到點的映射變換[7]?；谖⒎謳缀卫碚摰腄GM方法，采用Hodge變換（或即映射），建立了物理量原始序列（施加在原始網(wǎng)格上的物理量）與對偶序列（施加在對偶網(wǎng)格上的物理量）各自內(nèi)在的以及互相之間的相關(guān)關(guān)系。

根據(jù)定義在對偶網(wǎng)格上的微分形，物理場方程可以在原始網(wǎng)格和對偶網(wǎng)格上以代數(shù)形式表達出來。連續(xù)偏微分方程中的梯度、旋度和散度算子對映離散空間中的關(guān)聯(lián)矩陣G（點到棱，grad）、C（棱到面，curl）、D（面到體，div）。同理，在對偶網(wǎng)格上，定義關(guān)聯(lián)矩陣根據(jù)其拓撲結(jié)構(gòu)的對偶關(guān)系[7-8]有：

基于同一套網(wǎng)格信息的不同物理場方程完全可以使用相同的關(guān)聯(lián)矩陣，不同的地方僅在于本構(gòu)矩陣。

圖2 原始網(wǎng)格與對偶網(wǎng)格相互關(guān)聯(lián)示意圖，對偶網(wǎng)格的節(jié)點采用原始網(wǎng)格的外心構(gòu)成

1.3 本構(gòu)矩陣

物理變量根據(jù)各自的性質(zhì)，分別定義在原始和對偶網(wǎng)格的點（如標量電位）、線（如矢量點位，如梯度）、面（如電通量，如旋度）上[9-10]。聯(lián)接定義在原始和對偶網(wǎng)格上的物理量的方程稱為本構(gòu)方程，本構(gòu)方程離散形式的系數(shù)矩陣稱為本構(gòu)矩陣?；贒GM方法，本構(gòu)矩陣通過互相正交的原始網(wǎng)格與對偶網(wǎng)格系統(tǒng)中的長度、面積等幾何尺寸以及材料的物理特性即可計算得出，它們均是對角陣。在三維空間，兩者的維度分別等于原始網(wǎng)格上的棱邊數(shù)量及面片數(shù)量[11]。

以穩(wěn)恒電流場的本構(gòu)矩陣為例，其元素是定義在原始網(wǎng)格每條棱上的電導(dǎo)，如圖3所示，相當(dāng)于該棱相鄰的所有單元在該棱上的電導(dǎo)的并聯(lián)。四面體ijks在棱ij上的電導(dǎo)可表達成對偶網(wǎng)格中的面片的幾何尺度（即面積）與原始網(wǎng)格上的棱的幾何尺度（即長度）之間比值，它對該棱邊上總電導(dǎo)的貢獻等于：

圖3 單元網(wǎng)格棱上電導(dǎo)求解示意圖

故而可以得到單元本構(gòu)矩陣：

方程單元矩陣為KPσ=GTMPeσG，將每個單元矩陣中的元素組裝起來就構(gòu)成整體矩陣，從而計算出未知量v。

穩(wěn)態(tài)熱傳導(dǎo)可以類比穩(wěn)恒電流場，電導(dǎo)網(wǎng)絡(luò)對應(yīng)熱阻網(wǎng)絡(luò)，即電導(dǎo)率對應(yīng)熱導(dǎo)率。同理求得熱場的本構(gòu)矩陣。

對于應(yīng)力場來說，應(yīng)力和應(yīng)變均為張量，為了獲得本構(gòu)矩陣，需要做一些中間處理，相對位移矢量h與機械應(yīng)變量通過如下定律聯(lián)接：

其中，P是位移梯度矩陣的對稱部分。單元內(nèi)應(yīng)力-應(yīng)變關(guān)系為：

其中，D為彈性矩陣。

其中，為溫變的楊氏模量，為泊松比。最后，應(yīng)力和面積力的柯西關(guān)系為。綜上所示，力學(xué)場的本構(gòu)矩陣。

1.4 基于DGM的電-熱-力耦合分析

電流場產(chǎn)生焦耳熱，焦耳熱可以成為溫度場場源；同時，溫度場會引起熱應(yīng)變（力），作為應(yīng)力場的場源。另一方面，溫度場通過模型材料的電阻率溫變特性而影響電場分布。換句話說，三個物理場通過彼此提供物理場場源或改變材料特性的方式進行強耦合。對于隨溫度變化不大的材料參數(shù)，如材料密度、熱容、泊松比等，通常取其室溫時的值；需考慮溫變的材料參數(shù)包括電導(dǎo)率σ（T）、熱導(dǎo)率λ（T）、楊氏模量E（T）和熱膨脹系數(shù)α（T）等。

首先，對于電-熱耦合，焦耳熱公式為

其次，對于熱-力耦合，熱量載荷作為外力源表示為：

fTH=其中，T為達到穩(wěn)態(tài)平衡時的溫度，T0為物體初始溫度，αTH為材料的熱膨脹系數(shù)。整理后，得到電-熱-力三種物理場的耦合方程組為：

聯(lián)立求解，即可得到各個節(jié)點上的標量電勢、溫度，以及形變，進一步則可以推導(dǎo)出其他諸如焦耳損失、應(yīng)力分布等物理特性[13]。

2 算法驗證

鑒于三維集成電路中的熱管理問題非常突出，利用典型的三維集成電路結(jié)構(gòu)——硅通孔（Though silicon via，TSV）陣列作為算例，研究其電-熱-力耦合問題，以便對DGM做初步驗證。熱源來自通孔內(nèi)部電流產(chǎn)生的焦耳熱。模型被剖分為100 787個四面體（節(jié)點個數(shù)為 18 990），總的未知量個數(shù)為93174，COMSOL軟件和DGM程序使用同一套網(wǎng)格信息。通孔上銅襯墊加0.05V電壓，底部的銅襯墊固定無位移。計算的模型和結(jié)果如圖4所示。從圖中可知，陣列的最大應(yīng)力值出現(xiàn)在通孔與二氧化硅交界面處。由于不同TSV發(fā)生相互作用，使得在相鄰?fù)字g的熱應(yīng)力也較密集。對相同算例，采用多物理耦合仿真軟件COMSOL（基于有限元法）進行對比驗證，對應(yīng)的結(jié)果非常吻合，計算得到的結(jié)論與文獻[16]也相一致，DGM的可靠性得到了驗證。表征三個物理場的參數(shù)的計算結(jié)果如表1所示。從表1的對比結(jié)果可知，文中所述的DGM結(jié)果與先進的多物理場仿真軟件COMSOL具有相同的精度。

圖4 （a）單個TSV結(jié)構(gòu)示意圖1至5的結(jié)構(gòu)材料分別是銅（Cu）、二氧化硅（SiO2）、硅（Si）、二氧化硅（SiO2）、銅（Cu）。（b）3×3TSV陣列應(yīng)力示意圖，硅襯底和上下的銅襯墊已省略（c）3×3TSV陣列沿Z軸形變量的俯視圖（d）沿（c）中黑線上應(yīng)力大小曲線圖

表1 3×3TSV陣列多物理場參數(shù)COMSOL和DGM計算結(jié)果對比

3 結(jié)束語

文中探討了求解物理場的離散幾何法，將該方法成功地用于三維模型的穩(wěn)態(tài)電-熱-力耦合問題之中。相對于有限元法使用Whitney棱單元構(gòu)建本構(gòu)矩陣[17-18]，離散幾何法直接使用網(wǎng)格的幾何尺度來獲得本構(gòu)矩陣，計算形式簡單直觀，且計算結(jié)果與仿真結(jié)果的精度一致。

文中所述的方法在考慮時變項后可以擴展到瞬態(tài)的多物理場耦合問題中去。利用Tonti圖的對偶循環(huán)路徑，及對偶空間的能量互補特性，求解基于原始和對偶路徑的兩套方程，能在粗網(wǎng)格條件下得出精確解[7]。

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

[2]謝德鑫，楊仕友.工程電磁場數(shù)值分析與綜合[M].北京：機械工業(yè)出版社，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，F(xiàn)reschi 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，F(xiàn)reschi 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]王勖成 ，邵敏.有限單元法基本原理和數(shù)值方法[M].2版.北京:清華大學(xué)出版社，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，F(xiàn)rance）

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

國家自然科學(xué)基金項目（61574167）

高 展（1991—），女，河南平頂山人，碩士研究生。研究方向：電磁及多物理耦合仿真，IC設(shè)計。