繆紅林，張 良

（浙江大學(xué) 熱工與動力系統(tǒng)研究所，浙江 杭州 310027）

0 引 言

接觸熱阻作為制約傳熱過程的關(guān)鍵熱阻之一，在航天[1]、微電子[2-3]、超導(dǎo)[4]、能源[5]等領(lǐng)域的熱設(shè)計和熱管理中具有重要影響。 準(zhǔn)確預(yù)測和分析導(dǎo)熱材料間的接觸熱阻形成機理和規(guī)律，并實現(xiàn)對接觸熱阻特征的主動控制，對于傳熱器件散熱結(jié)構(gòu)的優(yōu)化設(shè)計，提升傳熱性能和能效水平具有重要的意義。

現(xiàn)有研究表明，材料的性質(zhì)，接觸面的形貌特征，間隙介質(zhì)、接觸面載荷及環(huán)境等是影響接觸熱阻的主要因素[6]。 目前大量有關(guān)接觸熱阻的理論研究[7-8]和實驗研究[9-12]工作得以開展，其中，實驗研究主要圍繞具體材料在不同溫度、壓力條件下的接觸熱阻特征進行測量。 實驗測量法大致可以分為適用于測量較厚材料的穩(wěn)態(tài)測量法[9]和適合測量薄膜材料的瞬態(tài)測量法[10]。 如馮飆等[11-12]提出的改進版穩(wěn)態(tài)測量法，可用于測量亞毫米厚度材料間的接觸熱阻，并將測量誤差控制在5%以內(nèi)。

相較于實驗研究，部分研究者致力于采用數(shù)值方法對不同因素的影響規(guī)律進行更為系統(tǒng)和深入的分析，以期望實現(xiàn)接觸熱阻的準(zhǔn)確預(yù)測：Xiang等[13]采用非平衡的分子動力學(xué)方法模擬接觸熱阻，Gou 等[14]將材料真實物理形貌導(dǎo)入ANSYS 軟件，使用有限元方法求解接觸熱阻。 其中，在數(shù)值模擬方法中，由于介觀尺度的格子Boltzmann 方法在處理具有復(fù)雜邊界的傳熱模型上具有一定的優(yōu)勢[15]，因此采用該方法對接觸截面上的具有復(fù)雜形貌特征的接觸熱阻的計算更為有效。 在基于格子Boltzmann 方法的研究中，根據(jù)模型中熱阻值是否已知，可以將之分為兩類。 一種是在已知接觸熱阻值的情況下，在模型的接觸面處添加接觸熱阻影響的部分反彈法[16-17]和粒子圖像法[18]，這兩種方法都可以在傳熱模型中便捷地添加接觸熱阻的影響。 另一種則是考慮到接觸面形貌，形變模型和傳熱模型的建模方法：崔騰飛等[19]使用多尺度方法研究高斯表面間的接觸熱阻，在接觸面加密網(wǎng)格并使用格子Boltzmann方法求解傳熱問題，非接觸面則使用有限差分法加速計算；方文振等[20]提出的多重網(wǎng)格模擬方法，在接觸面加密網(wǎng)格捕捉平直面上的分形形貌，研究平直面間接觸熱阻的變化規(guī)律。 與此同時，當(dāng)前有關(guān)接觸熱阻的研究主要針對的是常見的平面上的接觸熱阻，而對于曲面間的接觸熱阻的研究鮮有報道，而曲面熱阻在以內(nèi)嵌管式換熱器[21]為代表的傳熱元件中普遍存在。 準(zhǔn)確預(yù)測固體曲面間的接觸熱阻，理清曲面熱阻與平面熱阻的共性規(guī)律與差異，對于相關(guān)傳熱元件的傳熱過程理解與優(yōu)化具有重要的科學(xué)和工程意義。

因此，本文建立基于格子Boltzmann 方法的單松弛數(shù)值模型對固-固曲面間的接觸熱阻特征進行數(shù)值計算，并研究了溫度、壓力、材料和表面粗糙度四種因素對曲面接觸熱阻特征的影響。

1 數(shù)值模型

1.1 物理模型

本文采用格子Boltzmann 方法對兩個緊密接觸的等壁厚圓環(huán)在恒溫邊界下（Tin，Tout）的傳熱過程模擬計算來研究固—固曲面間的接觸熱阻特性。 其中，外圓半徑rout為2000 μm，內(nèi)圓半徑rin為1500 μm，兩個圓環(huán)不相互接觸壁面設(shè)為絕對光滑。 該問題的計算域可以簡化為如圖1 所示的四分之一緊密接觸的扇形，扇形內(nèi)外壁面為恒溫，兩端為絕熱。

圖1 接觸熱阻模型示意圖

兩個圓環(huán)相接觸的圓弧表面形貌特征采用基于分形理論的W-M函數(shù)[22]進行描述，使用曲面上均勻遞增的弧長作為W-M函數(shù)的參數(shù)：

