相變材料熱控系統(tǒng)內(nèi)部接觸熱阻的辨識方法研究

石友安，桂業(yè)偉，杜雁霞，曾 磊，錢煒祺

（中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

0 引 言

兩種相互接觸的固體材料的交界面，由于表面粗糙度的存在，微觀上會呈現(xiàn)出“非一致接觸”。當(dāng)熱流從一種材料經(jīng)過界面向另一種傳遞時，由于“非一致接觸”的存在，傳熱方式將發(fā)生變化：導(dǎo)熱部分的熱流將向接觸點處收縮，空穴內(nèi)的氣體會帶來對流傳熱，高溫情況下，非接觸的壁面還會發(fā)生輻射換熱。這些變化都將對熱量的傳遞形成一種阻礙作用，表現(xiàn)為界面附近的溫度分布明顯的降低，如圖1所示。表征界面因非一致接觸而產(chǎn)生的阻礙熱流傳遞作用的物理量稱為接觸熱阻［1］，數(shù)學(xué)定義式為：Rc＝ΔT／q。

相變材料熱控系統(tǒng)是航天器及空間飛行器常用的一種熱控系統(tǒng)。在封裝過程中，由于材料凝固時的體積變化將導(dǎo)致相變材料與容器壁面之間形成非一致接觸，從而產(chǎn)生接觸熱阻。對于高溫條件下的熱控系統(tǒng)而言，接觸熱阻的存在，在一定程度上降低了傳熱效率，容易產(chǎn)生局部高溫，引發(fā)局部熱應(yīng)力，甚至帶來局部破壞［2］。接觸熱阻已逐漸成為高溫條件下熱量精確控制的一項瓶頸因素。圖2顯示了在不考慮接觸熱阻情況下，以正二十八烷為相變材料的熱控裝置的計算與實驗吸熱融化過程相界面運動速度的比較。

圖1 微觀視角下的界面?zhèn)鳠岱绞紽ig.1 Heat transfer in interface between solids in microscopic view

由于接觸界面處傳熱機理復(fù)雜、影響因素眾多，譬如接觸壓力、界面溫度等，且相變材料與容器之間的接觸熱阻又與材料凝固過程晶體的形成和生長特性密切相關(guān)，現(xiàn)有的理論計算方法大多過于簡化，無法直接應(yīng)用于工程實際。因此，采用反演熱傳導(dǎo)反問題的方法，將參數(shù)辨識中的靈敏度法和共軛梯度法應(yīng)用到實驗，建立了一種可用于穩(wěn)態(tài)和瞬態(tài)接觸熱阻評估的辨識方法。通過模擬辨識和初步應(yīng)用，得到了一些有意義的結(jié)論。

圖2 在不考慮接觸情況下，正二十八烷的計算與實驗相界面移動特性比較Fig.2 Comparison of phase interface of C28H58between computation and experiment without consideration of thermal contact resistance（Rc）

1 計算模型

考慮接觸熱阻時材料融化的計算模型如圖3所示。此時，材料內(nèi)部的傳熱過程分為熔化前的導(dǎo)熱過程和熔化開始后的固／液相變過程。

將傳熱簡化為一維問題，則控制方程如下：

不考慮接觸熱阻時，邊界條件為第一類邊界條件：

初值條件：t＝1，T＝T0；

考慮接觸熱阻時：

此時，固體內(nèi)部的傳熱仍然滿足一維傳熱方程，即式（1），但是，邊界條件不同：

初值條件：t＝0，T＝T0；

綜上所述，當(dāng)存在接觸熱阻時，內(nèi)部傳熱機制不變，但是邊界條件等效為從第一類轉(zhuǎn)變?yōu)榈谌悺?/p>

圖3 熱控裝置內(nèi)接觸熱阻的計算模型Fig.3 Calculation model of Rcin heat control equipment

2 參數(shù)辨識方法

辨識求解的思路是：以接觸點測量溫度Tm的時間歷程為依據(jù)，在數(shù)值求解熱傳導(dǎo)正問題的基礎(chǔ)上，根據(jù)輸出誤差原則，將反問題轉(zhuǎn)化為一個優(yōu)化問題來處理，等價于在函數(shù)空間中尋求合適的函數(shù)Rc或Rc（t）或Rc（xm，t），使如下目標(biāo)泛函達到極小的過程：

式中符號意義：測點的坐標(biāo)xm；測量溫度Tm；計算溫度T；測量誤差ε。

2.1 穩(wěn)態(tài)接觸熱阻的辨識

采用靈敏度法進行穩(wěn)態(tài)接觸熱阻［3］的辨識，其優(yōu)化策略是牛頓 拉夫遜（Newton－Raphson）算法［4］。下面介紹算法：

首先根據(jù)輸出誤差原則建立目標(biāo)泛函：

