李含靈 曹炳陽

(清華大學(xué)工程力學(xué)系,熱科學(xué)與動(dòng)力工程教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

體點(diǎn)導(dǎo)熱問題是電子器件散熱優(yōu)化方面的基礎(chǔ)問題之一,已有研究大多建立在傅里葉導(dǎo)熱定律的基礎(chǔ)上,但隨著電子器件的特征尺度降低到微納米量級(jí),導(dǎo)熱優(yōu)化需要考慮非傅里葉效應(yīng).本文結(jié)合聲子玻爾茲曼方程的數(shù)值解和對聲子平均自由程進(jìn)行插值的固體各向同性材料懲罰方法,發(fā)展了微納米尺度下彈道-擴(kuò)散導(dǎo)熱的拓?fù)鋬?yōu)化方法.在彈道-擴(kuò)散導(dǎo)熱機(jī)制下,體點(diǎn)導(dǎo)熱的拓?fù)鋬?yōu)化得到的材料分布明顯不同于擴(kuò)散導(dǎo)熱下的樹狀分布,且會(huì)隨努森數(shù)的變化而變化,與拓?fù)鋬?yōu)化的插值方式和聲子彈道輸運(yùn)有關(guān).隨著彈道效應(yīng)的增強(qiáng),尺寸效應(yīng)使得材料分布中微小結(jié)構(gòu)的等效熱導(dǎo)率低于粗壯結(jié)構(gòu),因而拓?fù)鋬?yōu)化結(jié)果朝著增多粗壯結(jié)構(gòu)、減少微小結(jié)構(gòu)的方向發(fā)展.彈道效應(yīng)足夠強(qiáng)時(shí),填充材料聚集在低溫邊界附近,主干和枝合并,呈團(tuán)狀分布.

1 引 言

以芯片為代表的電子器件,在當(dāng)代社會(huì)中扮演著不可替代的重要角色.隨著電子器件向小型化、集成化的方向發(fā)展,器件的特征尺寸已經(jīng)降低至微納米量級(jí)[1],器件內(nèi)的功率密度大幅增加.2012年,芯片內(nèi)的平均熱流密度便已達(dá)到250 W/cm2左右[2].高功率密度會(huì)讓電子器件溫度升高,產(chǎn)生局部熱點(diǎn),導(dǎo)致其可靠性和壽命顯著下降[3].因此,電子器件的熱管理已經(jīng)成為相關(guān)領(lǐng)域發(fā)展面臨的關(guān)鍵挑戰(zhàn)之一[4-7].根據(jù)電子器件內(nèi)熱點(diǎn)分布和傳熱結(jié)構(gòu)布局等因素,進(jìn)行散熱仿真和優(yōu)化改進(jìn),對實(shí)際產(chǎn)品的設(shè)計(jì)和制造有著重要意義.

在芯片內(nèi)部填充高導(dǎo)熱材料以加強(qiáng)導(dǎo)熱是提高芯片散熱能力的主要途徑,考慮到芯片內(nèi)空間寶貴,高導(dǎo)熱材料的添加數(shù)量必須受到限制,這樣的芯片散熱問題可以被抽象為體點(diǎn)(volume-to-point,VP)導(dǎo)熱問題[8],如圖1所示.給定體積的區(qū)域代表芯片,區(qū)域內(nèi)的均勻內(nèi)熱源代表焦耳加熱,產(chǎn)生的熱量通過邊界處一小段溫度為 T0的低溫?zé)岢亮鞒?其余邊界絕熱.現(xiàn)提供數(shù)量有限的高導(dǎo)熱填充材料,希望尋找合適的填充材料分布方式,構(gòu)造高導(dǎo)熱通道,讓整個(gè)系統(tǒng)的溫度(平均溫度或最高溫度)降到最低.關(guān)于體點(diǎn)導(dǎo)熱問題,人們采用多種優(yōu)化方法進(jìn)行了研究,包括構(gòu)型理論[8]、仿生優(yōu)化[9]、拓?fù)鋬?yōu)化[10-13]、模擬退火和遺傳算法[14]等.盡管這些方法都能有效降低系統(tǒng)溫度、提高散熱能力,但所得的高導(dǎo)熱填充材料分布方式不盡相同,這與優(yōu)化問題自身帶有的局部最優(yōu)性和多解性有關(guān).當(dāng)優(yōu)化問題改變時(shí),不同方法的性能優(yōu)劣也可能發(fā)生變化.在各類優(yōu)化方法中,拓?fù)鋬?yōu)化方法因其能夠改變結(jié)構(gòu)的拓?fù)錁?gòu)型、具備理論上的最高設(shè)計(jì)自由度而受到人們的青睞.利用拓?fù)鋬?yōu)化,設(shè)計(jì)者不僅可以實(shí)現(xiàn)設(shè)計(jì)對象某方面性能的顯著提升,而且所花費(fèi)的設(shè)計(jì)時(shí)間更少,能得到傳統(tǒng)經(jīng)驗(yàn)之外的具有啟發(fā)性的新結(jié)果[15].體點(diǎn)導(dǎo)熱問題的拓?fù)鋬?yōu)化會(huì)得到填充材料充分伸入內(nèi)部區(qū)域、整體呈樹狀的材料分布,其散熱性能明顯優(yōu)于其他方法所預(yù)測的結(jié)構(gòu)[11].

圖1 體點(diǎn)導(dǎo)熱問題示意圖Fig.1.The schematic diagram of the VP problem.

