二維瞬態(tài)熱傳導(dǎo)問題熱物性參數(shù)反演

周煥林， 徐興盛， 李秀麗， 胡 豪

（合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，安徽 合肥 230009）

在航空航天、冶金鑄造、化工制藥、材料冶金、機(jī)械制造、交通運(yùn)輸、核反應(yīng)堆、地?zé)崮芸碧?、生物傳熱、熱工測量、土木工程、無損探傷、冷凍儲(chǔ)藏等工程領(lǐng)域，熱傳導(dǎo)反問題是傳熱學(xué)研究的熱點(diǎn)之一，它是利用實(shí)驗(yàn)手段測得物體內(nèi)部或邊界上某些點(diǎn)的溫度、熱流及其隨時(shí)間的變化歷程，通過求解導(dǎo)熱微分方程來反演物體熱通量、材料熱傳導(dǎo)系數(shù)或物體內(nèi)部熱源分布等參數(shù)。

隨著反演問題理論研究的不斷深入，研究者提出了多種求解熱傳導(dǎo)反問題的方法。文獻(xiàn)［1］將材料的熱傳導(dǎo)系數(shù)值按溫度區(qū)間分段離散，基于遺傳算法和伴隨方程法通過材料邊界點(diǎn)的溫度測量來反演各溫度區(qū)間熱傳導(dǎo)系數(shù)值；文獻(xiàn)［2］建立了非線性穩(wěn)態(tài)熱傳導(dǎo)問題的有限元模型，對(duì)非線性熱物性和邊界條件進(jìn)行了反演；文獻(xiàn)［3］研究了非線性瞬態(tài)熱傳導(dǎo)問題中導(dǎo)熱系數(shù)和邊界條件的多宗量反演；文獻(xiàn)［4］基于無網(wǎng)格有限點(diǎn)法反演一維熱傳導(dǎo)問題的源參數(shù)；文獻(xiàn)［5］利用直接積分法研究了一維非線性熱傳導(dǎo)反問題，給出了熱傳導(dǎo)系數(shù)的反演結(jié)果；文獻(xiàn)［6－7］利用有限差分法與高斯－賽德爾迭代法，研究了一維瞬態(tài)熱傳導(dǎo)問題的熱參數(shù)；文獻(xiàn)［8］采用截?cái)嗥娈愔嫡齽t化方法反演二維各向同性彈性力學(xué)Cauchy問題的邊界條件；文獻(xiàn)［9］采用截?cái)嗥娈愔捣ǚ囱荻S各向同性材料Cauchy位勢問題的未知邊界條件；文獻(xiàn)［10］利用基本解法與移動(dòng)最小二乘法對(duì)未知邊界上的溫度值進(jìn)行反演計(jì)算；文獻(xiàn)［11］采用數(shù)值積分方法對(duì)熱通量進(jìn)行了反演識(shí)別，計(jì)算過程中使用了未來時(shí)刻的溫度值；文獻(xiàn)［12－13］利用共軛梯度法反演一維非線性瞬態(tài)熱傳導(dǎo)問題中的導(dǎo)熱系數(shù)和熱容量，且進(jìn)一步反演了二維非均質(zhì)材料的熱傳導(dǎo)系數(shù)；文獻(xiàn)［14］使用復(fù)變量求導(dǎo)法計(jì)算靈敏度矩陣。

近些年來，非線性反演方法迅猛發(fā)展，已經(jīng)成為反演理論方法的重中之重。非線性反演方法主要有：梯度法、牛頓法、共軛梯度法、蒙特卡洛法、變尺度法、模擬退火法、遺傳算法和人工神經(jīng)網(wǎng)格法等。共軛梯度法是介于梯度法和牛頓法之間的一種有效處理反問題的方法，并且集成了2種方法的優(yōu)點(diǎn)。共軛梯度法僅僅利用了目標(biāo)函數(shù)的一階偏導(dǎo)數(shù)的信息，既克服了梯度法收斂速度慢的缺點(diǎn)，又省去了牛頓法在每一次搜索時(shí)需要重復(fù)計(jì)算和儲(chǔ)存Hessian矩陣并求其逆矩陣的麻煩。并且共軛梯度法還具有收斂性和穩(wěn)定性高、所占存儲(chǔ)空間小、不需要引入外來參數(shù)以及具有如牛頓法那樣在極小值點(diǎn)附近收斂速度快的優(yōu)點(diǎn)。

本文基于邊界元法反演二維瞬態(tài)導(dǎo)熱問題的熱擴(kuò)散系數(shù)和導(dǎo)熱系數(shù)，并討論迭代初值、數(shù)據(jù)隨機(jī)偏差等對(duì)反演結(jié)果的影響。

1 正問題

1.1 導(dǎo)熱理論

對(duì)于二維非穩(wěn)態(tài)熱傳導(dǎo)問題，滿足方程：

