基于蒙特卡洛法和分形理論的接觸導(dǎo)熱計(jì)算

張 琳,馮妍卉,王智學(xué)

(1.中國(guó)航發(fā)四川燃?xì)鉁u輪研究院， 成都 610500;2.北京科技大學(xué) 能源與環(huán)境工程學(xué)院， 北京 100083)

隨著傳熱技術(shù)的不斷發(fā)展和廣泛應(yīng)用，接觸導(dǎo)熱逐漸成為一個(gè)不可忽略的影響傳熱的因素，成為航空航天、機(jī)械制造、電子、低溫超導(dǎo)等學(xué)科和動(dòng)力工程領(lǐng)域中的研究熱點(diǎn)。接觸導(dǎo)熱是一個(gè)受材料熱物性、機(jī)械特性、表面形貌及接觸環(huán)境(壓力、溫度和間隙填充介質(zhì)等)等眾多因素影響的非線性問(wèn)題。在工程應(yīng)用中，通常需要降低或提高接觸熱阻以控制接觸傳熱，在設(shè)計(jì)過(guò)程中需要采用理論、實(shí)驗(yàn)或數(shù)值分析的方法對(duì)接觸熱導(dǎo)率或接觸熱阻進(jìn)行預(yù)測(cè)。但在不同的應(yīng)用中，材料種類繁多，應(yīng)用環(huán)境也不同，用實(shí)驗(yàn)手段來(lái)獲取接觸熱導(dǎo)率通常需要設(shè)計(jì)復(fù)雜的實(shí)驗(yàn)裝置，具有較高的成本。近年來(lái)，由于計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值分析的方法越來(lái)越多地應(yīng)用于接觸導(dǎo)熱的研究中[1-3]，但此類方法通常需要進(jìn)行大量的前、后數(shù)據(jù)處理工作，時(shí)間成本較高。因此，工程領(lǐng)域通常采用理論模型和相應(yīng)的計(jì)算公式對(duì)接觸熱導(dǎo)率進(jìn)行預(yù)測(cè)。國(guó)內(nèi)外學(xué)者對(duì)接觸導(dǎo)熱進(jìn)行了大量的研究，產(chǎn)生了許多具有重要參考價(jià)值和工程應(yīng)用價(jià)值的理論和方法。接觸導(dǎo)熱的理論模型和計(jì)算方法的研究包括對(duì)粗糙表面的形貌分析、接觸形變分析和接觸導(dǎo)熱分析，研究者將不同的粗糙表面形貌模型、形變模型和導(dǎo)熱模型進(jìn)行組合，形成了不同的接觸導(dǎo)熱計(jì)算模型。王安良等[4]對(duì)接觸導(dǎo)熱的理論預(yù)測(cè)方法進(jìn)行了較為全面的綜述。

早期的接觸導(dǎo)熱計(jì)算模型主要以統(tǒng)計(jì)學(xué)參數(shù)描述固體粗糙表面形貌，認(rèn)為固體粗糙表面是由無(wú)數(shù)個(gè)均勻分布的微凸體堆積而成，這些微凸體的高度符合高斯分布。以此為基礎(chǔ)，許多學(xué)者建立了接觸導(dǎo)熱的計(jì)算公式，如Mikic[5]、Yovanovich等[6]、陳劍楠等[7]。此外，Kumar等[8]建立了接觸熱阻的蒙特卡洛模型，該模型根據(jù)粗糙表面粗糙峰分布符合高斯分布的特性，在粗糙表面上隨機(jī)模擬N個(gè)符合高斯分布的粗糙峰，并在不同的接觸界面間距下計(jì)算出發(fā)生接觸的粗糙峰數(shù)目，每個(gè)發(fā)生接觸的粗糙峰都形成一個(gè)單點(diǎn)接觸熱阻，界面間總的接觸熱阻即為這些單點(diǎn)接觸熱阻的并聯(lián)。