已有的導(dǎo)熱拓?fù)鋬?yōu)化方法,都建立在經(jīng)典的傅里葉導(dǎo)熱定律的基礎(chǔ)上[13],這在宏觀尺度下是合理的,但對微納米尺度下的導(dǎo)熱過程,傅里葉導(dǎo)熱定律不再成立[16-22],熱仿真及熱優(yōu)化必須考慮非傅里葉效應(yīng).聲子是電子器件半導(dǎo)體中的主要載熱子[23],聲子彈道輸運(yùn)和聲子相干是微納米尺度下發(fā)生非傅里葉導(dǎo)熱的兩個(gè)重要原因[5].聲子彈道輸運(yùn)的強(qiáng)度可由努森數(shù) Kn 來表征,其定義為材料聲子平均自由程l和系統(tǒng)特征尺寸L的比值,即Kn=l/L.宏觀尺度下,系統(tǒng)尺寸遠(yuǎn)大于聲子平均自由程(Kn→0),聲子間散射充分,聲子按擴(kuò)散方式傳遞,導(dǎo)熱過程遵循傅里葉導(dǎo)熱定律.但當(dāng)系統(tǒng)尺寸與聲子平均自由程相當(dāng)(Kn～1)甚至更小(Kn＜1)時(shí),聲子不經(jīng)歷散射而直接抵達(dá)邊界,這樣的過程稱為彈道輸運(yùn),會(huì)導(dǎo)致傅里葉導(dǎo)熱定律失效[24].在微納米結(jié)構(gòu)中,部分聲子發(fā)生彈道輸運(yùn),其余聲子仍按擴(kuò)散方式運(yùn)動(dòng),對應(yīng)的導(dǎo)熱機(jī)制稱為彈道-擴(kuò)散導(dǎo)熱[25,26].聲子彈道輸運(yùn)會(huì)引起熱導(dǎo)率的尺寸效應(yīng),表現(xiàn)為系統(tǒng)的等效熱導(dǎo)率與系統(tǒng)尺寸正相關(guān),尺寸減小會(huì)造成等效熱導(dǎo)率降低[27,28],系統(tǒng)熱點(diǎn)溫度升高[29].此外,在給定溫度邊界處,由于彈道輸運(yùn)的非平衡本質(zhì),邊界處的溫度并不連續(xù),會(huì)出現(xiàn)溫度跳躍現(xiàn)象[30,31].聲子相干則主要出現(xiàn)在特征尺寸和聲子波長相當(dāng)?shù)南到y(tǒng)中[32],例如石墨烯這樣的二維材料和超晶格這樣的特殊結(jié)構(gòu)[33-35].聲子相干通過改變聲子群速度、弛豫時(shí)間和態(tài)密度等因素來影響導(dǎo)熱過程[36].對常用的半導(dǎo)體材料硅,室溫下聲子平均自由程約為300 nm[37],聲子波長小于10 nm[17],而目前大多數(shù)電子器件的特征尺寸在100 nm量級(jí)[1],故本文主要關(guān)注聲子彈道輸運(yùn)引起的非傅里葉導(dǎo)熱.關(guān)于微納米尺度下彈道-擴(kuò)散導(dǎo)熱過程的優(yōu)化,前人曾結(jié)合動(dòng)力學(xué)理論和拓?fù)鋬?yōu)化來處理[38],但其優(yōu)化目標(biāo)只是特定點(diǎn)之間的溫差,無法考慮系統(tǒng)整體的溫度.目前,尚缺少適用于微納米尺度下彈道-擴(kuò)散導(dǎo)熱過程的一般性優(yōu)化方法,關(guān)于彈道效應(yīng)對優(yōu)化結(jié)果的影響也缺乏系統(tǒng)的認(rèn)識(shí).

本文結(jié)合聲子玻爾茲曼輸運(yùn)方程(Boltzmann transport equation,BTE)和固體各向同性材料懲罰(solid isotropic material with penalization,SIMP)方法,發(fā)展了微納尺度下彈道-擴(kuò)散導(dǎo)熱的拓?fù)鋬?yōu)化方法.用聲子BTE代替傅里葉導(dǎo)熱定律作為導(dǎo)熱本構(gòu)方程,以反映聲子彈道輸運(yùn); 用SIMP方法對聲子平均自由程的倒數(shù)進(jìn)行插值以引入拓?fù)鋬?yōu)化,并通過對相對密度變量全局梯度施加顯示約束的方法來確保解的適定性和網(wǎng)格無關(guān)性.將此方法用于體點(diǎn)導(dǎo)熱優(yōu)化,發(fā)現(xiàn)彈道-擴(kuò)散導(dǎo)熱機(jī)制下,拓?fù)鋬?yōu)化預(yù)測的填充材料分布明顯不同于傅里葉定律下經(jīng)典的樹狀分布,且會(huì)隨 Kn 的變化而變化,這與SIMP方法的插值方式和聲子彈道輸運(yùn)有關(guān).

2 優(yōu)化方法

2.1 聲子玻爾茲曼輸運(yùn)方程

聲子的運(yùn)動(dòng)規(guī)律符合聲子BTE,穩(wěn)態(tài)、弛豫時(shí)間近似下的聲子BTE寫作

其中f是聲子分布函數(shù); vg是群速度矢量; f0是平衡態(tài)聲子分布函數(shù),即玻色-愛因斯坦分布; τ 是弛豫時(shí)間.聯(lián)合求解(1)式和能量守恒方程,就可以獲得聲子的微觀運(yùn)動(dòng)規(guī)律和分布函數(shù),進(jìn)一步可計(jì)算得到溫度、熱導(dǎo)率等宏觀熱性質(zhì).

