陳品忠， 潘中良

（華南師范大學(xué)物理與電信工程學(xué)院， 廣州 510006）

1 引言

隨著集成電路特征尺寸的減小，作為集成電路重要組成的銅互連的延遲、噪聲和功耗卻不斷增加。 解決互連延遲（特別是全局互連）對(duì)性能的影響最終只能通過降低互連長度的途徑來實(shí)現(xiàn)。 三維集成技術(shù)在垂直方向進(jìn)行堆疊，相鄰層之間通過垂直互連實(shí)現(xiàn)電學(xué)連接，這為集成電路的互連提供了一個(gè)可行的技術(shù)方案[1]。 三維集成具有更高的器件集成度，芯片面積更小。 然而，三維集成電路的散熱仍舊需要通過芯片的上下表面來進(jìn)行，因?yàn)殡S著芯片面積的減小，散熱能力有所下降，這使得熱管理問題成為限制三維集成電路發(fā)展的主要因素之一[2]。

文獻(xiàn)[3]使用頂層均勻分布熱源，獲得了含有硅通孔等復(fù)雜模型的溫度場。 文獻(xiàn)[4]使用緊湊熱模型分析了3D IC 的橫向傳熱。 文獻(xiàn)[5]分析了在頂層含有非均勻分布熱源情況下的熱模型，并且計(jì)算出了等效的導(dǎo)熱系數(shù)。 文獻(xiàn)[6]討論了在模型材料為各項(xiàng)異性情況下的2 層3D IC 的穩(wěn)態(tài)熱傳導(dǎo)問題的解。 文獻(xiàn)[7]使用迭代模型計(jì)算了3 層3D IC 溫度場的穩(wěn)態(tài)和瞬態(tài)分布，雖然文中指出迭代次數(shù)少， 但是相比非迭代模型，計(jì)算量仍舊偏大。文獻(xiàn)[8]解決了迭代問題，得到多層3D IC 模型的穩(wěn)態(tài)解析解， 但是沒有計(jì)算出瞬態(tài)解析解。文獻(xiàn)[9]只用1 條解析式同時(shí)表示3 層模型穩(wěn)態(tài)和瞬態(tài)的非迭代解，瞬態(tài)解是直接的，但是穩(wěn)態(tài)解是間接的， 因?yàn)榉€(wěn)態(tài)解需要在瞬態(tài)解時(shí)間t 趨于無窮大的情況下得到。

本文分別給出了適用于N 層3D IC 的熱源非均勻分布的穩(wěn)態(tài)解和瞬態(tài)解，其中N≥2，當(dāng)確定源的分布和數(shù)值之后，解析解就可以直接得出；瞬態(tài)時(shí)熱源值的大小可以隨時(shí)間變化。

2 三維集成電路熱模型的建模

三維集成電路熱模型如圖1 所示，圖1（a）表示N層堆疊的熱模型，圖1（b）表示第i 層的矩形熱源分布。模型的長和寬分別為a、b，總厚度為zN。 在三維集成電路中，熱源q 是產(chǎn)生功耗的主要部分，放置在每一層的頂部，矩形形狀契合一般芯片形狀，而數(shù)量可以是任意的。 熱源q 厚度不考慮，因?yàn)樵谌S集成電路中，該部分厚度相對(duì)來說是非常薄的。 散熱片/散熱器在模型的最底部，最底部溫度假設(shè)為25 ℃。在實(shí)際應(yīng)用中，模型側(cè)面、頂面可以假設(shè)為熱絕緣的，因?yàn)樗鼈兊纳嵝Ч鄬?duì)底面會(huì)差很多。

圖1 三維集成電路熱模型

令U1代表穩(wěn)態(tài)溫度場分布函數(shù)，那么穩(wěn)態(tài)熱分析中該模型的熱擴(kuò)散方程可以為：