需要指出的是，表面粗糙高度的分布并不完全符合高斯分布，以上所述接觸導(dǎo)熱計(jì)算模型都建立在以統(tǒng)計(jì)學(xué)參數(shù)描述的表面形貌模型的基礎(chǔ)上，而研究表明，表面形貌的統(tǒng)計(jì)學(xué)參數(shù)很明顯地受到測(cè)量?jī)x器的分辨率和取樣長(zhǎng)度的影響[9]，不能唯一確定地表征一個(gè)粗糙表面，在此基礎(chǔ)上建立的接觸導(dǎo)熱計(jì)算模型通常具有較大的不確定性，計(jì)算結(jié)果的準(zhǔn)確性與工程應(yīng)用需求差距較大。因此，研究者們開(kāi)始尋求與儀器分辨率和取樣長(zhǎng)度無(wú)關(guān)的粗糙表面表征參數(shù)，而分形理論使得粗糙表面的唯一表征成為可能。分形理論采用反映表面不規(guī)則程度的分形維數(shù)D和反映表面輪廓幅值的比例參數(shù)G來(lái)確定一個(gè)表面輪廓，分形參數(shù)D和G都是與儀器分辨率和取樣長(zhǎng)度無(wú)關(guān)的參數(shù)。Majumdar等[9]采用修正的Weierstrass-Mandelbrot函數(shù)(W-M函數(shù))描述粗糙表面輪廓，建立了固體接觸界面的彈塑性分形接觸模型和接觸導(dǎo)熱的M-T分形模型。Warren等[10]采用Cantor集分形形貌模型模擬粗糙表面，建立了彈性-塑性接觸和完全塑性接觸的Cantor集接觸模型和接觸導(dǎo)熱模型。此外，趙蘭萍[11]、徐瑞萍[12]、Zou等[13]、Ji等[14]、馬麗娜[15]和李小彭等[16]也進(jìn)行了基于表面分形理論的接觸導(dǎo)熱計(jì)算模型研究。

一系列研究表明，基于表面分形理論的接觸導(dǎo)熱計(jì)算模型，不受儀器分辨率和取樣長(zhǎng)度的影響，對(duì)接觸導(dǎo)熱的預(yù)測(cè)具有確定性，能夠較為準(zhǔn)確地預(yù)測(cè)接觸導(dǎo)熱。近年來(lái)將分形理論應(yīng)用于接觸導(dǎo)熱預(yù)測(cè)的研究較少，現(xiàn)有的分形理論和模型較為復(fù)雜，與實(shí)際工程應(yīng)用存在很大差距[17]。因此，本文提出了一種適合于工程應(yīng)用的接觸導(dǎo)熱分形計(jì)算模型。模型基于固體粗糙表面形貌的隨機(jī)特性和分形特性，采用蒙特卡洛法和W-M分形函數(shù)模擬表面形貌和接觸，構(gòu)建了固體界面間接觸熱阻的并聯(lián)網(wǎng)絡(luò)模型，簡(jiǎn)化了接觸熱阻的計(jì)算過(guò)程，模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果基本一致，能夠較為準(zhǔn)確地預(yù)測(cè)固體界面間的接觸導(dǎo)熱，具有較高的工程應(yīng)用價(jià)值。

1 固體粗糙表面形貌的分形表征

固體粗糙表面形貌是一個(gè)非穩(wěn)定隨機(jī)過(guò)程，如圖1所示[9]。表面輪廓曲線具有隨機(jī)、多尺度和無(wú)序的特性以及統(tǒng)計(jì)自仿射和自相似的數(shù)學(xué)特性，分形幾何學(xué)中的W-M函數(shù)可以滿足表面輪廓曲線的所有數(shù)學(xué)特征，其表達(dá)式如下：

(1)

式中：D是分形維數(shù)，它反映的是輪廓z(x)在所有空間尺度上的不規(guī)則性；G是決定輪廓高度幅值的比例系數(shù)，決定了z(x)的具體尺寸；γn為空間頻率，曲線的最低頻率取決于取樣長(zhǎng)度L，并有γn1=1/L。對(duì)于任意取樣長(zhǎng)度，決定輪廓高度z(x)的是分形參數(shù)D和G。

圖1 表面形貌的統(tǒng)計(jì)自仿射和自相似特性

由式(1)可以得出

z(γx)=γ(2-D)z(x)

(2)

若位置坐標(biāo)x放大γ倍，則輪廓高度放大為原來(lái)的γ(2-D)倍。粗糙表面的分形參數(shù)D和G是尺度獨(dú)立的，可以提供存在于分形表面上的所有尺度范圍內(nèi)的全部粗糙度信息[10]。

圖2為某粗糙表面實(shí)測(cè)形貌與采用分形分析后得到分形參數(shù)并采用W-M分形函數(shù)模擬的形貌對(duì)比，可以看出二者表面形貌曲線波動(dòng)幅度和頻率具有較高的符合性，W-M函數(shù)對(duì)固體粗糙表面的形貌具有較好的模擬效果。