本文選擇離散坐標(biāo)法(discrete ordinate method,DOM)[39-41]來求解聲子BTE,此方法最早用于輻射輸運(yùn)方程的求解,相關(guān)計(jì)算方法的發(fā)展較為成熟.根據(jù)聲子和光子的類比,可定義聲子輻射強(qiáng)度其中ω表示聲子頻率,D(ω)表示態(tài)密度,“”表示對聲子極化的求和.利用此定義,可以將聲子BTE轉(zhuǎn)化為聲子的輻射輸運(yùn)方程(equation of phonon radiative transfer,EPRT)[42]

(3)式中包含N+1個(gè)等式,與未知的聲子輻射強(qiáng)度數(shù)量相對應(yīng),方程組封閉.

得到聲子輻射強(qiáng)度后便可計(jì)算溫度.聲子發(fā)生彈道輸運(yùn)時(shí),當(dāng)?shù)責(zé)崃W(xué)條件遠(yuǎn)離平衡態(tài),傳統(tǒng)的在熱平衡態(tài)下定義的溫度失去意義[43],而如何在這樣的非平衡態(tài)下定義溫度一直存在爭議[44].本文采用等效平衡溫度的概念[45,46],認(rèn)為空間內(nèi)某個(gè)離散單元在發(fā)射能量時(shí),產(chǎn)生的聲子服從玻色-愛因斯坦分布,根據(jù)發(fā)射能量反解出的分布函數(shù)中的溫度,便是等效平衡溫度.使用DOM離散輻射強(qiáng)度后,等效平衡溫度的計(jì)算式為其中 σphonon是聲子斯特藩-玻爾茲曼參數(shù)[47].使用這種等效平衡溫度的好處是,在接近擴(kuò)散導(dǎo)熱時(shí),解聲子BTE所得結(jié)果能夠與傅里葉導(dǎo)熱定律的結(jié)果相吻合[39].

2.2 拓?fù)鋬?yōu)化

拓?fù)鋬?yōu)化本質(zhì)上屬于大規(guī)模0-1整數(shù)規(guī)劃問題,區(qū)域內(nèi)的每個(gè)點(diǎn)要么是空洞,要么是實(shí)體.這樣的離散優(yōu)化問題在數(shù)學(xué)上是病態(tài)的,解的存在性、收斂性以及求解算法的實(shí)現(xiàn)都存在著巨大的困難[48].實(shí)現(xiàn)拓?fù)鋬?yōu)化的關(guān)鍵是對優(yōu)化問題進(jìn)行適當(dāng)?shù)恼齽t化處理,將設(shè)計(jì)變量由離散的變?yōu)檫B續(xù)的,從而能夠直接使用針對連續(xù)設(shè)計(jì)變量優(yōu)化問題的求解方法.在這方面,SIMP方法[49]是經(jīng)典而常用的技巧,其核心思想是將材料性質(zhì)設(shè)定為關(guān)于相對密度的冪函數(shù),從而將設(shè)計(jì)變量轉(zhuǎn)變?yōu)檫B續(xù)的相對密度.考慮到(2)式作為描述非傅里葉導(dǎo)熱的本構(gòu)方程,僅有聲子消光系數(shù) κ 是與材料類型有關(guān)的性質(zhì)[38],于是本文提出對 κ 進(jìn)行插值的SIMP方法來實(shí)現(xiàn)拓?fù)鋬?yōu)化.空間中任意位置的消光系數(shù)寫作

其中 ρ 是值在[0,1]區(qū)間上連續(xù)分布的相對密度變量; κs和 κf分別是基底材料和填充材料的消光系數(shù),亦即各自聲子平均自由程的倒數(shù); p是值大于1的懲罰因子.懲罰因子的作用是在優(yōu)化過程中對介于0和1之間的中間密度值進(jìn)行懲罰,使得 ρ 的值逐漸趨向于0和1兩個(gè)端點(diǎn)值,從而在優(yōu)化結(jié)果中得到清晰的空洞(ρ=0)和實(shí)體(ρ=1)分布.關(guān)于SIMP方法插值方式的選擇,將在下一節(jié)結(jié)合優(yōu)化結(jié)果做進(jìn)一步討論.

即使采用SIMP方法,拓?fù)鋬?yōu)化仍會(huì)面臨解的不存在性、不收斂性和不唯一性等問題,導(dǎo)致優(yōu)化結(jié)果會(huì)出現(xiàn)棋盤格、網(wǎng)格依賴性等問題[50].為了得到可靠的優(yōu)化結(jié)果,需要引入更多的正則化處理以抑制數(shù)值不穩(wěn)定性.本文采用對相對密度變量 ρ 的全局梯度施加顯式約束的方法[15],即在目標(biāo)函數(shù)中增加關(guān)于 ρ 的全局梯度的權(quán)重懲罰項(xiàng),以約束ρ的變化程度.至此,可寫出通過求解EPRT來獲得溫度分布的二維體點(diǎn)導(dǎo)熱拓?fù)鋬?yōu)化的數(shù)學(xué)模型為：

