基于有限元模型的三維集成電路熱分析

杜文雄，唐普英，王舒冰，張思曼

（電子科技大學(xué) 光電信息學(xué)院，四川 成都 610054）

基于有限元模型的三維集成電路熱分析

杜文雄，唐普英，王舒冰，張思曼

（電子科技大學(xué) 光電信息學(xué)院，四川 成都 610054）

基于有限元理論，對三維集成電路（3D-IC）進行了建模和仿真，研究了不同模型的熱分布和計算復(fù)雜度。通過Gmsh軟件創(chuàng)建3D-IC模型并生成網(wǎng)格化文件。利用Matlab軟件提取有限元參數(shù),獲到模型的剛度矩陣。用層次矩陣(Hierarchical matrix,H-matrix)表示剛度矩陣，得到了不同模型剛度矩陣的求逆所消耗的存儲空間和運算時間。結(jié)果表明：隨著模型剛度矩陣行列數(shù)目的增加，所需要的運算時間和存儲空間呈現(xiàn)線性變化關(guān)系。

Gmsh；H-matrix；三維集成電路；有限元熱分析；剛度矩陣

自1947年第一個半導(dǎo)體晶體管的發(fā)明，到1958年采用硅平面工藝的集成電路的誕生，半導(dǎo)體集成電路一直遵循著摩爾定律高速發(fā)展。隨著集成電路技術(shù)的持續(xù)發(fā)展，在單一芯片上集成更多的資源已面臨著重要的挑戰(zhàn)。三維集成電路(3D-IC)技術(shù)被視為延續(xù)摩爾定律的有效技術(shù)之一，它可以改善整個集成電路系統(tǒng)的功耗，提高器件的集成度。硅孔道(Through Silicon Via,TSV)是3D-IC中堆疊芯片實現(xiàn)互聯(lián)的一種新的解決方案。

1)集成電路熱分析的意義

一個集成電路實際上是由電學(xué)子系統(tǒng)和熱學(xué)子系統(tǒng)共同組成，兩個系統(tǒng)互相耦合：電學(xué)元件的功耗作為熱學(xué)網(wǎng)絡(luò)的熱源，而熱學(xué)網(wǎng)絡(luò)中不同溫度值作為參數(shù)會影響電學(xué)系統(tǒng)中元器件及其性能。隨著元器件密度，工作速度以及集成電路規(guī)模的不斷增加，集成電路的能耗密度也越來越大，導(dǎo)致片上溫度越來越高。如果不考慮集成電路的熱分布問題，電路設(shè)計和布局布線就有可能不合理而導(dǎo)致集成電路溫度過高以及器件之間溫度梯度過大[1]。因此，為保證電路性能和可靠性的要求，有必要在集成電路制造之前對其進行系統(tǒng)的熱分析。

2)3 D-IC的優(yōu)點和技術(shù)挑戰(zhàn)

3D-IC的優(yōu)點：①提高封裝密度。多層器件重疊結(jié)構(gòu)可成倍提高芯片集成度。②提高電路工作速度。重疊結(jié)構(gòu)使單元連線縮短，并使并行信號處理成為可能，從而實現(xiàn)電路的高速操作。③可實現(xiàn)新型多功能器件及電路系統(tǒng)。

3D-IC面臨技術(shù)的挑戰(zhàn)：①散熱問題：由于電路系統(tǒng)擁有了更高的集成度，熱功耗也隨之提升，表面積與體積比隨之下降，傳統(tǒng)的平面散熱技術(shù)不再能滿足立體集成電路的散熱要求。②測試問題：傳統(tǒng)測試技術(shù)只針對單層系統(tǒng)，而未提供針對多層芯片集成之后的整體系統(tǒng)測試技術(shù)。③制作工藝復(fù)雜，成本較高[2-3]。

3)本文的主要工作和意義

國外在3D-IC的熱模擬方面已進行了一些研究工作。國內(nèi)在這方面的報道很少，多是綜述性的文章，蘇州大學(xué)的李文石等人對二維和三維集成電路的熱阻進行了計算[4]。胡輝勇等人提出了一種基于 Si和SiGe材料的三維CMOS集成電路,與 Si三維 CMOS結(jié)構(gòu)相比，該電路結(jié)構(gòu)具有明顯的速度優(yōu)勢[5]。

