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

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

（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陣列的電-熱-力場的耦合問題，并將計算結果與基于有限元法的商業(yè)軟件COMSOL對比。實驗結果表明：在多物理耦合問題上，離散幾何法精度可靠，從而有望應用推廣至電子設計自動化工具之中。

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

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

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

1 離散幾何法算法及實施

1.1 Tonti圖

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

其中，未知量v為原始網格節(jié)點上的標量電勢，G為離散空間的關聯矩陣，Mσ為電導率相關的本構矩陣。

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

其中，未知量T為原始網格節(jié)點上的溫度值，G與（2）式一致，是相同的關聯矩陣，Mλ為熱傳導系數相關的本構矩陣。

表征彈性熱應力場的線性方程為：

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

圖1 Tonti圖

1.2 關聯矩陣

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

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

基于同一套網格信息的不同物理場方程完全可以使用相同的關聯矩陣，不同的地方僅在于本構矩陣。

圖2 原始網格與對偶網格相互關聯示意圖，對偶網格的節(jié)點采用原始網格的外心構成

1.3 本構矩陣

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

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

圖3 單元網格棱上電導求解示意圖

故而可以得到單元本構矩陣：

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

穩(wěn)態(tài)熱傳導可以類比穩(wěn)恒電流場，電導網絡對應熱阻網絡，即電導率對應熱導率。同理求得熱場的本構矩陣。

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

其中，P是位移梯度矩陣的對稱部分。單元內應力-應變關系為：

其中，D為彈性矩陣。

其中，為溫變的楊氏模量，為泊松比。最后，應力和面積力的柯西關系為。綜上所示，力學場的本構矩陣。

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

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

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

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

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

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

2 算法驗證

鑒于三維集成電路中的熱管理問題非常突出，利用典型的三維集成電路結構——硅通孔（Though silicon via，TSV）陣列作為算例，研究其電-熱-力耦合問題，以便對DGM做初步驗證。熱源來自通孔內部電流產生的焦耳熱。模型被剖分為100 787個四面體（節(jié)點個數為 18 990），總的未知量個數為93174，COMSOL軟件和DGM程序使用同一套網格信息。通孔上銅襯墊加0.05V電壓，底部的銅襯墊固定無位移。計算的模型和結果如圖4所示。從圖中可知，陣列的最大應力值出現在通孔與二氧化硅交界面處。由于不同TSV發(fā)生相互作用，使得在相鄰通孔之間的熱應力也較密集。對相同算例，采用多物理耦合仿真軟件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 結束語

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

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

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

[2]謝德鑫，楊仕友.工程電磁場數值分析與綜合[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，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設計。