其中，a為熱擴(kuò)散系數(shù)；λ為導(dǎo)熱系數(shù)。

由加權(quán)余量法得到二維瞬態(tài)導(dǎo)熱問題的邊界積分方程：

其中

1.2 邊界積分方程的離散

首先要對(duì)時(shí)間域離散，假設(shè)函數(shù)T、q隨時(shí)間變化，由于T和q比T＊和q＊的變化要慢得多，由此近似認(rèn)為在極短的時(shí)間段內(nèi)為常數(shù)，則（2）式轉(zhuǎn)變?yōu)椋?/p>

對(duì)時(shí)間內(nèi)層求積分得：

其中

其中，C為歐拉常數(shù)，C＝0.577 215 66。

將邊界Γ劃分成N個(gè)線性單元，空間域Ω劃分成M個(gè)四邊形單元，則

單元內(nèi)任意一點(diǎn)的溫度或熱通量可以由單元端點(diǎn)的值通過線性插值來確定。插值函數(shù)為：

整理（9）式可得：

其中

Pi為整個(gè)空間Ω上的初始溫度對(duì)源點(diǎn)i的影響。將（12）式寫成矩陣形式為：

其中，H為溫度系數(shù)矩陣；G為熱通量系數(shù)矩陣；Tt2為節(jié)點(diǎn)溫度；Qt2為節(jié)點(diǎn)熱通量。

引入邊界條件后，可將（16）式轉(zhuǎn)化為：

向量X包含邊界上的未知量。

1.3 系數(shù)矩陣的計(jì)算

（1）非對(duì)角元素的計(jì)算。如果源點(diǎn)不位于積分單元內(nèi)，（13）式和（14）式可采用四點(diǎn)高斯求積公式，即

（2）對(duì)角元素的計(jì)算。hii的計(jì)算，采用均溫場的概念，可避免Ci的計(jì)算，則

其中

其中，Nj（ξl，ηk）為四點(diǎn)線性插值公式；｜J｜為坐標(biāo)變換的雅戈比行列式。和的計(jì)算如下：

2 反問題

2.1 目標(biāo)函數(shù)

在反問題中，熱擴(kuò)散系數(shù)a或?qū)嵯禂?shù)λ是未知的，其他條件與正問題相同。本文反問題中需要的額外條件是熱通量。反演的優(yōu)化目標(biāo)函數(shù)可以寫為（26）式：

其中，x為待反演參數(shù)向量；L為熱通量測點(diǎn)的個(gè)數(shù)（x）為測點(diǎn)熱通量的計(jì)算值；（x）為測點(diǎn)熱通量的實(shí)際值。

2.2 復(fù)變量求導(dǎo)法

求導(dǎo)是反演計(jì)算的重要組成部分，求導(dǎo)的精度極大地影響反演的求解精度。通常求導(dǎo)采用有限差分法，對(duì)于復(fù)雜的系統(tǒng)或函數(shù)，差分法很難滿足計(jì)算精度的要求。復(fù)變量求導(dǎo)法最早由文獻(xiàn)［15］提出，它把偏導(dǎo)數(shù)的計(jì)算轉(zhuǎn)化為復(fù)域函數(shù)的計(jì)算，是一種使用方便、計(jì)算結(jié)果精確的函數(shù)偏導(dǎo)數(shù)數(shù)值計(jì)算方法。

對(duì)于任意一個(gè)實(shí)函數(shù)f（x），將所求導(dǎo)數(shù)的變量x施加一個(gè)很小的虛步h，并將其展開成泰勒級(jí)數(shù)形式：

當(dāng)h取極小值時(shí)，通常可以忽略三階以上的無窮小量，分別比較實(shí)部和虛部可得：

使用復(fù)變量求導(dǎo)法求導(dǎo)的精度遠(yuǎn)遠(yuǎn)大于普通的差分法，當(dāng)h取10－30時(shí)仍然適用，因此，復(fù)變量求導(dǎo)法是求解靈敏度的一個(gè)有效方法，在復(fù)雜函數(shù)的數(shù)值求導(dǎo)計(jì)算中非常有效。

2.3 共軛梯度法

共軛梯度法的基本思路是把共軛性與最速下降法相結(jié)合，利用已知點(diǎn)處的梯度構(gòu)造一組共軛方向，并沿這組方向進(jìn)行搜索，求出目標(biāo)函數(shù)的極小點(diǎn)。待反演參數(shù)x在第n＋1次迭代時(shí)的猜測值為：

其中，αn為搜索步長；dn為共軛搜索方向，是目標(biāo)函數(shù)的梯度方向與前一次搜索方向的線性組合。dn的計(jì)算公式為：

其中，βn為共軛系數(shù)，其計(jì)算公式為：