q1i表示穩(wěn)態(tài)問題中第i 層熱源；Δ 是拉普拉斯算符；系數(shù)k 代表熱導(dǎo)率，盡管多數(shù)材料的熱導(dǎo)率k 都是溫度的函數(shù)，由于芯片工作溫度范圍較小，為了簡單，在這里將每層k 視為同一常數(shù)， 并且是均勻各項(xiàng)同性的。 令U2代表瞬態(tài)溫度場分布函數(shù)，瞬態(tài)熱分析中該模型的熱擴(kuò)散方程可以為：

ρ 代表材料密度，C 為比定壓熱容， 將每層ρ 和C也分別統(tǒng)一為相同的常數(shù)。 q2i表示瞬態(tài)問題中第i 層熱源，t 為時(shí)間。 在穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)問題中，都假設(shè)側(cè)面(x=0,a 以及y=0,b)和最頂面(z=zN)為熱絕緣的，那么函數(shù)U1、U2對(duì)應(yīng)的第二類邊界條件全是齊次的，底面對(duì)應(yīng)的第一類邊界條件非齊次，因此：

層間溫度值連續(xù)，熱流密度值相差qi，由于瞬態(tài)熱擴(kuò)散方程還有U2對(duì)t 求偏導(dǎo)部分，需要添加初始條件才能定解：U2(x,y,z,0)=25。

3 解析解

察邊界條件可以發(fā)現(xiàn)式(3)是非齊次的，考慮將邊界條件齊次化。 對(duì)于穩(wěn)態(tài)問題，構(gòu)造輔助函數(shù)W1(x,y,z)=25，令U1(x,y,z)=T1(x,y,z)+W1(x,y,z)，根據(jù)疊加原理可以得到：

函數(shù)T1(x,y,z)和U1(x,y,z)有類似的齊次邊界條件，即在x=0,a 或y=0,b 或z=zN時(shí)第二類邊界條件是齊次的。在輔助函數(shù)的影響下，函數(shù)T1(x,y,z)的第一類邊界條件也是齊次的，即T1(x,y,0)=0。 因此，偏微分方程式(4)的定解條件全是齊次的。

對(duì)于瞬態(tài)問題，構(gòu)造輔助函數(shù)W2(x,y,z,t)=25，令U2(x,y,z,t)=T2(x,y,z,t)+W2(x,y,z,t)，根據(jù)疊加原理可以得到：

和穩(wěn)態(tài)問題類似，偏微分方程式(5)的定解條件也全是齊次的，不僅包括邊界條件T2(x,y,0,t)=0，還包括初始條件T2(x,y,z,0)=0。

因?yàn)檩o助函數(shù)W1、W2是已知的， 所以對(duì)函數(shù)U1、U2的求解也就轉(zhuǎn)變?yōu)閷?duì)函數(shù)T1、T2的求解。 以下主要討論關(guān)于函數(shù)T1、T2的求解過程。

3.1 穩(wěn)態(tài)解析解

文獻(xiàn)[9]對(duì)復(fù)雜源的溫度場提供了格林函數(shù)的思想，可以將格林函數(shù)方法運(yùn)用到這里的求解中。 格林函數(shù)一般討論點(diǎn)源對(duì)場的影響，但上述模型熱源分布相對(duì)簡單，為了方便計(jì)算，令函數(shù)G1i(x,y,z；zi)表示在穩(wěn)態(tài)問題中z=zi時(shí)該層頂面熱源對(duì)模型溫度場產(chǎn)生的影響，根據(jù)式(4)可以得到一個(gè)含有單位沖激函數(shù)δ 的等式：ΔG1i=-qi(x,y)/kδ(z-zi)。 那么函數(shù)T1和函數(shù)G1i的關(guān)系則為：

下面來求解函數(shù)G1i(x,y,z；zi)，根據(jù)分離變量法、本征函數(shù)展開法和本征函數(shù)正交性得到：

其中n、m 為任意非負(fù)整數(shù)，α= [(nπ)/a]2，β=[(mπ)/b]2，積分微元dσ 的積分區(qū)域是第i 層頂面區(qū)域(如圖1（b） 所示)，δn0、δm0是kronecker 函數(shù)，