其中 Tave是系統(tǒng)平均溫度,Tave,ref是人為選取的參考溫度,α 是相對密度變量全局梯度的權(quán)重系數(shù),是無量綱的相對密度變量梯度,A=a2是二維區(qū)域的面積,? 是預(yù)先給定的高導(dǎo)熱填充材料占整個(gè)區(qū)域的比例.對平均溫度和設(shè)計(jì)變量梯度進(jìn)行無量綱化的目的是使二者通過線性求和構(gòu)成實(shí)際的目標(biāo)函數(shù)g.現(xiàn)在的模型以降低系統(tǒng)平均溫度為優(yōu)化目標(biāo),如果目標(biāo)是減小系統(tǒng)的最高溫度,則可以將優(yōu)化對象變更為文獻(xiàn)[51]中提到的幾何平均溫度,其他設(shè)置不變.需要說明的是,因?yàn)楸疚牡年P(guān)注點(diǎn)是如何在導(dǎo)熱優(yōu)化中考慮聲子彈道輸運(yùn)并分析彈道輸運(yùn)對優(yōu)化結(jié)果的影響,所以優(yōu)化模型中并未考慮不同材料之間的界面熱阻,這樣可以更直觀地反映彈道效應(yīng)的影響.界面熱阻對現(xiàn)有優(yōu)化結(jié)果的影響將在第3節(jié)中進(jìn)行討論.

求解(5)式的拓?fù)鋬?yōu)化的流程圖如圖2所示,包括以下步驟：1)初始化,對區(qū)域進(jìn)行離散并任意給定一組 ρ 的初始值; 2)系統(tǒng)重分析,根據(jù)SIMP方法插值得到當(dāng)前設(shè)計(jì)點(diǎn)對應(yīng)的消光系數(shù),再求解EPRT和穩(wěn)態(tài)能量守恒方程獲得聲子輻射強(qiáng)度,進(jìn)而計(jì)算平均溫度和目標(biāo)函數(shù); 3)靈敏度分析,計(jì)算目標(biāo)函數(shù)和約束函數(shù)對設(shè)計(jì)變量的導(dǎo)數(shù); 4)優(yōu)化求解,根據(jù)導(dǎo)數(shù)信息,確定滿足約束且減小目標(biāo)函數(shù)的設(shè)計(jì)變量演化方向(可行下降方向),并計(jì)算最佳步長,更新當(dāng)前設(shè)計(jì)變量值; 5)收斂判斷,滿足收斂條件時(shí)結(jié)束迭代,輸出結(jié)果; 否則重復(fù)步驟2)-4).

圖2 體點(diǎn)導(dǎo)熱拓?fù)鋬?yōu)化的流程圖Fig.2.The flow chart of topology optimization for the VP problem.

3 結(jié)果與分析

在體點(diǎn)導(dǎo)熱問題中,設(shè)置低溫邊界與區(qū)域邊長的比值為 c/a=0.1 ,低溫邊界的溫度為 T0=300K ,填充材料體積占比為 ?=0.1 ,設(shè)計(jì)變量初始值取為 ρ=? ,即材料均勻分布.定義無量綱空間坐標(biāo)為 X=x/a ,Y=y/a ,參考文獻(xiàn)[12]定義無量綱溫度其中 ks是基底材料熱導(dǎo)率,這樣無量綱化的好處是消除了內(nèi)熱源功率、熱源面積、熱導(dǎo)率和邊界溫度具體值對優(yōu)化結(jié)果的影響,使填充材料和基底材料間導(dǎo)熱性能差異成為優(yōu)化結(jié)果的決定因素.在數(shù)值計(jì)算中,用數(shù)量為n×n的方形網(wǎng)格離散計(jì)算域,網(wǎng)格尺寸為 h=a/n.DOM的離散方向數(shù)設(shè)定為 N=24 ,此值既能避免求解溫度場時(shí)出現(xiàn)過強(qiáng)的“射線”效應(yīng)[41],又不至于讓計(jì)算量過大.參考文獻(xiàn)[52],相對密度變量梯度的無量綱化方式為方法的懲罰系數(shù)取經(jīng)典值 p=3.為了提高求解效率,對聲子BTE使用灰體近似,迭代過程中的靈敏度分析使用伴隨方法[53],優(yōu)化求解則采用構(gòu)造原問題局部近似模型的移動(dòng)漸進(jìn)方法(method of mo ving asymptotes,MMA)[54].收斂準(zhǔn)則設(shè)定為其中上標(biāo)j表示迭代次數(shù),認(rèn)為滿足此條件時(shí)設(shè)計(jì)變量、目標(biāo)函數(shù)值都達(dá)到穩(wěn)定.考慮到大部分區(qū)域?yàn)榛撞牧?用基底材料的平均自由程來計(jì)算努森數(shù),即 Kn=ls/a.

3.1 擴(kuò)散導(dǎo)熱的體點(diǎn)優(yōu)化