應(yīng)用復(fù)變量求導(dǎo)法的（28）式計(jì)算目標(biāo)函數(shù)的梯度：

搜索步長αn可以通過優(yōu)化目標(biāo)函數(shù)J（xn－αndn）獲得：

其中，▽fk為行向量，可由復(fù)變量求導(dǎo)法獲得。

2.4 反演問題的求解流程

（1）給定某一小正數(shù)ε作為收斂精度，選擇待反演參數(shù)的初始猜想值x0，令n＝0。

（2）求解 （17）式得到測點(diǎn)處的熱通量，并判斷是否滿足‖dn‖＜ε，如果滿足則停止迭代，否則繼續(xù)。

（3）分別按（31）式和（32）式計(jì)算共軛系數(shù)βn及梯度▽J（xn）。

（4）分別按（30）式及（33）式計(jì)算搜索方向dn及搜索步長αn。

（5）令n＝n＋1，按（29）式對(duì)待反演參數(shù)進(jìn)行更新，返回步驟（2）。

3 結(jié)果與討論

3.1 算例

考慮一個(gè)二維瞬態(tài)導(dǎo)熱問題，區(qū)域Ω是由x＝0，x＝10，y＝0，y＝20構(gòu)成的一個(gè)二維空間。在空間內(nèi)，x＝0和x＝10的兩邊是絕熱的，空間內(nèi)的初始溫度為0，邊界y＝0和y＝20的溫度突然變?yōu)?0。將邊界劃分為30個(gè)線性單元，域內(nèi)劃分為50個(gè)四邊形單元，如圖1所示。熱擴(kuò)散系數(shù)和導(dǎo)熱系數(shù)的精確解a＝5.0，λ＝10.0。在已知測點(diǎn)熱通量的條件下，反演a和λ。

圖1 單元和節(jié)點(diǎn)

3.2 雙參數(shù)反演

在本節(jié)反問題中，材料的熱擴(kuò)散系數(shù)a和導(dǎo)熱系數(shù)λ都是未知的，其他的條件不變。令迭代收斂精度ε為10－4，分別取a和λ的初始參考值為（2，4），選取圖1中節(jié)點(diǎn)2～5和節(jié)點(diǎn)17～20共8個(gè)測點(diǎn)，用測點(diǎn)在t＝2s時(shí)的熱通量反演熱擴(kuò)散系數(shù)a和導(dǎo)熱系數(shù)λ。該算例的熱擴(kuò)散系數(shù)a、導(dǎo)熱系數(shù)λ及目標(biāo)函數(shù)J（x）的收斂曲線，如圖2所示，開始時(shí)反演計(jì)算的收斂速度較快，之后收斂趨于平緩，在經(jīng)過多次迭代后，導(dǎo)熱系數(shù)λ和熱擴(kuò)散系數(shù)a分別收斂于精確解。本文的方法對(duì)熱擴(kuò)散系數(shù)a和導(dǎo)熱系數(shù)λ的同時(shí)反演是有效的。

在雙參數(shù)反演中，選取幾組不同的初始參考值，熱擴(kuò)散系數(shù)a的反演計(jì)算結(jié)果如圖3所示，導(dǎo)熱系數(shù)λ的反演計(jì)算結(jié)果如圖4所示。由圖可見，選取不同的初始值，最后都能收斂到精確解。

圖2 反演參數(shù)及目標(biāo)函數(shù)

圖3 熱擴(kuò)散系數(shù)a的計(jì)算結(jié)果

圖4 導(dǎo)熱系數(shù)λ的計(jì)算結(jié)果

3.3 隨機(jī)偏差對(duì)雙參數(shù)反演結(jié)果的影響

在雙參數(shù)反演的過程中，為了研究測量數(shù)據(jù)的隨機(jī)偏差對(duì)反演結(jié)果的影響，選取迭代收斂精度ε為10－4，分別取熱擴(kuò)散系數(shù)a和導(dǎo)熱系數(shù)λ的初始迭代值為（2，4），選取圖1中節(jié)點(diǎn)2～5和節(jié)點(diǎn)17～20共8個(gè)測點(diǎn)，分別對(duì)測點(diǎn)熱通量施加2%和4%的隨機(jī)偏差，0%是無隨機(jī)偏差的情況。熱擴(kuò)散系數(shù)a和導(dǎo)熱系數(shù)λ各迭代步的計(jì)算結(jié)果如圖5、圖6所示。從圖中可見，在施加一定范圍的隨機(jī)偏差后，本文反演方法仍具有很好的收斂性，隨機(jī)偏差越小，計(jì)算結(jié)果越趨近于精確解。

圖5 隨機(jī)偏差影響下熱擴(kuò)散系數(shù)a的計(jì)算結(jié)果

圖6 隨機(jī)偏差影響下導(dǎo)熱系數(shù)λ的計(jì)算結(jié)果