式中：z為分形形貌在曲面法線方向偏離曲面的高度；l為曲面均勻遞增的弧長；G為影響特征曲面高度的分形特征尺度系數(shù)；D為影響表面復(fù)雜度的分形維數(shù)；γ為常數(shù)（對于正態(tài)分布的隨機輪廓，一般取γ=1.5）；n 為自然序列。本文中選取D=1.4，G=1.0 ×10-9，n =20 ～40 生成表面形貌[23]。

接觸面的粗糙度使用算術(shù)平均高度Ra表征：

式中：L 為積分總弧長；l為積分微元；z為表面高度。

為了進一步考慮接觸面上的形變特征，在接觸點的形變模型上，使用CMYTCR接觸模型[24]進行描述。 該模型假設(shè)接觸形變是純塑性形變，接觸點的受力平衡方程為：

式中：Ar為實際接觸面積；Aa為名義接觸總面積；Hc為材料的微觀硬度。 在模型中由于實際接觸面積占總面積之比很少，因此不考慮塑性形變部分的體積變化。 故壓力的計算方法：

1.2 格子Boltzmann 方法求解

根據(jù)上述物理模型及邊界條件特征，其接觸熱阻Rtcr的表達式為：

式中：?為曲面上的熱通量（采用二維模管，熱通量單位為W/m）。 為了求得熱通量，需要建立格子Boltzmann 方法對如下二維熱擴散方程進行求解：

式中：T為溫度；t為時間；λ為材料的導(dǎo)熱系數(shù)；ρCp為體積比熱。

本文采用單松弛的格子Boltzmann 方法求解溫度擴散問題的演化方程[25]如下：

如圖2 所示，速度離散使用D2Q5[26]的離散方法。 離散后的速度矢量c為：

圖2 D2Q5 模型

式中：q1與q2為沿分布函數(shù)遷移方向的已知熱流密度，γ和β 分別為兩已知熱流密度與法向熱流密度之間的夾角。

圖3 熱流密度示意圖

在接觸面間的傳熱過程中，以分形曲線為邊界，將空間中的點劃分為固體材料的點與間隙介質(zhì)的點，熱流需要經(jīng)過固體材料，流經(jīng)間隙介質(zhì)到另一邊的固體材料。 為保證溫度與熱流密度在耦合界面連續(xù)，需假設(shè)間隙介質(zhì)與固體材料的體積比熱相等[27]，即：

其中，（ρCp）f為間隙介質(zhì)的體積比熱，（ρCp）s為固體材料的體積比熱。 上述假設(shè)，并不影響最終穩(wěn)態(tài)的溫度場。

1.3 曲面邊界的處理

圖4 曲邊邊界示意圖

2 模型驗證

為了驗證本文建立的格子Boltzmann 模型計算曲面上接觸熱阻的可靠性，本文通過對圖1 中的表面粗糙度進行理想光滑處理后進行模擬計算。 其中接觸間隙取20 μm，外壁面與內(nèi)壁面溫度分別取303 K和293 K。 材料和間隙介質(zhì)的導(dǎo)熱系數(shù)分別取210 W/（m·K）和10 W/（m·K）。其中，熱阻理論解可由以下公式直接求得：

不同網(wǎng)格分辨率的計算結(jié)果與理論解之間的誤差如表1 所示。 而且，隨著網(wǎng)格分辨率提高，網(wǎng)格對計算結(jié)果的準(zhǔn)確性影響遞減。 表明本文建立的格子Boltzmann 模型能夠準(zhǔn)確求解曲面邊界的傳熱問題。

表1 理想光滑間隙下不同網(wǎng)格分辨率的計算模型準(zhǔn)確性驗證

與此同時，如表1 所示，隨著網(wǎng)格分辨率的增加，計算誤差越小。 因此，為了進一步驗證實際模型中的網(wǎng)格無關(guān)性，本文對圖1 中初始條件下（選取D=1.4，G=1.0 ×10-9，n =20 ～40 生成表面形貌）的粗糙表面下的接觸熱阻計算的網(wǎng)格無關(guān)性進行驗證，固體材料為內(nèi)外壁溫分別為293 K和303 K的不銹鋼圓環(huán)，其結(jié)果如表2 所示。

表2 分形形貌下網(wǎng)格分辨率測試

從表2 中可以看出，當(dāng)網(wǎng)格分辨率從4 μm降低至2 μm時，接觸熱阻提高了9.06%。 而當(dāng)分辨率從1 μm降低至0.5 μm時，熱阻值僅提高0.3%。 網(wǎng)格對計算結(jié)構(gòu)的影響可以忽略。 考慮計算效率，最優(yōu)網(wǎng)格分辨率為1 μm。

3 結(jié)果與討論

3.1 溫度的影響