在進(jìn)行求解聲子BTE的拓?fù)鋬?yōu)化之前,先在傅里葉定律適用的條件下實(shí)現(xiàn)體點(diǎn)導(dǎo)熱的拓?fù)鋬?yōu)化,以檢驗(yàn)所發(fā)展的拓?fù)鋬?yōu)化方法.盡管凡是在實(shí)體和空材料間構(gòu)造連續(xù)函數(shù)并對中間密度值進(jìn)行懲罰的方法都可稱為SIMP方法,但通常都是對本構(gòu)方程的系數(shù)進(jìn)行插值[15].對基于傅里葉定律的拓?fù)鋬?yōu)化,本構(gòu)方程中與材料性質(zhì)有關(guān)的系數(shù)是熱導(dǎo)率,對應(yīng)SIMP方法的插值公式為k(ρ)=ks+(kf-ks)ρp,其中 kf是填充材料熱導(dǎo)率.此時(shí),填充材料和基底材料導(dǎo)熱能力的差異通過熱導(dǎo)率比kf/ks來體現(xiàn),kf/ks越大,傳導(dǎo)同樣熱量所需要的高導(dǎo)熱材料橫截面積越小,因而在相同的填充量下填充材料能夠更遠(yuǎn)地延伸到區(qū)域內(nèi)部,形成的樹狀結(jié)構(gòu)的“樹枝”更為細(xì)長,優(yōu)化效果也更好[12].對本文提出的優(yōu)化模型,本構(gòu)方程中與材料性質(zhì)有關(guān)的系數(shù)是聲子消光系數(shù) κ ,對應(yīng)的插值公式為(4)式.根據(jù)動(dòng)力學(xué)理論,熱導(dǎo)率和聲子平均自由程間的關(guān)系為其中C是材料比熱,v是聲子群g速度的大小.材料的比熱、群速度通常是確定的,于是有 kf/ks∝lf/ls,意味著不同材料導(dǎo)熱性能的差別主要源自聲子平均自由程的差別[55].此時(shí)(4)式中關(guān)于聲子消光系數(shù) κ (=l-1)的插值,對應(yīng)于傅里葉導(dǎo)熱定律下對熱導(dǎo)率倒數(shù) k-1的插值,而非傳統(tǒng)的對熱導(dǎo)率k的插值.在 kf/ks=500 ,α=0.05 ,n=80的條件下,分別對k和 k-1進(jìn)行插值,基于傅里葉導(dǎo)熱定律的拓?fù)鋬?yōu)化得到的結(jié)果如圖3所示.設(shè)定 kf/ks=500 是為了讓填充材料和基底材料的導(dǎo)熱能力存在較大差異,更直觀地體現(xiàn)優(yōu)化的效果,且實(shí)際電子器件中絕緣材料和導(dǎo)電材料的熱導(dǎo)率通常也相差2個(gè)數(shù)量級(jí)[1].文獻(xiàn)中關(guān)于傅里葉導(dǎo)熱定律下體點(diǎn)導(dǎo)熱拓?fù)鋬?yōu)化的研究也大多采用此熱導(dǎo)率比值.圖3左側(cè)表示 ρ 在區(qū)域內(nèi)的分布情況,白色代表基體材料(ρ=0),黑色表示填充材料(ρ=1),顏色越深意味著該位置的 ρ 值越接近于1;右側(cè)則是無量綱溫度 T?的分布.插值k的拓?fù)鋬?yōu)化得到了充分伸入內(nèi)部區(qū)域、含多個(gè)分岔和枝的樹狀填充材料分布,如圖3(a)所示,符合文獻(xiàn)[10-13]中的結(jié)果,對應(yīng)的無量綱溫度平均值是插值 k-1的優(yōu)化結(jié)果則由圖3(b)給出,填充材料集中在區(qū)域下半平面,與低溫邊界相接觸,向區(qū)域內(nèi)伸展的過程中只形成了2處分岔,此時(shí)與圖3(a)相比,圖3(b)中填充材料的主干更短、更寬,枝數(shù)量更少,對應(yīng)的平均溫度和最大溫度更高.在擴(kuò)散導(dǎo)熱條件下,不同形狀、尺寸的填充材料熱導(dǎo)率相等,具有相同的構(gòu)建高導(dǎo)熱通道的能力.為了降低溫度,填充材料要設(shè)法伸入內(nèi)部區(qū)域,縮短高導(dǎo)熱通道和區(qū)域內(nèi)熱源間的距離,所以圖3(a)的結(jié)果要比圖3(b)更好.受制于SIMP方法的插值方式,圖3(b)得到的是體點(diǎn)導(dǎo)熱優(yōu)化的某個(gè)局部最優(yōu)解,但填充材料仍然有向區(qū)域內(nèi)伸展的趨勢.圖3的結(jié)果表明,SIMP方法的插值方式會(huì)影響拓?fù)鋬?yōu)化的效果,對本構(gòu)方程的系數(shù)進(jìn)行插值是更合理的選擇.以傅里葉導(dǎo)熱定律為例,盡管插值 k-1的拓?fù)鋬?yōu)化也能收斂,但對應(yīng)的是優(yōu)化問題的某個(gè)局部最優(yōu)解,優(yōu)化效果不如插值k的結(jié)果好.類似地,當(dāng)導(dǎo)熱的本構(gòu)方程變成(2)式時(shí),SIMP方法應(yīng)當(dāng)對方程的系數(shù) κ 進(jìn)行插值.

圖3 擴(kuò)散導(dǎo)熱下不同插值方式對材料分布和溫度分布的影響 (a)插值k; (b)插值k-1Fig.3.The effect of interpolation ways on the material distributions and temperature distributions in diffusive heat conduction：(a) Interpolating k; (b) interpolating k-1.

圖4 擴(kuò)散導(dǎo)熱下拓?fù)鋬?yōu)化得到的材料分布隨 α 和n的變化Fig.4.The material distributions obtained by topology optimization in diffusive heat conduction varying with α and n.