圖2 W-M函數(shù)模擬形貌與真實(shí)表面形貌對(duì)比

2 基于蒙特卡洛法和分形理論的接觸熱阻計(jì)算

2.1 基本假設(shè)

本文研究在一定壓力下相接觸的固體表面間的接觸熱阻，基本假設(shè)如下：

1) 粗糙表面輪廓具有隨機(jī)性和分形特性，固體表面之間的接觸可等效為一個(gè)當(dāng)量粗糙表面和剛性光滑平面之間的接觸。

2) 固體表面間的接觸和傳熱僅僅發(fā)生在離散分布的一系列大小不同的接觸點(diǎn)上，接觸點(diǎn)發(fā)生塑性變形，忽略間隙介質(zhì)的導(dǎo)熱和熱對(duì)流以及非接觸部分的輻射換熱。

3) 表面間各接觸點(diǎn)的傳熱互不影響，每個(gè)接觸點(diǎn)相當(dāng)于一個(gè)熱流通道，并形成一個(gè)單點(diǎn)接觸熱阻，總的接觸熱阻是所有單點(diǎn)接觸熱阻的并聯(lián)。

2.2 接觸界面形貌的模擬

接觸界面形貌模擬的目的是確定粗糙表面接觸點(diǎn)的數(shù)量、尺寸及分布等參數(shù)。粗糙表面輪廓本質(zhì)上是一個(gè)隨機(jī)過(guò)程，因此本文將分形理論與蒙特卡洛法相結(jié)合，對(duì)表面粗糙峰的分布進(jìn)行模擬。

首先將2個(gè)粗糙表面等效成一個(gè)當(dāng)量粗糙表面，將粗糙表面的接觸等效為當(dāng)量粗糙表面與剛性光滑平面的接觸，從而求出當(dāng)量表面分形參數(shù)。一條數(shù)字化輪廓曲線可以視為一個(gè)時(shí)間序列，在z(x)曲線上取時(shí)間間隔為Δτ的N個(gè)采樣點(diǎn)，令τ=nΔτ，則該曲線的結(jié)構(gòu)函數(shù)為

E(τ)=〈[z(x+nΔτ)-z(x)]2〉=

(3)

等價(jià)于

E(τ)=CG2(D-1)τ4-2D

(4)

其中

(5)

式中函數(shù)Γ為Gamma函數(shù)，當(dāng)1

在雙對(duì)數(shù)坐標(biāo)上，結(jié)構(gòu)函數(shù)E(τ)與τ呈線性關(guān)系，即

lgE(τ)=B+klgτ

(6)

式中：

k=4-2D

(7)

B=lgCG2(D-1)

(8)

對(duì)于當(dāng)量粗糙表面，其結(jié)構(gòu)函數(shù)為2個(gè)粗糙表面結(jié)構(gòu)函數(shù)之和，即

E(τ)=E1(τ)+E2(τ)

(9)

當(dāng)量硬度H、當(dāng)量彈性模量E和當(dāng)量熱導(dǎo)率k通過(guò)以下各式求得：

H=min(H1,H2)

(10)

(11)

(12)

式中ν為材料的泊松比，下標(biāo)1和2分別代表相接觸的2個(gè)表面。

研究表明，表面粗糙峰均方根高度σ與分形參數(shù)具有如下關(guān)系[9]：

(13)