本文通過對含TSV的3D-IC進行建模和仿真，研究不同模型的熱分布、存儲空間復(fù)雜度和運算時間復(fù)雜度。仿真結(jié)果對3D-IC的設(shè)計有指導(dǎo)作用。

1 基于有限元理論的熱模型

有限元分析是使用有限元方法來分析靜態(tài)或動態(tài)的物理物體或物理系統(tǒng)。在這種方法中一個物體或系統(tǒng)被分解為由有限多個相互聯(lián)結(jié)的、簡單、獨立的點組成的幾何模型。實際的物理模型中推導(dǎo)出來的平衡方程式被使用到每個點上，由此產(chǎn)生了一個方程組。這個方程組可以用線性代數(shù)的方法來求解。有限元分析的基本流程：1)節(jié)點編號和單元劃分。2)計算各單元的單元剛度方程。3)組裝各單元剛度方程。4)處理邊界條件并求解[6]。

有限元熱分析基于能量守恒原理的熱平衡方程,用有限元法計算各節(jié)點的溫度。下面僅介紹系統(tǒng)的溫度場不隨時間變化的穩(wěn)態(tài)熱傳遞。

傳熱現(xiàn)象由傅里葉定律決定，即物體的瞬態(tài)溫度場T(x, y,z,t)應(yīng)滿足以下方程：

其中ρ為材料密度kg/m3；CT為材料比熱(J/(kg·K)）；k是熱傳導(dǎo)系數(shù)（W/(m·K)）；g為物體內(nèi)部的熱量產(chǎn)生率（W/m3）。

考慮方程服從第三類邊界條件，即：

其中ni為邊界外法線的方向余弦；hi為物體與周圍介質(zhì)的熱交換系數(shù)（W/(m2·K)）；Ta為環(huán)境溫度。對于絕熱系統(tǒng)（系統(tǒng)與外界環(huán)境無熱交換），hi=0。一般情況下，對于不同物質(zhì)材料和溫度，熱傳導(dǎo)系數(shù)k不同。不失一般性，假設(shè)對于每一種材料，k是一個常數(shù)。因此，(1)式就可以寫成：

滿足下面的邊界條件：

初始溫度T=T0，q為在邊界S2上的給定熱流量，h為物體與周圍介質(zhì)的熱交換系數(shù)。運用有限元熱分析方法求解(3)式。通過在對應(yīng)節(jié)點添加相關(guān)數(shù)值，可以得到下面方程來表示整個區(qū)域。

其中f是整個區(qū)域的溫度載荷，C是比熱矩陣，G是剛度矩陣?？紤]穩(wěn)態(tài)問題，即溫度不隨時間變化。

因此，(7)式可以改寫為：

若求出剛度矩陣G的逆矩陣G-1，就可以得到整個模型的熱分布T。

2 Gm sh和H-matrix概述

2.1 Gm sh

Gmsh[7]軟件是一個帶有內(nèi)置在計算機輔助設(shè)計和后期處理器的自動三維有限元網(wǎng)格生成器。Gmsh有4個模塊：幾何模塊，網(wǎng)格化模塊，Solver模塊，后處理模塊。

1)幾何模塊

Gmsh軟件通過邊界表示法來描述幾何實體。通過定義點、有向曲線、有向曲面、體，自下而上的創(chuàng)建模型?；谶@些基礎(chǔ)的幾何實體可以定義幾何實體組。Gmsh軟件的程序描述語言允許所有的幾何實體通過完整的參數(shù)定義來表示。

2)網(wǎng)格化模塊

有限單元網(wǎng)格化將給定的模型細分。網(wǎng)格的生成與幾何模型的創(chuàng)建一樣是自下而上的。首先細分曲線，網(wǎng)格化的曲線被用于網(wǎng)格化曲面，網(wǎng)格化的曲面被用于網(wǎng)格化體。在這個過程中，網(wǎng)格化的實體僅僅被網(wǎng)格化的邊界所限制，這樣可以保證網(wǎng)格化的一致性。

3)Solver模塊

外部解決器可以通過Unix或者傳輸控制協(xié)議及網(wǎng)際協(xié)議與Gmsh相連接，允許修改Solver參數(shù)，啟動外部計算，處理直接從Gmsh后處理模塊得到的結(jié)果。

4)后處理模塊

