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

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

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

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

需要指出的是，表面粗糙高度的分布并不完全符合高斯分布，以上所述接觸導(dǎo)熱計算模型都建立在以統(tǒng)計學(xué)參數(shù)描述的表面形貌模型的基礎(chǔ)上，而研究表明，表面形貌的統(tǒng)計學(xué)參數(shù)很明顯地受到測量儀器的分辨率和取樣長度的影響[9]，不能唯一確定地表征一個粗糙表面，在此基礎(chǔ)上建立的接觸導(dǎo)熱計算模型通常具有較大的不確定性，計算結(jié)果的準確性與工程應(yīng)用需求差距較大。因此，研究者們開始尋求與儀器分辨率和取樣長度無關(guān)的粗糙表面表征參數(shù)，而分形理論使得粗糙表面的唯一表征成為可能。分形理論采用反映表面不規(guī)則程度的分形維數(shù)D和反映表面輪廓幅值的比例參數(shù)G來確定一個表面輪廓，分形參數(shù)D和G都是與儀器分辨率和取樣長度無關(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]也進行了基于表面分形理論的接觸導(dǎo)熱計算模型研究。

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

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

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

(1)

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

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

由式(1)可以得出

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

(2)

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

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

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

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

2.1 基本假設(shè)

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

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

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

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

2.2 接觸界面形貌的模擬

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

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

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

(3)

等價于

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

(4)

其中

(5)

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

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

lgE(τ)=B+klgτ

(6)

式中：

k=4-2D

(7)

B=lgCG2(D-1)

(8)

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

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

(9)

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

H=min(H1,H2)

(10)

(11)

(12)

式中ν為材料的泊松比，下標1和2分別代表相接觸的2個表面。

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

(13)

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

N=L2η

(14)

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

(15)

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

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

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

δi=zi-d

(16)

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

ai=δi/m

(17)

(18)

(19)

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

(20)

(21)

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

圖3 塑性圓錐形變模型

2.4 單點接觸熱阻和總熱阻的計算

研究表明，半無限長流管換熱模型更適用于接觸熱阻的計算，Yovanovich[6]采用的單點接觸熱阻計算公式為

(22)

式中

ψ=(1-ε)1.5

(23)

并且有

(24)

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

Aa=L2

(25)

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

(26)

(27)

3 結(jié)果與分析

3.1 計算方法的驗證

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

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

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

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

圖4 接觸熱導(dǎo)率理論計算值與實驗值對比(試樣1/試樣2)

圖5 接觸熱導(dǎo)率理論計算值與實驗值對比(試樣3/試樣4)

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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