溫晶晶, 李磊濤, 劉承騖, 吳斌

(西北工業(yè)大學(xué) 航天學(xué)院, 陜西 西安 710072)

兩固體表面的接觸實際只發(fā)生在一些離散點或微小面上。當熱量從一個界面向另一個界面?zhèn)鲗?dǎo)時，接觸面附近的溫度會發(fā)生突變。把這種由界面不一致接觸引起的接觸換熱的附加阻力稱為接觸熱阻[1]。接觸熱阻研究已經(jīng)深入到機械制造、微電子、航空航天等諸多工程領(lǐng)域[2]。特別是航天領(lǐng)域：在真空、高低溫交替的環(huán)境下，航天器部件之間的接觸熱阻對傳熱效率起著至關(guān)重要的作用，進而直接影響航天器的工作狀態(tài)和壽命[3]。

目前，關(guān)于接觸熱阻計算方法的研究主要包括理論研究[4]、仿真計算[5-6]及試驗三方面。考慮到接觸傳熱機理復(fù)雜，影響因素多，個別物性參數(shù)定義不準確等原因，試驗依然是求解接觸熱阻不可或缺的方法。而這種根據(jù)試驗測得溫度場分布，結(jié)合數(shù)據(jù)處理方法辨識出接觸熱阻的方法本質(zhì)上屬于熱傳導(dǎo)反問題范疇。其中，李建平等人[7]通過試驗測得鋁卷之間溫度分布，再結(jié)合差分方法辨識出了鋁卷間的接觸熱阻；劉冬歡等人[8]通過試驗測得C/C復(fù)合材料和GH600高溫合金之間的軸向溫度值和熱流密度，利用Fourier定律外推出接觸界面溫差，進而求得接觸熱阻；石友安等人[9]將靈敏度法和共軛梯度法應(yīng)用到試驗中，對相變材料之間的接觸熱阻進行了反辨識研究；El-Sabbagh、Gill等人[10-11]也分別利用共軛梯度法和靈敏度法對接觸熱阻進行了反辨識研究；王超[12]將順序函數(shù)法應(yīng)用到試驗中，對高強度鋼熱成形過程中板料與模具之間的接觸熱阻進行了反辨識研究?？偨Y(jié)來看，目前國內(nèi)外關(guān)于接觸熱阻的反辨識研究都是將問題簡化為兩層介質(zhì)之間的一維導(dǎo)熱問題，并利用參數(shù)辨識的方法對接觸熱阻進行求解。該方法存在以下不足：①只能辨識出接觸熱阻的單一值，然而接觸面附近的實際傳熱方式是三維的，真實的接觸熱阻是隨接觸面位置變化而變化的，特別在接觸區(qū)域較大或接觸面拓撲形貌變化較大的情況下必須考慮接觸熱阻的空間變化；②溫度測點必須放置在介質(zhì)內(nèi)部軸線上，以最大限度模擬出一維溫度場梯度變化，這就要求熱電偶必須嵌入試件內(nèi)部，而這會給溫度測點的選擇帶來不便，同時也會破壞試件本身的溫度場并帶來測量誤差，并且很難應(yīng)用于尺寸較小的試件的測試[13]。

針對上述問題，本文采用邊界元法對溫度場進行離散求解，采用共軛梯度法進行參數(shù)辨識，二者相結(jié)合對接觸熱阻進行反辨識研究。本方法將問題抽象為兩層介質(zhì)之間的二維導(dǎo)熱問題，可以求出隨界面接觸線位置變化的接觸熱阻；充分發(fā)揮邊界元算法只需計算邊界溫度與熱流不需計算內(nèi)部溫度與熱流的優(yōu)勢，給溫度測點的選取帶來了很大的便利；算例分析表明：該方法可以有效求解出接觸熱阻，但作為接觸熱阻反辨識方法的一種，同樣存在對溫度測量誤差敏感的不適定性問題，并且溫度測點距離接觸界面越遠誤差越大，采用最小二乘法對計算結(jié)果進行優(yōu)化可以有效提高辨識的精度和穩(wěn)定性。