表面粗糙峰高度具有隨機(jī)性，并符合W-M函數(shù)，因此對(duì)于一個(gè)邊長(zhǎng)為L(zhǎng)的方形接觸表面，可隨機(jī)選取N個(gè)x值(0

N=L2η

(14)

式中η為表面粗糙峰密度。根據(jù)Hsieh[20]的結(jié)果，表面粗糙峰密度

(15)

式中m為表面粗糙峰平均斜率。

2.3 表面接觸參數(shù)的計(jì)算

本文采用的接觸點(diǎn)形變模型為塑性圓錐模型，如圖3所示。假設(shè)當(dāng)量粗糙表面與剛性光滑平面的距離為d，則高度為zi的粗糙峰的峰頂與剛性光滑平面的距離δi為

δi=zi-d

(16)

當(dāng)δi>0時(shí)，粗糙峰與剛性光滑平面發(fā)生接觸，其接觸半徑ai、接觸面積Aci和接觸載荷Fi分別為：

ai=δi/m

(17)

(18)

(19)

通過(guò)對(duì)所有粗糙峰的計(jì)算，可以得到發(fā)生接觸的粗糙峰的個(gè)數(shù)n，進(jìn)而計(jì)算出表面距離為d時(shí)的實(shí)際接觸面積Ar和總接觸壓力Fc：

(20)

(21)

上述計(jì)算過(guò)程是在已知d的情況下求界面間接觸載荷Fc等參數(shù)，然而在進(jìn)行實(shí)驗(yàn)研究或?qū)嶋H工程應(yīng)用中很難知道接觸界面的實(shí)際距離，但很容易知道接觸載荷Fc，因此在接觸熱阻的理論計(jì)算中可迭代求得d，進(jìn)而求出其他參數(shù)。

圖3 塑性圓錐形變模型

2.4 單點(diǎn)接觸熱阻和總熱阻的計(jì)算

研究表明，半無(wú)限長(zhǎng)流管換熱模型更適用于接觸熱阻的計(jì)算，Yovanovich[6]采用的單點(diǎn)接觸熱阻計(jì)算公式為

(22)

式中

ψ=(1-ε)1.5

(23)

并且有

(24)

式(24)中Aa為名義接觸面積，

Aa=L2

(25)

根據(jù)本文的基本假設(shè)可知，總熱阻為各單點(diǎn)接觸熱阻的并聯(lián)，因此可以得出總接觸熱阻Rc和接觸熱導(dǎo)率hc：

(26)

(27)

3 結(jié)果與分析

3.1 計(jì)算方法的驗(yàn)證

本文采用文獻(xiàn)[11]中的2組實(shí)驗(yàn)數(shù)據(jù)對(duì)接觸熱阻的計(jì)算模型進(jìn)行了驗(yàn)證，并與Mikic模型、Yovanovich模型和MT分形模型進(jìn)行了對(duì)比。實(shí)驗(yàn)材料為鋁合金Al5052，實(shí)驗(yàn)溫度為155 K，表1和表2分別為材料Al5052的機(jī)械性能和熱力學(xué)性能參數(shù)以及4個(gè)Al5052試樣的表面分形參數(shù)。

圖4和圖5分別為理論計(jì)算結(jié)果與文獻(xiàn)[11]中的實(shí)驗(yàn)結(jié)果的對(duì)比，可以看出：Yovanovich模型的接觸熱導(dǎo)率計(jì)算值高于實(shí)驗(yàn)值，并有很大偏差，而Mikic模型的接觸熱導(dǎo)率計(jì)算值小于實(shí)驗(yàn)值。本文的分形模型和MT分形模型的接觸熱導(dǎo)率計(jì)算值與實(shí)驗(yàn)值比較接近，這表明分形模型能夠較為準(zhǔn)確地預(yù)測(cè)接觸熱阻。在Yovanovich模型中，認(rèn)為粗糙表面微凸體的形變是塑形形變，但并沒(méi)有考慮變形后體積的損失，所以在計(jì)算過(guò)程中，每對(duì)微凸體接觸點(diǎn)的接觸半徑都要小于實(shí)際變形后的接觸半徑，接觸點(diǎn)的對(duì)數(shù)就會(huì)相應(yīng)地增加，因此得出的接觸熱導(dǎo)率計(jì)算值偏大。Mikic模型認(rèn)為微凸體的變形為彈性形變，而在實(shí)際接觸過(guò)程中，隨著壓力的增大，必然會(huì)引起塑形形變，所以接觸熱導(dǎo)率求解的結(jié)果會(huì)偏低。

表1 鋁試樣Al5052的熱力學(xué)和機(jī)械性能

表2 鋁試樣Al5052的表面分形參數(shù)

圖4 接觸熱導(dǎo)率理論計(jì)算值與實(shí)驗(yàn)值對(duì)比(試樣1/試樣2)

圖5 接觸熱導(dǎo)率理論計(jì)算值與實(shí)驗(yàn)值對(duì)比(試樣3/試樣4)

表面形貌是接觸導(dǎo)熱的主要影響因素，傳統(tǒng)的粗糙度參數(shù)受到儀器分辨率和取樣長(zhǎng)度的影響，使用尺度獨(dú)立的分形參數(shù)模擬粗糙表面形貌建立的接觸熱阻計(jì)算模型更有利于對(duì)接觸傳熱的預(yù)測(cè)。