除了SIMP方法的插值方式,α 的取值也很重要,因?yàn)?α 顯著響著拓?fù)鋬?yōu)化解的數(shù)值穩(wěn)定性[15].保持 kf/ks=500 不變,并使用對k的插值,擴(kuò)散導(dǎo)熱條件下拓?fù)鋬?yōu)化得到的填充材料分布隨 α 和n的變化如圖4所示.整體來看,填充材料的形狀都呈樹狀,與圖3(a)的結(jié)果大致相同,但細(xì)節(jié)上存在不少差異.當(dāng) α=0 時(shí),隨著網(wǎng)格密度的增大,高導(dǎo)熱通道的主干周圍會(huì)出現(xiàn)越來越多狹長的末梢,這是因?yàn)橥負(fù)鋬?yōu)化的數(shù)值解天然地存在著網(wǎng)格依賴性.當(dāng) α=0.05 時(shí),由于引入了對設(shè)計(jì)變量全局梯度的約束,拓?fù)鋬?yōu)化結(jié)果的網(wǎng)格依賴性得到有效抑制,網(wǎng)格細(xì)化并未導(dǎo)致太多狹長末梢.進(jìn)一步增大到 α=0.1 ,不同網(wǎng)格密度下拓?fù)鋬?yōu)化結(jié)果間的差異更小,網(wǎng)格無關(guān)性更好.但是,比較相同網(wǎng)格密度的拓?fù)鋬?yōu)化結(jié)果發(fā)現(xiàn),增大 α 的值后,中間密度值對應(yīng)的灰色區(qū)域增多,填充材料和基底材料之間的邊界變得模糊,這是控制設(shè)計(jì)變量全局梯度的代價(jià).綜合考慮后,本文的拓?fù)鋬?yōu)化選擇 α=0.05 ,n=80,此時(shí)所得結(jié)果具有較好的網(wǎng)格無關(guān)性和較清晰的材料邊界,且計(jì)算量相對較少.

3.2 彈道-擴(kuò)散導(dǎo)熱的體點(diǎn)優(yōu)化

根據(jù)上文確定的參數(shù)設(shè)置,對(5)式進(jìn)行數(shù)值迭代,獲得了體點(diǎn)導(dǎo)熱在彈道-擴(kuò)散導(dǎo)熱機(jī)制下的拓?fù)鋬?yōu)化解.設(shè)定 lf/ls=500 以便和傅里葉導(dǎo)熱定律下的結(jié)果對比,不同 Kn 下優(yōu)化得到的材料分布和溫度分布如圖5所示.總體來看,填充材料集中在區(qū)域下半平面,與低溫邊界充分接觸,其形狀完全不同于圖3(a)所示的伸入內(nèi)部區(qū)域且含多個(gè)枝的樹狀結(jié)構(gòu),而與圖3(b)的材料分布更為接近.如3.1節(jié)所述,本文采用的插值 κ 的方法,對應(yīng)于傅里葉定律下對 k-1的插值,因而優(yōu)化結(jié)果與插值k-1的結(jié)果相似.可以預(yù)見的是,在 Kn→0 時(shí),導(dǎo)熱過程接近純擴(kuò)散導(dǎo)熱,優(yōu)化會(huì)趨向于傅里葉定律下插值 k-1的結(jié)果.但受聲子彈道輸運(yùn)的影響,圖5與圖3(b)仍然存在明顯的差異.當(dāng) Kn=0.002 時(shí),低溫邊界處幾乎無溫度跳躍,區(qū)域內(nèi)無量綱溫度的最大值略高于圖3(b)中的0.23,表明大部分聲子以擴(kuò)散方式導(dǎo)熱,但仍有部分聲子發(fā)生了彈道輸運(yùn),這是因?yàn)樘畛洳牧系穆曌悠骄杂沙桃押蛥^(qū)域尺寸相當(dāng)(lf～a).此時(shí)填充材料在低溫邊界處連通,主干呈“U”型,兩側(cè)共有6處外凸的細(xì)枝.增大 Kn 至0.01,靠近低溫邊界處 T?=0.06 ,發(fā)生了一定的溫度跳躍,則是0.44,約為圖3(b)結(jié)果的2倍,表明彈道效應(yīng)增強(qiáng).此時(shí)填充材料分布的“U”形主干開口變大、寬度變厚,枝數(shù)量減少為2個(gè).當(dāng) Kn=0.1 時(shí),填充材料的體聲子平均自由程已比區(qū)域尺寸大一個(gè)數(shù)量級(jí),邊界處無量綱溫度的跳躍值接近0.6,區(qū)域內(nèi)峰值溫度幾乎是圖3(b)結(jié)果的6倍,彈道效應(yīng)已十分顯著.此時(shí),填充材料聚集在低溫邊界附近,除了頂部有凹陷外,其它方向的伸展長度基本一致,主干和枝合并形成團(tuán)狀分布.

為了檢驗(yàn)所提出的拓?fù)鋬?yōu)化方法的效果,分別在不同努森數(shù)下,對圖3和圖5中出現(xiàn)的5種材料分布進(jìn)行熱仿真,得到的值列于表1,括號(hào)內(nèi)的數(shù)據(jù)表示的值,即相對于材料均勻分布時(shí)的變化程度.表1顯示,不同努森數(shù)下,本文發(fā)展的拓?fù)鋬?yōu)化方法所對應(yīng)的值最低,如表中加粗?jǐn)?shù)據(jù)所示,表明所預(yù)測的材料分布方式確實(shí)是最優(yōu)的,驗(yàn)證了聲子BTE和拓?fù)鋬?yōu)化相結(jié)合的優(yōu)化方法的正確性.以 Kn=0.1 時(shí)為例,解傅里葉定律的拓?fù)鋬?yōu)化所得的材料分布,對應(yīng)的值大概為80%,而解聲子BTE的拓?fù)鋬?yōu)化得到的材料分布,對應(yīng)的值在70%左右,進(jìn)一步下降了10個(gè)百分點(diǎn),其中又以Kn=0.1條件下的拓?fù)鋬?yōu)化所得的構(gòu)型效果最好,此外,對相同的材料分布方式,隨著 Kn 的增大,的值增大,表明添加高導(dǎo)熱材料所能帶來的優(yōu)化效果變差.例如對圖3(a)所示的樹狀結(jié)構(gòu)(擴(kuò)散導(dǎo)熱條件下插值k的拓?fù)鋬?yōu)化結(jié)果),Kn=0.002 時(shí),是33.7%,但Kn=0.01時(shí)此值上升到53.2%,Kn=0.1 時(shí)只有86.0%.圖5和表1的結(jié)果共同表明,彈道-擴(kuò)散導(dǎo)熱機(jī)制下,體點(diǎn)導(dǎo)熱問題的最佳材料分布會(huì)隨著彈道效應(yīng)強(qiáng)度的變化而變化,基于傅里葉定律的拓?fù)鋬?yōu)化所得的材料分布不再是最優(yōu)的.因此,對微納米尺度下的散熱優(yōu)化問題,在設(shè)計(jì)階段必須要考慮聲子彈道輸運(yùn)的影響.值得一提的是,如果將現(xiàn)有優(yōu)化模型中對 κ 的插值改為對其倒數(shù)l的插值,優(yōu)化模型的收斂性會(huì)嚴(yán)重惡化,材料分布會(huì)出現(xiàn)嚴(yán)重的棋盤格問題和大量的中間密度值.