優(yōu)化算法描述如下：

式中，參數(shù)的上標(biāo)k，k＋1分別表示迭代前、后的參數(shù)值，rex是松弛因子。式中?T／?Rc為靈敏度，可以通過求解靈敏度方程得出。將式（1）對參數(shù)求導(dǎo)后即得到靈敏度方程。具體形式為（U定義為?T／?Rc）：

初值條件：t＝0，U＝0；

至此，穩(wěn)態(tài)接觸熱阻辨識基本建立，依據(jù)式（5）進行辨識即可得到Rc。

2.2 瞬態(tài)接觸熱阻的辨識

相較于穩(wěn)態(tài)而言，瞬態(tài)接觸熱阻［5］更能細(xì)致地刻畫接觸熱阻的變化關(guān)系。采用共軛梯度法［6－7］進行辨識。算法描述：

式中下標(biāo)i表示接觸熱阻的時間離散次數(shù)；上標(biāo)l、l＋1表示迭代層次，β為步長，P為共軛梯度，具體優(yōu)化過程表述如下：

式中：J′l＝ （?J／?Rc，i）l

J′l－1和J′l是第l－1步和l步迭代的梯度下降方向，可以采用靈敏度分析或求解伴隨方程來獲取。由于滿足目標(biāo)泛函的溫度場T，也必須滿足主控方程，因此可根據(jù)拉格朗日乘數(shù)原理對目標(biāo)函數(shù)進行擴展，得到如下擴展泛函：

式中λ為伴隨變量。對上式后半部分積分再做變分，得到伴隨方程：

初值條件：t＝tmax，λ（x，tmax）＝0

邊值條件：λ（0，t）｜x＝0＝0；

式中δ為克羅內(nèi)克符號。整理得到，目標(biāo)函數(shù)對Rc（T）的梯度為：

步長由下式計算：

式中ΔT是由ΔR＝λn引起的溫度場的變化值，可通過解靈敏度方程獲得。

在優(yōu)化過程中，收斂準(zhǔn)則根據(jù)輸出誤差原則獲得，即：

σ為測量誤差的均方差，M為測點數(shù)。

至此，共軛梯度法已完全建立，依據(jù)式（7）進行辨識即可得到Rc（t）。

3 辨識方法的初步應(yīng)用

采用封裝二十八烷長方體熱控單元進行融化實驗，實驗系統(tǒng)見圖4，測點位置為x＝0.006m，y＝0.040m。采用電加熱器作為恒溫?zé)嵩?，溫度控制采用ANV TF100PID自動演算控制器，通過加熱功率的調(diào)節(jié)實現(xiàn)實驗過程所需的加熱邊界條件，溫度控制精度為±1K。溫度信號由Agilent 34970A多點數(shù)據(jù)采集器采集。

圖4 實驗系統(tǒng)圖Fig.4 Scheme of experiment system

實驗測量溫度數(shù)據(jù)見圖5，從圖中可以看出，溫升斜率的增長趨勢是先增加后減小，由此可以判定加熱邊界在測量記錄中已經(jīng)出現(xiàn)了融化。因此，辨識數(shù)據(jù)選取在0～300s內(nèi)邊界未融化的一段時間。

圖5 測點溫度數(shù)據(jù)Fig.5 Data of temperature at measured point

根據(jù)測點的溫升歷程測量參數(shù)進行辨識。分別取測量時間為100、140s進行穩(wěn)態(tài)接觸熱阻的辨識，得到的結(jié)果依次為0.0799和0.0665。圖6是利用穩(wěn)態(tài)辨識值計算得到的測點溫度與實測值的對比?？梢园l(fā)現(xiàn)：采用穩(wěn)態(tài)辨識結(jié)果的計算溫度值與實測值存在一定的差異，且隨著時間的增長，差異增大，分析認(rèn)為隨著相變材料溫度的升高，接觸熱阻可能隨時間發(fā)生了變化，穩(wěn)態(tài)接觸熱阻的假設(shè)可能帶來了原理性誤差。

圖6 測量溫度和計算溫度的對比（穩(wěn)態(tài)結(jié)果）Fig.6 Comparison of temperature between calculation and experiment in steady Rc

為了提高界面接觸熱阻的辨識精度，進一步開展了瞬態(tài)接觸的辨識。圖7是瞬態(tài)接觸熱阻的辨識值。圖8是測點的實測溫度和計算溫度的對比?？梢园l(fā)現(xiàn)：利用瞬態(tài)辨識值得到的測點計算溫度能夠更好的吻合實測溫度。分析表明，隨著加熱時間的延長，材料溫度不斷升高，壁面和材料之間的空穴對熱流的阻礙作用將減小，當(dāng)材料出現(xiàn)融化后，融化的液體將填充空穴，接觸熱阻將繼續(xù)減小，直至湮滅。這進一步證明了進行瞬態(tài)接觸熱阻辨識的可行性和可信性。