根據(jù)式(6)、(7)，可獲得函數(shù)T1(x,y,z)的解。

3.2 瞬態(tài)解析解

根據(jù)格林函數(shù)思想，令函數(shù)G2i(x,y,z,t；zi)表示在瞬態(tài)問題中z=zi時(shí)該層頂面熱源對(duì)模型溫度場產(chǎn)生的影響，根據(jù)式(5)可以得到以下等式：ΔG2i+q2i(x,y,t)/kδ(z-zi)=1/μ?G2i/?t。 那么函數(shù)T2和函數(shù)G2i的關(guān)系則為：

其中μ=k/(ρC)。 假設(shè)熱源和時(shí)間的關(guān)系為q2i(x,y,t)=從1 到Ji，Ji表示 第i 層 熱源數(shù) 量，p2ij(x,y)表示第i 層第j 個(gè)熱源，功率大小為1 W/mm2，qij(t)表示第i 層第j 個(gè)熱源的值隨時(shí)間變化函數(shù)，當(dāng)qij(t)為常量時(shí)表示熱源功率不隨時(shí)間變化。 使用分離變量法、本征函數(shù)展開法和本征函數(shù)正交性可以得到：

其中n、m、p 為任意非負(fù)整數(shù)，α= [(nπ)/a]2，β=[(mπ)/b]2，γ=[(p+1/2)π/zN]2，λ=α+β+γ，h(t)是與時(shí)間t 有關(guān)的單位階躍函數(shù)積分微元dv 的積分區(qū)域是整個(gè)模型，如圖1(a)所示，δn0、δm0是kronecker函數(shù)。

根據(jù)式(8)、(9)就能夠求出T2(x,y,z,t)的解。

4 實(shí)驗(yàn)結(jié)果與分析

為了驗(yàn)證解的準(zhǔn)確性，在這部分將給出穩(wěn)態(tài)解和瞬態(tài)解在MATLAB 中計(jì)算得到的結(jié)果， 并與COMSOL 得到的結(jié)果進(jìn)行比較。首先說明實(shí)驗(yàn)統(tǒng)一使用的參數(shù)： 模型長a 和寬b 都設(shè)為10 mm；z1、z2、z3、z4、z5、z6分別為3 mm、4 mm、5 mm、7 mm、8 mm、9 mm，3層堆疊時(shí)模型最高為7 mm (z4)，5 層堆疊時(shí)模型最高為9 mm (z6)；k 為1 W/(mm·K)，ρ 為1 kg/mm3，C 為1 J/(kg·K)；使用大中小3 種尺寸的熱源，大熱源1 mm×8 mm，2 W/mm2，中熱源1 mm×2 mm，10 W/mm2，小熱源1 mm×1 mm，20 W/mm2。

4.1 穩(wěn)態(tài)熱傳導(dǎo)問題

穩(wěn)態(tài)問題第1～5 層熱源位置分布如圖2 所示。 第1 層小熱源幾何中心分別在(2 mm,6 mm)、(7 mm，5.5 mm)，中熱源幾何中心分別在(3 mm,2 mm)、(5 mm,4 mm)，大熱源幾何中心在(5 mm,8.5 mm)；第2 層小熱源幾何中心在(7 mm,7 mm)，中熱源幾何中心在(3 mm,5.5 mm)，大熱源幾何中心在(5 mm,2.5 mm)；第3 層小熱源幾何中心分別在(5 mm,8 mm)、(3 mm,5.5 mm)，中熱源幾何中心在(6 mm,4 mm)；第4 層大熱源幾何中心在(7.5 mm,4.5 mm)；第5 層小熱源幾何中心分別在(5 mm,5 mm)、(8.5 mm,1.5 mm)。 3 層堆疊時(shí)只提供第1、2、3 層熱源，頂層(z=7 mm)不布置熱源；5 層堆疊時(shí)包含1～5 層全部熱源。

圖2 穩(wěn)態(tài)問題第1～5 層熱源位置分布