3.2 分形參數(shù)對(duì)接觸熱阻的影響分析

計(jì)算模擬了2個(gè)具有相同性能的尺寸為2.5 mm×2.5 mm的材料表面的接觸(材料參數(shù)見(jiàn)表3)，計(jì)算中設(shè)定的分形參數(shù)為當(dāng)量分形參數(shù)。

表3 模擬材料的熱力學(xué)和機(jī)械性能

圖6為G=1.0×10-5m時(shí)，不同壓力下的接觸熱導(dǎo)率與分形維數(shù)D的關(guān)系?？梢钥闯觯簤毫υ酱蠼缑骈g的接觸熱導(dǎo)率越大，這反映了壓力增大使界面間的接觸面積增大。在定壓下看，接觸熱導(dǎo)率隨著分形維數(shù)D的增大而增大，這是因?yàn)楫?dāng)D增大時(shí)，表面輪廓的精細(xì)結(jié)構(gòu)增多，精細(xì)結(jié)構(gòu)增加導(dǎo)致整個(gè)輪廓的整體高度降低，減低了輪廓粗糙度，導(dǎo)致接觸導(dǎo)熱的增強(qiáng)，壓力越大接觸導(dǎo)熱的增強(qiáng)越明顯。

圖7為在接觸壓力p=2 MPa時(shí)、不同分形參數(shù)G下，接觸熱導(dǎo)率與分形維數(shù)D的關(guān)系?？梢钥闯觯航佑|熱導(dǎo)率隨分形維數(shù)D的增大而增大，分形比例參數(shù)G越小，增大越明顯。圖8為接觸壓力p=2 MPa時(shí)、不同分形維數(shù)D下接觸熱導(dǎo)率與分形參數(shù)G的關(guān)系，可以看出：分形維數(shù)D一定時(shí)，分形參數(shù)G越大，接觸熱導(dǎo)率越小。這是因?yàn)榉中螀?shù)G越大，表面粗糙度越大，增加了界面間的傳熱阻力。分形維數(shù)D越大，G的變化對(duì)接觸熱導(dǎo)率的影響越明顯。

圖6 不同接觸壓力下分形維數(shù)D對(duì)接觸熱導(dǎo)率的影響(G=1.0×10-5 m)

圖7 不同分形比例參數(shù)G條件下分形維數(shù)D對(duì)接觸熱導(dǎo)率的影響(P=2 MPa)

圖8 不同分形維數(shù)D的條件下分形比例參數(shù)G對(duì)接觸熱導(dǎo)率的影響(P=2 MPa)

綜上所述，在一定的接觸壓力下，接觸界面間的熱導(dǎo)率隨著分形維數(shù)D的增大而增大，隨著分形參數(shù)G的增大而減小。這是因?yàn)樵谙嗤姆中螀?shù)G下，D越大表面粗糙度越小，在相同的分形參數(shù)D下，G越小表面粗糙度越小。表面粗糙度越小，材料表面越光滑，表面之間的接觸越充分，接觸熱導(dǎo)率越大。

4 結(jié)束語(yǔ)

本文采用蒙特卡洛法和W-M分形函數(shù)模擬表面形貌，建立了一種能夠反映固體粗糙表面形貌的隨機(jī)特性和分形特性的固體接觸界面間接觸導(dǎo)熱的計(jì)算模型，其計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果基本一致，能夠較為準(zhǔn)確地預(yù)測(cè)固體界面間的接觸導(dǎo)熱。分析了分形參數(shù)對(duì)固體界面間接觸導(dǎo)熱的影響。分析結(jié)果表明：在一定的接觸壓力下，分形參數(shù)G一定時(shí)，接觸熱導(dǎo)率隨著分形維數(shù)D的增大而增大，G越小，增大越明顯；分形參數(shù)D一定時(shí)，接觸熱導(dǎo)率隨著分形參數(shù)G的增大而減小，分形維數(shù)D越大，G的變化對(duì)接觸熱導(dǎo)率的影響越明顯。由于分形參數(shù)能夠唯一確定地表征固體粗糙表面的形貌，本文提出的模型對(duì)接觸熱阻的預(yù)測(cè)具有確定性，為接觸傳熱的工程設(shè)計(jì)提供了一種有效的方法。