圖5 彈道-擴(kuò)散導(dǎo)熱下拓?fù)鋬?yōu)化的材料分布和溫度云圖隨 Kn 的變化 (a) Kn=0.002; (b) Kn=0.01; (c) Kn=0.1Fig.5.The topology optimization obtained material distributions and temperature distributions varying with Kn in ballistic-diffusive heat conduction：(a) Kn=0.002; (b) Kn=0.01; (c) Kn=0.1.

進(jìn)一步分析圖5中拓?fù)鋬?yōu)化得到的最優(yōu)材料分布,發(fā)現(xiàn)隨著 Kn 的增大,彈道效應(yīng)增強(qiáng),填充材料從低溫邊界向內(nèi)部區(qū)域發(fā)展的趨勢減弱,且分布形狀中微小的枝所占的比例減少,粗壯的主干所占的比例增多,這可以由彈道輸運(yùn)產(chǎn)生熱導(dǎo)率的尺寸效應(yīng)來解釋.不同于擴(kuò)散導(dǎo)熱,彈道-擴(kuò)散導(dǎo)熱機(jī)制下,結(jié)構(gòu)的等效熱導(dǎo)率會(huì)隨尺寸的減小而降低.因?yàn)槲⑿≈Φ某叽绫却謮阎鞲傻某叽绺?所以前者的等效熱導(dǎo)率比后者更低,前者在構(gòu)建高導(dǎo)熱通道方面的性能不如后者; 且隨著彈道效應(yīng)的增強(qiáng),這種能力上的差別會(huì)增大.在 Kn 較小時(shí),不同尺寸的填充材料間導(dǎo)熱性能差異較小,影響優(yōu)化效果的主要還是熱源和高導(dǎo)熱通道的距離,所以圖5(a)的結(jié)果像圖3(b)那樣有一些分岔,且伸入內(nèi)部區(qū)域較多.隨著 Kn 的增大,彈道效應(yīng)增強(qiáng),結(jié)構(gòu)等效熱導(dǎo)率的降低對優(yōu)化效果的影響加強(qiáng),拓?fù)鋬?yōu)化傾向于產(chǎn)生等效熱導(dǎo)率相對較大的粗壯結(jié)構(gòu),而非伸入內(nèi)部區(qū)域的微小分散結(jié)構(gòu).Kn=0.1 時(shí),盡管基底材料聲子平均自由程僅為系統(tǒng)特征尺寸的0.1倍,但填充材料的聲子平均自由程已經(jīng)是系統(tǒng)特征尺寸的50倍,整個(gè)系統(tǒng)的彈道效應(yīng)已十分顯著,團(tuán)狀的填充材料分布在構(gòu)建高導(dǎo)熱通道方面的性能最優(yōu),但由于填充材料比例限制,通道只能覆蓋低溫邊界附近的區(qū)域.彈道導(dǎo)熱極限下,拓?fù)鋬?yōu)化的結(jié)果中填充材料聚集在低溫邊界附近,沿各個(gè)方向的伸展長度基本相同,呈半圓形團(tuán)狀分布.圖5(c)的結(jié)果已經(jīng)較為接近彈道極限下的結(jié)果.溫度分布中低溫區(qū)域的輪廓反映了填充材料所形成的高導(dǎo)熱通道的效果.在圖3和圖5中,低溫區(qū)域的輪廓都與填充材料的形狀大體一致,表明所有填充材料都有效地形成了高導(dǎo)熱通道.根據(jù)上述分析,在彈道效應(yīng)較強(qiáng)時(shí),如果填充材料的形狀中含有較多的微小結(jié)構(gòu),由于這些小尺寸結(jié)構(gòu)在構(gòu)建高導(dǎo)熱通道方面是低效的,對應(yīng)區(qū)域的溫度將無法得到有效降低,溫度輪廓將勢必不同于填充材料的形狀.作為檢驗(yàn),給出擴(kuò)散導(dǎo)熱條件下插值 k-1的拓?fù)鋬?yōu)化所預(yù)測的材料分布(圖3(b))在 Kn=0.1 時(shí)的溫度場如圖6所示.盡管材料分布與圖3(b)相同,但圖6的溫度輪廓明顯不同于圖3(b),反而與圖5(c)的溫度輪廓更為接近,說明填充材料中粗壯的主干形成了高效的導(dǎo)熱通道,但細(xì)小的枝因?yàn)榈刃釋?dǎo)率更低,所以在導(dǎo)熱方面是相對低效的,未能有效降低所在區(qū)域的溫度.

表1 不同材料分布在不同努森數(shù)下對應(yīng)的無量綱溫度平均值Table 1.The average dimensionless temperature of different material distributions at different Knudsen numbers.