由于解析解是無窮級(jí)數(shù)，MATLAB 在工作時(shí)需要確定n、m 的值。 為了使結(jié)果誤差在一定范圍內(nèi)而n、m的值又不會(huì)太大，需要進(jìn)行調(diào)試。在使用上述3 層堆疊并且n、m 值相同的情況下，級(jí)數(shù)上界采用5、11、15、21共4 組數(shù)據(jù)進(jìn)行比較： 在穩(wěn)態(tài)解情況下，COMSOL 得到的峰值溫度約為39.030 ℃，MATLAB 得到的峰值溫度分別為36.186 ℃、38.786 ℃、39.395 ℃、39.640 ℃，所以誤差分別為7.29%、0.63%、0.94%、1.56%。 提取模型中峰值溫度所在的x-y 平面如圖3 所示，圖3(a)是n=m=15 時(shí)的MATLAB 溫度分布圖， 圖3 (b)是COMSOL 溫度分布圖，在圖中過峰值溫度所在點(diǎn)做平行于x 軸的虛線。 圖4 除“ COMSOL”以外的4 條溫度曲線對(duì)應(yīng)圖3(a)虛線上使用不同求和上界求得的溫度值，圖4“ COMSOL”溫度曲線對(duì)應(yīng)圖3(b)虛線上的溫度值。 “ n=m=15”和“ n=m=11”峰值溫度誤差都低于1%， 而“ n=m=15” 與“ n=m=21” 這2 條溫度曲線和“ COMSOL”溫度曲線吻合度都非常高，綜合以上2 點(diǎn)因素，穩(wěn)態(tài)解選用求和上界n=m=15。

圖3 穩(wěn)態(tài)問題3 層3DIC 峰值溫度所在的x-y 平面溫度分布

圖4 y=5.5 mm,z=5 mm, 沿x 方向3 層3D IC 不同求和上界的溫度曲線

穩(wěn)態(tài)問題沿x 方向的溫度曲線如圖5 所示。 在3層3D IC 模型中令z=4 mm，y=5 mm，沿x 方向提取溫度值，得到“ 3 層z=4 mm”溫度曲線；在5 層3D IC 模型中令z=7 mm，y=5 mm，沿x 方向提取溫度值，得到“ 5 層z=7 mm”溫度曲線；在5 層3D IC 模型中找到峰值溫度所在x-y 平面，即z=8 mm 時(shí)的x-y 平面，根據(jù)峰值溫度所在的幾何點(diǎn)固定y 值(通過仿真確定y=1.5 mm)，沿x 方向提取溫度值，得到“ 5 層z=8 mm”溫度曲線，“ 5 層峰值溫度” 指的是5 層堆疊情況下模型最大溫度值。 可以發(fā)現(xiàn)，“ 3 層z=4 mm”、“ 5 層z=7 mm” MATLAB 溫度曲線和COMSOL 溫度曲線基本是重合的， 通過MATLAB 計(jì)算得到這2 組曲線最大誤差分別為0.30%、0.08%?！?5 層峰值溫度”橢圓標(biāo)記處溫差比較大， 在MATLAB 計(jì)算得到該處誤差為1.07%。

圖5 穩(wěn)態(tài)問題沿x 方向的溫度曲線

穩(wěn)態(tài)問題沿z 方向溫度曲線如圖6 所示。 在3 層3D IC 模型中令x=5 mm、y=5 mm， 沿z 方向提取溫度值，得到“ 3 層x=5 mm”溫度曲線；在5 層3D IC 模型中令x=5 mm、y=5 mm，沿z 方向提取溫度值，得到“ 5層x=5 mm”溫度曲線。3 層堆疊模型高度只有7 mm，5層堆疊模型高度為9 mm，所以“ 3 層x=5 mm”溫度曲線在7≤z≤9 區(qū)間沒有溫度值。 觀察圖6 可以發(fā)現(xiàn)，MATLAB 曲線和COMSOL 曲線吻合度非常高，在MATLAB 中計(jì)算“ 3 層x=5 mm”和“ 5 層x=5 mm”曲線最大誤差分別為0.33%、1.03%。 對(duì)比“ 3 層x=5 mm”曲線和“ 5 層x=5 mm”曲線可以發(fā)現(xiàn)，在3 層3D IC 模型基礎(chǔ)上再增加2 層熱源， 這會(huì)導(dǎo)致第1、2、3 層溫度整體性提高。 圖6 中z=8 mm 時(shí)COMSOL 和MATLAB曲線都出現(xiàn)了尖點(diǎn)，即一階不可導(dǎo)點(diǎn)，這是熱源在z 方向上沒有厚度導(dǎo)致的。