1 二維熱傳導(dǎo)反問題的計算模型

1.1 建立二維熱傳導(dǎo)模型

建立如圖1所示的二維熱傳導(dǎo)模型：導(dǎo)熱系數(shù)分別為k1和k2的2塊材料在y=l處接觸;在接觸界面附近選取一系列溫度測點并測取溫度值;大多數(shù)研究忽略了接觸熱阻在接觸區(qū)域內(nèi)沿接觸面空間的變化,本文假設(shè)接觸熱阻在二維結(jié)構(gòu)中沿接觸線位置變化,記為R(x)。建立穩(wěn)態(tài)熱傳導(dǎo)控制方程為:

圖1 二維熱傳導(dǎo)模型

(1)

在不考慮熱對流和熱輻射的情況下,由第一類邊界條件得:

(1)式～(5)式中，T=T(x,y)為溫度變量;?T/?x=0為兩側(cè)(x=L或x=0處)的邊界絕熱條件;Tu,Td分別為上、下表面(y=L、y=0處)的給定溫度值;Ti(x)為y=Li處的溫度函數(shù)值,可由上述方程組解出。

1.2 建立目標函數(shù)和迭代停止標準

(6)

式中,Tl(x)為假定的接觸界面的溫度(以下簡稱界面溫度);Ti(x)是根據(jù)T(x)由控制方程(1)計算得到的測點處的溫度函數(shù),Ti(xj)為其離散溫度值。

用迭代法搜索求解界面溫度Tl(x),并規(guī)定迭代停止標準為:

(7)

若考慮溫度測量誤差,則ε也可以寫成:

(8)

式中,σ為溫度測量的方均根誤差。

2 二維熱傳導(dǎo)反問題的參數(shù)辨識方法

采用共軛梯度法進行接觸熱阻的反辨識。共軛梯度法的求解策略分為:靈敏度算法求解迭代步長、伴隨算法求解迭代方向和共軛系數(shù)、邊界元法求解正問題。其迭代格式為:

(9)

2.1 靈敏度算法求解迭代步長

靈敏度是指當界面溫度Tl(x)有一個增量ΔTl(x)時,區(qū)域內(nèi)某點溫度T(x,y)的變化量ΔT(x,y)如何變化。將T+ΔT,Tl(x)+ΔTl(x)代入(1)式～(5)式中得靈敏度控制方程及相應(yīng)邊界條件為:

(11)式表示兩側(cè)(x=L或x=0處)的熱流變化量邊界條件;(12)式表示上、下表面(y=L、y=0處)的溫度變化量邊界條件。可由上述方程組求解出y=Li處的溫度變化量ΔTi(x)。

將(9)式代入(6)式得:

(14)

(15)

對(15)式中的βk求導(dǎo)使其導(dǎo)數(shù)為0,并將其沿x坐標離散得:

(16)

2.2 伴隨算法求解迭代方向

(17)

(18)

式中，γk為共軛系數(shù),規(guī)定γ0=0。

(19)

式中，λ為Lagrange算子;δ(x-xj)為Dirac函數(shù)。

用Ti(x)+ΔTi(x)代替Ti(x),用Tl(x)+ΔTl(x)代替Tl(x),經(jīng)過一系列變換可得泛函增量為:

(20)

將方程(20)右邊第二項分部積分,并利用靈敏度問題的邊界條件,同時令ΔJ→0,可得伴隨問題:

(22)式表示側(cè)面(x=L或x=0處)的伴隨問題邊界條件;(23)式、(24)式表示上、下表面(y=L,y=0處)的伴隨問題邊界條件。同樣,可由上述方程組求出伴隨函數(shù)λ(x,y)的相應(yīng)值。

根據(jù)Alifanov的定義[14],由泛函增量可以求得:

(25)

2.3 邊界元法求解正問題

正問題是根據(jù)已知的邊界條件求解微分方程,得到未知的熱流或溫度,方程組(1)～(5)、(10)～(13)、(21)～(24)作為偏微分方程組的邊界值問題,其求解過程均可歸結(jié)為正問題求解的范疇。本文采用邊界元方法求解正問題,以方程組(1)～(5)為例,首先將控制方程化為邊界積分形式:

(26)

式中,Γ為邊界;Ti為場點(xi,yi)處的溫度值;Tj為源點處的溫度值;qj為源點處的熱流值;q*=-k?G*/?n,G*為控制方程的基本解,在二維問題中,G*=ln(1/R)/(2πk),R為場點到源點的距離;Ci為自由項系數(shù),當(xi,yi)點位于區(qū)域內(nèi)時Ci=1,當(xi,yi)點位于邊界時Ci=1/2。

如圖2所示,將邊界Γ劃分成N個單元,每個單元記為Γj(j=1,2,…,N)。

圖2 二維熱傳導(dǎo)問題的邊界元模型

采用常單元插值,可將(26)式寫成:

(27)

(28)

將(28)式中未知量移到等式左邊,已知量移到等式的右邊,寫成通常的方程組形式為:

[A][X]=[F]

(29)

采用四點高斯積分公式，得到數(shù)值積分的系數(shù)。式中，[A]為數(shù)值積分系數(shù)Hij、Gij組成的矩陣；[X]為未知量組成的向量；[F]為已知量組成的向量。結(jié)合已知邊界條件,按(29)式即可求解出邊界未知溫度和熱流。

如果要求區(qū)域內(nèi)的溫度,可以將(27)式寫成:

(30)

結(jié)合已知的邊界處的熱流值和溫度值,由(30)式即可得到區(qū)域內(nèi)部的溫度值。

2.4 求解接觸熱阻

(31)

具體計算流程如圖3所示。

圖3 二維接觸熱阻計算流程圖

3 算例分析

3.1 算法驗證及反問題不適定性分析

如圖4所示:2塊1 m×1 m的鋼板在y=0處接觸;側(cè)面(x=0 m或x=1 m處)為絕熱;上、下表面恒溫且Tu=500 K、Td=0 K;2個區(qū)域材料的導(dǎo)熱系數(shù)為k1=k2=14.9 W/(m·K);在y=0.05 m處有一系列溫度測點,設(shè)Ti(x)=282.87+0.25sin(6πx)(單位為K),其中0.25sin(6πx)為模擬的測量誤差,可見測量誤差不超過0.09%;假設(shè)界面接觸熱阻為:R(x)=0.001 5-0.004x+0.004x2+0.000 05·

sin(4πx)(單位為K·m2/W),稱為準確接觸熱阻。

在不考慮溫度測點誤差的情況下,利用本文第2節(jié)的參數(shù)辨識方法可以計算出上界面溫度,再通過界面接觸熱阻反推出下界面溫度;同理,在考慮溫度測點誤差的情況下,也可以計算出上、下界面溫度。利用溫度測點無誤差情況下計算出的上界面溫度和溫度測點有誤差情況下計算出的下界面溫度可以計算出考慮溫度測點誤差的接觸熱阻,稱為計算接觸熱阻,與準確接觸熱阻的對比如圖5所示。

圖4 計算模型

圖5 準確接觸熱阻和計算接觸熱阻的對比圖

在離散化的準確接觸熱阻和計算接觸熱阻中各取n個值,求取平均相對誤差,計算公式為:

(32)

式中,R*(xi)為計算接觸熱阻離散值;R(xi)為準確接觸熱阻離散值。

再由Bessel公式可得計算接觸熱阻的方均根誤差計算公式為:

(33)

計算可得:v=6.15%,σ=5.31×10-4Km2/W。該算例表明:本文提出的邊界元法和共軛梯度法相結(jié)合的參數(shù)辨識方法可以解決接觸熱阻反問題;但反問題的求解結(jié)果對溫度測點處的誤差比較敏感,很小的輸入誤差(不到0.09%)會導(dǎo)致較大的輸出誤差(6.15%),并且輸出結(jié)果的穩(wěn)定性也不高(方均根誤差較大),即反問題的求解存在不適定性。

3.2 測點位置對求解精度的影響