對體材料熱導(dǎo)率比值和界面熱阻對拓?fù)鋬?yōu)化結(jié)果的影響做簡單說明.盡管上述分析建立在kf/ks=500的條件下,但因?yàn)閺椀垒斶\(yùn)主要通過熱導(dǎo)率的尺寸效應(yīng)來影響拓?fù)鋬?yōu)化結(jié)果,所以改變kf/ks的值不影響體點(diǎn)導(dǎo)熱拓?fù)鋬?yōu)化結(jié)果隨彈道效應(yīng)強(qiáng)度的定性變化.定量上,體材料熱導(dǎo)率比的值越大,發(fā)生彈道輸運(yùn)時(shí)不同尺寸結(jié)構(gòu)等效熱導(dǎo)率的差異越大,相同努森數(shù)下彈道輸運(yùn)對拓?fù)鋬?yōu)化結(jié)果的影響越顯著.界面熱阻會(huì)使得材料間界面處出現(xiàn)溫度不連續(xù),導(dǎo)致系統(tǒng)總熱阻增大,平均溫度升高[4].雙層材料間界面熱阻的影響會(huì)隨著材料尺寸的減小而增大[56],如果在拓?fù)鋬?yōu)化中考慮界面熱阻,為了降低界面熱阻對導(dǎo)熱的阻礙作用,優(yōu)化結(jié)果會(huì)朝著增大結(jié)構(gòu)尺寸的方向發(fā)展.圖5中隨著彈道效應(yīng)的增強(qiáng),填充材料構(gòu)型也明顯表現(xiàn)出粗壯結(jié)構(gòu)增多、微小結(jié)構(gòu)減少的趨勢,這表明雖然現(xiàn)在的優(yōu)化模型未考慮界面熱阻,但所得的結(jié)果有利于抑制界面熱阻的作用.

圖6 擴(kuò)散導(dǎo)熱下插值 k-1 的拓?fù)鋬?yōu)化所得的材料分布在 Kn=0.1 時(shí)的溫度分布Fig.6.The temperature distribution for the material distribution obtained by topology optimization which interpolates k-1in diffusive heat conduction at Kn=0.1.

4 結(jié) 論

針對微納米電子器件的散熱優(yōu)化設(shè)計(jì),結(jié)合聲子玻爾茲曼方程和固體各向同性材料懲罰方法,發(fā)展了適用于微納尺度下彈道-擴(kuò)散導(dǎo)熱的拓?fù)鋬?yōu)化方法.聲子玻爾茲曼方程被轉(zhuǎn)化為聲子輻射輸運(yùn)方程,然后用離散坐標(biāo)法求解以獲得某個(gè)設(shè)計(jì)方案下的溫度場; 拓?fù)鋬?yōu)化則通過對聲子平均自由程的倒數(shù)進(jìn)行插值的固體各向同性材料懲罰方法來實(shí)現(xiàn),并采用對相對密度變量全局梯度施加顯示約束的方法來確保解的適定性和網(wǎng)格無關(guān)性.

成功實(shí)現(xiàn)了體點(diǎn)導(dǎo)熱問題在彈道-擴(kuò)散導(dǎo)熱機(jī)制下的拓?fù)鋬?yōu)化解,發(fā)現(xiàn)彈道-擴(kuò)散導(dǎo)熱機(jī)制下,拓?fù)鋬?yōu)化預(yù)測的最佳材料分布明顯不同于傅里葉定律下經(jīng)典的樹狀結(jié)構(gòu),且會(huì)隨努森數(shù)Kn的變化而變化,表明微納米尺度下的導(dǎo)熱優(yōu)化問題必須要考慮聲子彈道輸運(yùn)的影響.在Kn較小時(shí),解聲子玻爾茲曼方程的拓?fù)鋬?yōu)化的結(jié)果會(huì)趨向傅里葉定律下于對熱導(dǎo)率倒數(shù)進(jìn)行插值的拓?fù)鋬?yōu)化結(jié)果,而非傳統(tǒng)的插值熱導(dǎo)率所對應(yīng)的樹狀結(jié)構(gòu).隨著Kn的增大,拓?fù)鋬?yōu)化結(jié)果中填充材料逐漸向低溫邊界附近聚集,其構(gòu)型中粗壯的主干結(jié)構(gòu)所占比例增大,微小的枝結(jié)構(gòu)所占的比例降低,且優(yōu)化后系統(tǒng)平均溫度和優(yōu)化前系統(tǒng)平均溫度的比值增大.

拓?fù)鋬?yōu)化結(jié)果和Kn的相關(guān)性與聲子彈道輸運(yùn)產(chǎn)生的熱導(dǎo)率尺寸效應(yīng)有關(guān).擴(kuò)散導(dǎo)熱機(jī)制下,不同尺寸的填充材料熱導(dǎo)率相同,均能有效構(gòu)建高導(dǎo)熱通道; 但隨著彈道效應(yīng)的增強(qiáng),尺寸效應(yīng)使得材料分布中微小結(jié)構(gòu)的等熱導(dǎo)率比粗壯結(jié)構(gòu)低,微小結(jié)構(gòu)在提升系統(tǒng)導(dǎo)熱性能方面的作用不如粗壯結(jié)構(gòu),因而拓?fù)鋬?yōu)化會(huì)向增多粗壯結(jié)構(gòu)、減少微小結(jié)構(gòu)的方向發(fā)展.當(dāng)彈道效應(yīng)足夠強(qiáng)時(shí),填充材料會(huì)聚集在低溫邊界附近,主干和枝合并,呈團(tuán)狀分布.