沒有偏差和施加2%、4%的隨機(jī)偏差得到的反演結(jié)果、其相對(duì)誤差及終止迭代步見表1所列。熱擴(kuò)散系數(shù)a計(jì)算結(jié)果的相對(duì)誤差分別0.03%、0.72%、4.48%，導(dǎo)熱系數(shù)λ計(jì)算結(jié)果的相對(duì)誤差分別為0.04%、0.78%、5.68%，迭代終止步數(shù)分別為31、54、81。由此可見，施加的隨機(jī)偏差越大，產(chǎn)生的相對(duì)誤差也越大，迭代次數(shù)也越多。因此，為保證反演的正確性，隨機(jī)偏差要控制在一定范圍內(nèi)。

表1 隨機(jī)偏差對(duì)反演結(jié)果的影響

4 結(jié)束語

本文基于邊界元法反演了二維瞬態(tài)熱傳導(dǎo)問題的熱擴(kuò)散系數(shù)和導(dǎo)熱系數(shù)。邊界元法用于構(gòu)建二維瞬態(tài)熱傳導(dǎo)問題的數(shù)值分析模型，引入復(fù)變量求導(dǎo)法求解目標(biāo)函數(shù)的梯度矩陣，共軛梯度法用于優(yōu)化目標(biāo)函數(shù)獲得反演結(jié)果。算例表明二維瞬態(tài)熱傳導(dǎo)問題的熱擴(kuò)散系數(shù)和導(dǎo)熱系數(shù)的反演方法有效，具有較高的精度和較好的穩(wěn)定性。另外選取不同的初始迭代值都能得到精確的反演結(jié)果。測量數(shù)據(jù)的隨機(jī)偏差越小，計(jì)算結(jié)果越精確。

［1］ 唐中華，錢國紅，錢煒祺.材料熱傳導(dǎo)系數(shù)隨溫度變化函數(shù)的反演方法［J］.計(jì)算力學(xué)學(xué)報(bào)，2011，28（3）：377－382.

［2］ 楊海天，薛齊文.兩級(jí)敏度分析求解非線性穩(wěn)態(tài)多宗量熱傳導(dǎo)反問題［J］.工程熱物理學(xué)報(bào)，2003，24（3）：463－465.

［3］ 薛齊文，魏 偉，楊海天.多宗量瞬態(tài)熱傳導(dǎo)反演識(shí)別［J］.固體力學(xué)學(xué)報(bào)，2009，30（1）：65－69.

［4］ 程榮軍，程玉民.帶源參數(shù)的熱傳導(dǎo)反問題的無網(wǎng)格方法［J］.物理學(xué)報(bào)，2007，56（10）：5569－5574.

［5］ Kim S，Kim M C，Kim K Y.Non－iterative estimation of temperature－dependent thermal conductivity without internal measurements ［J］.International Journal of Heat and Mass Transfer，2003，46（10）：1801－1810.

［6］ Yang C Y.A linear inverse model for the temperature－dependent thermal conductivity determination in one－dimensional problems［J］.Applied Mathematical Modelling，1998，22（1／2）：1－9.

［7］ Yang C Y.Estimation of the temperature－dependent thermal conductivity in inverse heat conduction problems［J］.Applied Mathematical Modelling，1999，23（6）：469－478.

［8］ 周煥林，江 偉，胡 豪，等.二維彈性力學(xué)邊界條件反識(shí)別TSVD正則化方法［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，36（9）：1076－1081.

［9］ 卞歩喜，周煥林，程長征，等.二維位勢邊界條件反識(shí)別TSVD正則化法［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2014，37（9）：1097－1101.

［10］ Dong C F，Sun F Y，Meng B Q.A method of fundamental solutions for inverse heat conduction problems in all anisotropic medium［J］.Engineering Analysis with Boundary Elements，2007，31（1）：75－82.

［11］ Beck J V.Surface Heat flux determination using an integral method［J］.Nuclear Engineering and Design，1968，7（2）：170－178.

［12］ Huang C H，Yan J Y.An inverse problem in simultaneously measuring temperature－dependent thermal conductivity and heat capacity［J］.International Journal of Heat and Mass Transfer，1995，38（18）：3433－3441.

［13］ Huang C H，Chin S C.A two－dimensional inverse problem in imaging the thermal conductivity of a non－h(huán)omogeneous medium［J］.International Journal of Heat and Mass Transfer，2000，43（22）：4061－4071.

［14］ Gao X W，He M C.A new inverse analysis approach for multi－region heat conduction BEM using complex－variabledifferentiation method［J］.Engineering Analysis with Boundary Elements，2005，29（8）：788－795.

［15］ Lyness N，Moler C B.Numerical differentiation of analytic functions［J］.SIAM Journal on Numerical Analysis，1967，4（2）：202－210.