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

（中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川 綿陽(yáng) 621000）

0 引 言

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

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

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

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

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

1 計(jì)算模型

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

將傳熱簡(jiǎn)化為一維問(wèn)題，則控制方程如下：

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

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

考慮接觸熱阻時(shí)：

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

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

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

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

2 參數(shù)辨識(shí)方法

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

式中符號(hào)意義：測(cè)點(diǎn)的坐標(biāo)xm；測(cè)量溫度Tm；計(jì)算溫度T；測(cè)量誤差ε。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步長(zhǎng)由下式計(jì)算：

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

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

σ為測(cè)量誤差的均方差，M為測(cè)點(diǎn)數(shù)。

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

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

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

圖4 實(shí)驗(yàn)系統(tǒng)圖Fig.4 Scheme of experiment system

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

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

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

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

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

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

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

圖9 測(cè)點(diǎn)與邊界處的靈敏度的比較Fig.9 Comparison of sensitivity between points in boundary and points in interior

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

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

4 小 結(jié)

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

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

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

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

［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］ 錢(qián)煒祺，蔡金獅.再入航天飛機(jī)表面熱流密度辨識(shí)［J］.宇航學(xué)報(bào)，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］ 石友安，曾磊，錢(qián)煒祺，等.熱參數(shù)辨識(shí)在測(cè)熱試驗(yàn)中的應(yīng)用研究［J］.工程熱物理學(xué)報(bào)，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.