結(jié)合上述模型,分別在距離上界面0.05 m,0.10 m,0.15 m,0.20 m,0.25 m,0.30 m處布置溫度測點,并計算接觸熱阻及其相應(yīng)的平均相對誤差,如圖6所示。

圖6 接觸熱阻平均相對誤差隨測點位置y變化的曲線

由圖6結(jié)果可知:測點位置距離接觸界面越遠,接觸熱阻計算誤差越大,并且相對誤差增加幅度也越大,該結(jié)論與文獻[13]中結(jié)論一致。因此,在靠近接觸界面處布置溫度測點并讀取溫度測量值可以提高接觸熱阻計算精度。

3.3 基于最小二乘法的不適定性處理方法

利用區(qū)域Ⅰ中測點溫度的準確值計算出上界面溫度,利用區(qū)域Ⅱ中測點溫度的準確值計算出下界面溫度,結(jié)合計算出的熱流,即可計算出準確接觸熱阻值R(x);同理,可以計算出考慮測量誤差的計算接觸熱阻值R*(x)。下面采用最小二乘法對計算接觸熱阻進行優(yōu)化,得到修正接觸熱阻值Rcr(x):

在不失一般性的情況下,不妨假設(shè):

(34)

式中,αn為未知系數(shù)。定義:

(35)

分別求(35)式對α1,…,αn的偏導(dǎo),并令結(jié)果為0,可得到如下方程組:

(36)

求解(36)式即可解出未知的系數(shù)αn,進而代入(34)式求出Rcr(x)。本文取n=6,求得:

(37)

對比R(x)、R*(x)、Rcr(x),如圖7所示。

圖7 R(x),R*(x),Rcr(x)對比圖

利用(32)式、(33)式計算可得:R*(x),Rcr(x)的平均相對誤差分別為5.06%,2.96%;二者的方均根誤差分別為6.11×10-5Km2/W,5.69×10-5Km2/W。可見經(jīng)過最小二乘法處理后的接觸熱阻的精度得到了提高,并且計算結(jié)果的穩(wěn)定性也有所提高。

4 結(jié) 論

1) 考慮到大多數(shù)研究忽略了接觸熱阻隨接觸面空間的變化,本文采用邊界元法和共軛梯度法相結(jié)合的方法反辨識出二維熱傳導(dǎo)問題中沿接觸線位置變化的接觸熱阻。對于接觸區(qū)域較大或接觸面拓撲形貌變化較大等不宜將模型簡化為一維熱傳導(dǎo)問題的情況,本方法可以計算出更準確的接觸熱阻;

2) 考慮到邊界元算法不需要對內(nèi)部區(qū)域進行離散的特點,本方法可以全域選擇溫度測點;

3) 算例分析表明,本方法可以有效辨識出接觸熱阻,但計算結(jié)果對溫度測量誤差非常敏感,并且溫度測點離接觸界面越遠敏感程度越高;

4) 利用最小二乘法對接觸熱阻的辨識結(jié)果進行優(yōu)化處理,提高了接觸熱阻的辨識精度和穩(wěn)定性;

5) 本方法基于二維兩層介質(zhì)模型建立,具有一定的代表性,可以用于類似的接觸熱阻測試工作,為計算接觸熱阻提供了新的思路;

6) 可以更進一步將本方法推廣至三維接觸熱阻的分析情況,計算出更精確的隨接觸界面位置變化的接觸熱阻。

[1] Madhusudana C V, Fletcher L S. Contact Heat Transfer——The Last Decade[J]. AIAA Journal, 1986, 24(3):510-523

[2] 王安良, 趙劍鋒. 接觸熱阻預(yù)測的研究綜述[J]. 中國科學(xué): 技術(shù)科學(xué), 2011, 41(5):545-556

Wang Anliang, Zhao Jianfeng. The Research Overview of Prediction of Thermal Contact Resistance[J]. Scientia Sinica Techologica, 2011, 41(5):545-556 (in Chinese)

[3] Ramamurthi K, Kumar S S, Abilash P M. Thermal Contact Conductance of Molybdenum-Sulphide-Coated Joints at Low Temperature[J]. Thermophys Heat Transf, 2007, 21(4):811-813