圖6 穩(wěn)態(tài)問題沿z 方向溫度曲線

4.2 瞬態(tài)熱傳導(dǎo)問題

瞬態(tài)問題第1～5 層熱源位置分布如圖7 所示。 在瞬態(tài)解中各層模型尺寸和3 種熱源尺寸都和穩(wěn)態(tài)解的一樣，只是熱源個(gè)數(shù)和分布不一樣。第1 層小熱源的幾何中心分別為(4 mm,7.5 mm)、(5 mm,5 mm)、(3 mm，3 mm)，中熱源的幾何中心為(6 mm,3 mm)，大熱源的幾何中心為(8 mm,5 mm)；第2 層中熱源幾何中心分別為(6.5 mm,7.5 mm)、(3 mm,5 mm)；第3 層中熱源幾何中心為(5 mm,6.5 mm)，大熱源幾何中心為(8 mm，5 mm)；第4 層小熱源幾何中心為(5 mm,1.5 mm)，中熱源幾何中心為(3 mm,8 mm)，大熱源幾何中心為(5 mm,5 mm)；第5 層小熱源幾何中心分別為(6.5 mm,8 mm)、(3.5 mm,3.5 mm)，中熱源幾何中心為(3.5 mm,6.5 mm)。5 層3D IC 包括上述所有熱源，3 層3D IC 沒有第4、5層熱源。

圖7 瞬態(tài)問題第1～5 層熱源位置分布

MATLAB 使用“ n=m=5，p=11”、“ n=m=11，p=21”、“ n=m=15，p=31”、“ n=m=21，p=41”4 組求和上界計(jì)算得到在t=1000 s 時(shí)3 層3D IC 模型峰值溫度分別為32.245 ℃、34.357 ℃、35.182 ℃、35.501 ℃，COMSOL仿真峰值溫度為35.423 ℃， 對(duì)應(yīng)誤差分別為8.97%、3.01%、0.68%、0.22%。3 層3D IC 模型峰值溫度所在的x-y 平面(z=3 mm)如圖8 所示。 在這些平面上令y=5 mm，作平行于x 軸的虛線，提取虛線上的溫度值得到MATLAB 中不同求和上界的溫度曲線和COMSOL 中的溫度曲線，結(jié)果如圖9 所示?？梢园l(fā)現(xiàn)，當(dāng)求和上界值越大，MATLAB 溫度曲線和COMSOL溫度曲線的吻合度越高，并且“ n=m=21，p=41”曲線峰值溫度誤差最小，所以瞬態(tài)解選用n=m=21，p=41。

圖8 t=1000 s，3 層3D IC 峰值溫度所在的x-y 平面溫度分布

圖9 t=1000 s,y=5 mm,z=3 mm, 沿x 方向使用不同求和上界的3 層3D IC 溫度曲線