圖8 測量溫度和計算溫度的對比（瞬態(tài)接觸熱阻）Fig.8 Comparison of temperature between calculation and experiment in unsteady Re

另外，對影響辨識質(zhì)量的因素做了初步分析，發(fā)現(xiàn)：測溫數(shù)據(jù)的信噪比是一個影響辨識質(zhì)量的重要因素。測溫數(shù)據(jù)是反問題求解的基礎(chǔ)，其信噪比的高低直接決定求解過程能否克服非穩(wěn)態(tài)傳熱固有的擴散性，影響解的穩(wěn)定性，可參閱文獻［6－8］。保障測量數(shù)據(jù)具有較高信噪比的措施之一是合理地安排測點的位置。圖8是不同測點處的計算靈敏度的對比。從圖中可以得到：越靠近加熱邊界，靈敏度越高。即在實驗允許的情況下，測點位置宜盡可能地靠近擾動邊界，以獲得較高信噪比的測溫數(shù)據(jù)。

圖9 測點與邊界處的靈敏度的比較Fig.9 Comparison of sensitivity between points in boundary and points in interior

為進一步研究接觸熱阻對材料熔化的延遲作用，初步應(yīng)用穩(wěn)態(tài)接觸熱阻辨識值，對75℃等溫加熱條件下，考慮和不考慮接觸熱阻效應(yīng)的材料融化計算速率與相同條件下實驗熔化速率進行了比較，如圖9所示。當(dāng)計及穩(wěn)態(tài)接觸熱阻時，實驗與計算熔化速率之間的吻合度有所提高，如圖10所示。特別是融化延遲段很好地吻合試驗，說明考慮界面接觸熱阻有利于提高熔化時間的預(yù)測準(zhǔn)度。

圖10 穩(wěn)態(tài)接觸熱阻對熔化速率的影響Fig.10 Influence of melt rate with consideration of steady Rc

4 小 結(jié)

將參數(shù)辨識方法引入相變過程界面接觸熱阻的辨識中，將數(shù)值計算與實驗相結(jié)合，通過參數(shù)辨識獲得了界面接觸熱阻的定量參數(shù)，對有效評估界面熱阻對相變傳熱的影響作用提供了一種新的思路，得到了一些有意義的結(jié)果：

（1）所建立的穩(wěn)態(tài)和瞬態(tài)接觸熱阻的兩種辨識算法具有較高的精度、穩(wěn)定性好、抗噪性強，算法有效；

（2）由于界面熱阻的獲得具有較大難度，采用基于傳熱反問題的參數(shù)辨識方法，是獲得相變材料與容器之間界面熱阻的一種有效途徑。該方法對于獲取復(fù)合相變材料等不易直接測量的物性參數(shù)也具有一定的參考價值；

（3）雖然初步應(yīng)用驗證了辨識算法的可行性，但是，由于接觸熱阻值較小，影響辨識質(zhì)量的因素較多，且相變傳熱實驗中還存在一些不確定因素，例如物性隨溫度變化。因此將辨識算法應(yīng)用于實際工程，還需開展進一步的深入研究。

［1］ MADHUSUDANA C V，F(xiàn)LETCHER L S.Contact heat transfer－the last decade［J］.AIAA Journal，1986，24（3）：510－523.

［2］ WEI D，SUN X S，ZHU J M，et al.Structural response and behavior analysis of steel joint components at elevated temperatures［M］.Key Engineering Materials.2007，353－358：2672－2675.

［3］ JENNIFER Gill，EDUARDO Divo，ALAIN J Kassab.Estimating thermal contact resistance using sensitivity analysis and regularization［J］.Engineering Analysis with Boundary Elements，2009，33：54－62.

［4］ 錢煒祺，蔡金獅.再入航天飛機表面熱流密度辨識［J］.宇航學(xué)報，2000，21（4）：1－6.

［5］ FIEBERG C，KNEER R.Determination of thermal contact resistance from transient temperature measurements［J］.International Journal of Heat and Mass Transfer，2008，51：1017－1023.

［6］ 石友安，曾磊，錢煒祺，等.熱參數(shù)辨識在測熱試驗中的應(yīng)用研究［J］.工程熱物理學(xué)報，2010，31（9）：1555－1558.

［7］ OLIVEIRAY A P D，ORLANDE H R B.Estimation of the heat flux at the surface of ablating materials by using temperature and surface position measurements［J］.Inverse Problems in Science and Engineering，2004，20（5）：563－577.

［8］ BECK J V，BLACKWELL B，CLAIR C R S.Inverse heat conduction ill－posed problems［M］.New York：John Wiley ＆Sons，1985.