如圖5 所示為，內(nèi)外壁溫分別為293 K和303 K不銹鋼圓環(huán)，選取D=1.4，G=1.0 ×10-9，n =20 ～40生成粗糙接觸面，在18.3 MPa壓力下，穩(wěn)態(tài)條件下得到的溫度分布云圖。 通過局部放大可以看出，在接觸間隙附近存在較大的溫度梯度。 這是由于間隙介質(zhì)空氣的導(dǎo)熱系數(shù)遠小于接觸材料，是主要傳熱熱阻，顯著抑制接觸界面間的傳熱。

圖5 接觸熱阻的溫度分布云圖

圖6 為不銹鋼接觸熱阻值隨溫度（上下壁面溫度的平均值，其中上下壁面溫差為10 K）的變化規(guī)律，接觸熱阻值單調(diào)降低，這是因為間隙介質(zhì)空氣的導(dǎo)熱系數(shù)隨著溫度的升高而升高。

圖6 接觸熱阻隨溫度的變化

上述模擬結(jié)果表明，接觸界面區(qū)域的傳熱過程主要受間隙介質(zhì)的導(dǎo)熱系數(shù)抑制。

3.2 接觸壓力及接觸材料的影響

如圖7 所示，分別計算了鋁和不銹鋼兩種固體材料環(huán)形曲面在特征粗糙度為1.81 μm（選取D=1.4，G=1.0 ×10-9，n =20 ～40 生成表面形貌）的具有分形特征的接觸表面下接觸熱阻隨壓力的變化特征。 從圖中可以看出，兩種材料的接觸熱阻隨壓力的變化規(guī)律表現(xiàn)一致，即隨著壓力的增大接觸熱阻呈現(xiàn)出降低的趨勢，且其降低速率逐漸趨于緩和。

圖7 接觸熱阻隨壓力的變化

為了進一步定量分析壓力對接觸面積影響，本文對不銹鋼材料的接觸面上間隙介質(zhì)節(jié)點數(shù)隨壓力的變化特征進行了統(tǒng)計分析，其結(jié)果如表3所示。 當(dāng)壓力從0.8 MPa提高至13 MPa時，間隙介質(zhì)點平均每兆帕減少898.52 個。 而壓力從40 MPa提高到88 MPa，間隙介質(zhì)點平均每兆帕減少57.29 個。 對應(yīng)的，其接觸面的接觸特征如圖8所示，可以看出間隙寬度隨著壓力升高明顯減小。這表明隨著接觸面間的壓力增加，微觀接觸面積會增加，且其增加速率會逐漸變小，壓力對接觸熱阻的影響變得不再明顯。 這一結(jié)果與平面下的結(jié)論較為一致[21]。 此外，由于不同材料的硬度特征不一樣，導(dǎo)致在同一壓力下其接觸面積不一致，從而導(dǎo)致其接觸熱阻的絕對值上存在差異。 如圖7所示，高硬度的不銹鋼，在同樣壓力工況下，接觸熱阻略大于鋁。

圖8 不同壓力下粗糙表面的接觸情況

表3 不同壓力下間隙介質(zhì)厚度量化分析

3.3 表面粗糙度的影響

為研究表面粗糙度對曲面間接觸熱阻的影響規(guī)律，本文選取內(nèi)外壁面溫度為293 K與303 K，間隙介質(zhì)為空氣，接觸材料為不銹鋼，選取D=1.4，n =20 ～40 生成分形表面形貌，并通過改變分形特征尺度系數(shù)G來控制材料表面粗糙度，分別計算幾種粗糙度下不同壓力工況的接觸熱阻，具體模擬參數(shù)如表4 所示。

表4 分形特征尺度系數(shù)與表面粗糙度

接觸熱阻隨表面粗糙度的變化規(guī)律如圖9 所示。 對于三種不同形貌的表面，接觸熱阻均隨壓力的增加而降低，并且降低速率變小，這與3.2 節(jié)的模擬結(jié)果吻合。 與此同時，接觸熱阻值在各壓力工況下，均隨表面粗糙度的增加而增加，可見高表面粗糙度是接觸熱阻存在的重要原因。

圖9 粗糙度對接觸熱阻的影響

4 結(jié) 論

（1）本文提出一種基于格子Boltzmann 方法預(yù)測曲面間接觸熱阻的模型。 本模型使用曲面坐標(biāo)作為W-M 函數(shù)的參數(shù)生成接觸面分形形貌，使用塑性形變模型分析接觸面的形變特征，使用插值方法處理曲面邊界。 在計算有解析解的理想光滑曲面間隙熱阻模型時，得到較為準(zhǔn)確的計算結(jié)果。

（2）研究結(jié)果表明：不同影響因子主要通過改變接觸面的間隙特征與間隙介質(zhì)的導(dǎo)熱系數(shù)最終影響接觸熱阻的大小。 其中溫度主要影響間隙介質(zhì)的導(dǎo)熱系數(shù)，壓力、材料和粗糙度影響接觸面的接觸特征。

（3）當(dāng)間隙介質(zhì)為空氣時，接觸熱阻隨溫度的升高緩緩降低。 高壓力，低硬度且表面光滑的材料，接觸更充分，接觸面間隙更小，其接觸熱阻往往更小。