t=100 s 時(shí)沿x 方向的溫度曲線如圖10 所示。“ 3層z=4 mm” 曲線是3 層3D IC 模型z=4 mm、y=5 mm時(shí)沿x 方向的溫度曲線，最大誤差為0.76%；“ 5 層z=7 mm”曲線是5 層3D IC 模型z=7 mm、y=5 mm 時(shí)沿x方向的溫度曲線，最大誤差為2.22%。t=100 s 時(shí)沿z 方向的溫度曲線如圖11 所示?！?3 層x=5 mm” 曲線是3層3D IC 模型x=5 mm、y=5 mm 時(shí)沿z 方向的溫度曲線，最大誤差為1.68%；“ 5 層x=5 mm”曲線是5 層3D IC 模型x=5 mm、y=5 mm 時(shí)沿z 方向的溫度曲線，最大誤差為2.49%。 因?yàn)? 層3D IC 在z 方向上最大只有7 mm，所以“ 3 層x=5 mm”溫度曲線在7＜z≤9 時(shí)沒有溫度值。 由于熱源設(shè)置為沒有厚度，COMSOL 曲線z=3 mm 時(shí)出現(xiàn)尖點(diǎn)，但是從式(9)可知，瞬態(tài)解無法產(chǎn)生和z 有關(guān)的尖點(diǎn)，因?yàn)榭蓪?dǎo)，瞬態(tài)解在MATLAB 中只能通過不斷增加求和上界的方式來無限逼近尖點(diǎn)形狀。 圖11 和圖6 相比，可以發(fā)現(xiàn)穩(wěn)態(tài)解卻能夠在MATLAB 中成功產(chǎn)生尖點(diǎn)，根據(jù)式(7)可知，穩(wěn)態(tài)解是以分段函數(shù)的方式來產(chǎn)生尖點(diǎn)的。

圖10 t=100 s 時(shí)沿x 方向的溫度曲線

以上討論的熱源都是與時(shí)間無關(guān)的，現(xiàn)在使用瞬態(tài)問題中的5 層3D IC 模型， 全部熱源幅值比例設(shè)為a， 并且只將第1 層熱源的值改為與時(shí)間有關(guān)的，將2 W/mm2的熱源改為a [2+sin (bt))] W/mm2，將10 W/mm2的熱源改為a[10+ sin(bt)] W/mm2，將20 W/mm2的熱源改為a[20+sin(bt)]W/mm2。在這樣的設(shè)置前提下得到的實(shí)驗(yàn)結(jié)果如圖12 所示。 MATLAB和COMSOL 在幾何中心點(diǎn) (x=5 mm，y=5 mm, z=4.5 mm)的溫度變化曲線是十分相似的：當(dāng)a=1、b=0.01時(shí)， 計(jì)算結(jié)果和仿真結(jié)果之間最大誤差為2.30%；當(dāng)a=1.5、b=0.015 時(shí)，最大誤差為2.60%；當(dāng)a=2、b=0.02時(shí)，最大誤差為2.50%。

圖11 t=100 s 時(shí)沿z 方向的溫度曲線

圖12 x=5 mm, y=5 mm,z=4.5 mm 幾何點(diǎn)溫度隨時(shí)間變化曲線

5 結(jié)論

本文基于一個(gè)多層3D IC 模型，給定適當(dāng)?shù)倪吔鐥l件，通過分離變量法來解決模型的熱傳導(dǎo)問題，并且分別給出穩(wěn)態(tài)解析解和瞬態(tài)解析解。 該模型可以任意設(shè)置熱源的尺寸、擺放位置，功率大小可以隨時(shí)間變化。用3 層和5 層3D IC 模型驗(yàn)證解析解的正確性，解析解的溫度場用MATLAB 表示出來， 并且與COMSOL 仿真得到的溫度場進(jìn)行比較。 結(jié)果表明，穩(wěn)態(tài)解析解計(jì)算3 層3D IC 模型得到的溫度結(jié)果和COMSOL 仿真結(jié)果之間的誤差小于1%，計(jì)算5 層3D IC 模型誤差在1%左右； 瞬態(tài)解析解計(jì)算3 層3D IC模型得到的溫度結(jié)果和COMSOL 仿真結(jié)果誤差小于2%，計(jì)算5 層3D IC 模型誤差小于3%。這些研究結(jié)果可以應(yīng)用于三維集成電路在預(yù)測(cè)溫度分布方面的設(shè)計(jì)過程中。