Gmsh可以加載和操縱多個后處理標(biāo)量、矢量或張量映射。標(biāo)量場由等值線（面）或者顏色圖表示，而矢量場由三維箭頭或者位移圖表示。后處理功能包括部分計算,補償，正面圖，邊界，組件提取，彩色圖，范圍修改，動畫，矢量圖形，輸出等。所有后處理選項都可以通過交互界面或者輸入腳本文件來存取。腳本文件允許所有后處理操作自動化，用戶定義的操作也可以進行后處理。

2.2 H-matrix

H-matrix[8](Hierarchical matrix)，其方法主要目的是使矩陣運算滿足線性復(fù)雜度。將一個高維的復(fù)雜矩陣分層，得到的許多矩陣子集又分解為行列數(shù)更低的矩陣，通過模型的幾何網(wǎng)格化信息，建立大塊的簇樹。H-matrix可以沒有任何誤差地表示剛度矩陣。一個H-matrix H∈CI×I通常與一個可采納條件相關(guān)聯(lián)，其中I={i1,i2,...,in}包含整個有限元模型的形狀函數(shù)φi和節(jié)點信息。子矩陣Ht×s可以通過因式分解表示為Ht×s≈ABT，其中A是一個t×k矩陣，B是一個s×k矩陣。為了構(gòu)建合理的H矩陣結(jié)構(gòu)，基于有限元模型幾何信息構(gòu)建的簇樹TI×I是必須的。應(yīng)用大塊的簇樹構(gòu)造算法，整體模型的剛度矩陣G就可以沒有任何誤差地用H-matrix表示。在確定一定精度下，可以用H-matrix表示剛度矩陣的逆矩陣。

3 建模與仿真

本文的建模與仿真是在Intel i5處理器和4 GB內(nèi)存的仿真環(huán)境下進行的。

3.1 Gm sh創(chuàng)建含TSV的三維集成電路模型

通過Gmsh軟件創(chuàng)建了含有一個TSV的3D-IC模型（如圖1所示）。該模型包括：熱沉，集成電路塊IC1和IC2，TSV，焊球和襯底。其中TSV實現(xiàn)垂直互聯(lián)。

圖1 未網(wǎng)格化的3D-IC模型Fig.1 3D-IC model without grids

3.2 M atlab軟件提取有限元參數(shù),得到剛度矩陣

由上述基于有限元理論的熱模型，利用Matlab軟件提取Gmsh建模的有限元參數(shù)，得到剛度矩陣。圖2是6233行(列)數(shù)目的剛度矩陣截圖。

圖2 6233行（列）數(shù)目剛度矩陣Fig.2 Stiffness matrix with line(column)6233

3.3 可視化的三維集成電路模型熱分布

設(shè)初始溫度為T0=300 K，模型的熱量分布從300～339 K，其中熱源在集成電路塊IC2。圖3～圖5分別表示行（列）數(shù)目為4504、5440、6233的剛度矩陣對應(yīng)模型的熱分布圖。

圖3 4504行(列)數(shù)目剛度矩陣對應(yīng)模型熱分布Fig.3 The thermal distribution of stiffness matrix with line(column) 4504 corresponding to the model

圖4 5440行(列)數(shù)目剛度矩陣對應(yīng)模型熱分布Fig.4 The thermal distribution of stiffness matrix with line(column) 5440 corresponding to the model

圖5 6233行(列)數(shù)目剛度矩陣對應(yīng)模型熱分布Fig.5 The thermal distribution of stiffness matrix with line(column) 6233 corresponding to the model

3.4 不同模型剛度矩陣求逆的計算復(fù)雜度

對行(列)數(shù)較大的矩陣求逆是很耗時的，因此模型的計算復(fù)雜度主要體現(xiàn)在剛度矩陣求逆所消耗的存儲空間和運算時間。表1列出了不同行(列)數(shù)目剛度矩陣求逆耗費的存儲空間和運算時間。

表1 不同模型剛度矩陣求逆耗費的存儲空間和運算時間Tab.1 Different model stiffness matrix inversion cost of storage space and computing time

進一步分析表1中的數(shù)據(jù)，分別用圖6和圖7表示對應(yīng)的存儲空間復(fù)雜度和時間復(fù)雜度。

圖6 求解剛度矩陣的逆耗費的存儲空間Fig.6 To solve the stiffness matrix of the inverse cost of storage space

4 結(jié)果分析

對上述仿真結(jié)果進行分析，得到以下結(jié)論：

1)圖3至圖5較好地體現(xiàn)了模型的熱分布。