[4] Bahrami M, Yovanovich M M, Culham J R. Thermal Contact Resistance at Low Contact Pressure: Effect of Elastic Deformation[J]. International Journal of Heat and Mass Transfer, 2005, 48(16):3284-3293

[5] 劉冬歡, 王飛, 曾凡文,等. 高溫接觸熱阻的有限元模擬方法[J]. 工程力學(xué), 2012, 29(9):375-379

Liu Donghuan, Wang Fei, Zeng Fanwen, et al. Finite Element Simulation Method of High Temperature Thermal Contact Resistance[J]. Engineering Mechanics, 2012, 29(9):375-379 (in Chinese)

[6] 李磊濤, 鄭曉亞. 一種瞬態(tài)接觸熱導(dǎo)數(shù)值模擬方法[J]. 機械工程學(xué)報, 2016, 52(6):174-180

Li Leitao, Zheng Xiaoya. A Numerical Simulation Method of Transient Thermal Contact Conductance[J]. Journal of Mechanical Engineering, 2016, 52(6):174-180 (in Chinese)

[7] 李建平, 劉濤, 王伯長,等. 鋁卷徑向接觸熱阻的仿真研究[J]. 金屬熱處理, 2009, 34(4):91-95

Li Jianping, Liu Tao, Wang Bochang, et al. Simulation of the Radial Contact Thermal Resistance for the Aluminum Coil[J]. Heat Treatment of Metals, 2009, 34(4):91-95 (in Chinese)

[8] 劉冬歡, 鄭小平, 黃拳章,等. C/C復(fù)合材料與高溫合金GH600之間高溫接觸熱阻的試驗研究[J]. 航空學(xué)報, 2010, 31(11):2189-2194

Liu Donghuan, Zheng Xiaoping, Huang Quanzhang, et al. Experimental Investigation of High-Temperature Thermal Contact Resistance between C/C Composite Material and Superalloy GH600[J]. ACTA Aeronautica et Astronautica Sinica, 2010, 31(11): 2189-2194 (in Chinese)

[9] 石友安, 桂業(yè)偉, 杜雁霞,等. 相變材料熱控系統(tǒng)內(nèi)部接觸熱阻的辨識方法研究[J]. 實驗流體力學(xué), 2012, 26(4):54-58

Shi You′an, Gui Yewei, Du Yanxia, et al. Study on Thermal Contact Resistance Estimation in the Phase-Change Material Thermal Control Device[J]. Journal of Experiments in Fluid Mechanics, 2012, 26(4):54-58 (in Chinese)

[10] El-Sabbagh A M, Fieberg C, El-Marg E, et al. Modeling of Transient Thermal Contact Resistance out of Conjugate Gradient Method[J]. Materialwissenschaft und Werkstofftechnik, 2007, 38(4):288-293

[11] Gill J, Divo E, Kassab A J. Estimating Thermal Contact Resistance Using Sensitivity Analysis and Regularization[J]. Engineering Analysis with Boundary Elements, 2009, 33(1):54-62

[12] 王超. 高強鋼熱成形接觸導(dǎo)熱和零件力學(xué)性能及工藝優(yōu)化研究[D]. 武漢:華中科技大學(xué), 2014

Wang Chao. Investigation on Thermal Contact Behavior and Prediction of Mechanical Properties and Process Optimization in Hot-Stamping[D]. Wuhan, Huazhong University of Science and Technology, 2014 (in Chinese)

[13] 張平, 宣益民, 李強. 界面接觸熱阻的研究進展[J]. 化工學(xué)報, 2012, 63(2):335-349

Zhang Ping, Xuan Yimin, Li Qiang. Development on Thermal Contact Resistance[J]. CIESC Journal, 2012, 63(2):335-349 (in Chinese)

[14] Alifanov O M. Solution of an Inverse Problem of Heat Conduction by Iteration Methods[J]. Journal of Engineering Physics and Thermophysics, 1974, 26(4): 471-476 (in Chinese)