2)圖6表明：對剛度矩陣求逆所需的存儲空間隨剛度矩陣行(列)數(shù)目的增加呈線性變化關(guān)系。

3)圖7表明：對剛度矩陣求逆所需的運算時間隨剛度矩陣行(列)數(shù)目的增加呈線性變化關(guān)系。

圖7 求解剛度矩陣的逆耗費的CPU運算時間Fig.7 To solve the stiffness matrix of the inverse consuming CPU operation time

5 結(jié)論

通過Gmsh軟件創(chuàng)建三維集成電路模型并生成網(wǎng)格化文件，利用Matlab軟件提取有限元參數(shù),得到剛度矩陣，利用層次矩陣(H-matrix)表示剛度矩陣。獲得了不同模型的熱分布圖，得到了不同模型剛度矩陣求逆所消耗的存儲空間和運算時間。進一步分析表明，隨剛度矩陣行列數(shù)目的增加，求對應(yīng)逆矩陣所需的運算時間和存儲空間呈線性變化關(guān)系。仿真結(jié)果對三維集成電路的設(shè)計有指導(dǎo)作用。

[1]王乃龍，劉森，周潤德.用有限差分法實現(xiàn)集成電路的電熱耦合模擬[J].微電子學(xué)，2004,34(3):295-297.WANG Nai-long,LIU Miao,ZHOU Run-de.Electro-Thermal simulation of integrated circuits using finite difference method[J].Microelectronics，2004,34(3):295-297.

[2]王鳳娟，朱樟明，王寧，等.考慮硅通孔的三維集成電路最高層溫度模型[J].計算物理，2012(4)：109-113.WANG Feng-juan，ZHU Zhang-ming，WANG Ning，et al.A thermal model for top layer of three-dimensional integrated circuits with through silicon via[J].Chinese Journal of Computational Physics，2012(4)：109-113.

[3]杜秀云，唐禎安.三維集成電路工作熱載荷工況的有限元仿真[J].系統(tǒng)仿真學(xué)報，2012(2)：49-52.DU Xiu-yun,TANG Zhen-an.Finite element simulation for 3D-IC of working process under thermal loads[J].Journal of System Simulation，2012(2)：49-52.

[4]李文石.二維和三維集成電路的熱阻計算 [J].微電子學(xué)，2005(5):40-43.LI Wen-shi.Computation of thermal resistance for 2D-and 3D-IC’s[J].Microelectronics，2005(5):40-43.

[5]胡輝勇，張鶴鳴，賈新章，等.Si-SiGe材料三維CMOS集成電路技術(shù)研究[J].半導(dǎo)體學(xué)報，2007(5):49-53.HU Hui-yong,ZHANG He-ming,JIA Xin-zhang,et al.Study on Si-SiGe three-dimensional CMOS integrated circuits[J].Chinese Journal of Semiconductors，2007(5):49-53.

[6]曾攀.有限元分析基礎(chǔ)教程[M].北京：清華大學(xué)出版社，2008.

[7]Geuzaine C,Jean-Francois Remacle.Gmsh:a three-dimensional finite element mesh generator[J].International Journal for Numerical Methods in Engineering,2009,79(11)：1309-1331.

[8]Mario Bebendorf.Hierarchical Matrices[M].Germany:Springer, 2008.

3D-IC thermal analysis based on the finite element model

DU Wen-xiong,TANG Pu-ying,WANG Shu-bing,ZHANG Si-man
（School of Optoelectronic Information,University of Electronic Science and Technology of China,Chengdu 610054,China）

Three-dimensional integrated circuit(3D-IC)was modeled and simulated based on the finite element model.The thermal distribution and the computation complexity for different model were studied.The 3D-IC models were created and the mesh files were generated by Gmsh.The finite element parameters were extracted and the stiffness matrices were obtained using Matlab.The stiffness matrix was represented by the hierarchical matrix(H-matrix).The storage space and computing time in?solving inverse matrix of stiffness matrix were generated for different model.The results show that the computation time and the storage space for solving different complexity of 3D-IC model present linear relationship as the increase of the number of the stiffness matrix ranks.

Gmsh;H-matrix;3D-IC;finite element thermal analysis;stiffness matrix

TN407

A

1674－6236（2015）10-0079-04

2014-09-24 稿件編號：201409210

杜文雄(1993—)，男，四川廣安人。研究方向：電子技術(shù)。