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

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

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

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

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

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

1 正問題

1.1 導熱理論

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

其中，a為熱擴散系數(shù)；λ為導熱系數(shù)。

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

其中

1.2 邊界積分方程的離散

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

對時間內(nèi)層求積分得：

其中

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

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

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

整理（9）式可得：

其中

Pi為整個空間Ω上的初始溫度對源點i的影響。將（12）式寫成矩陣形式為：

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

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

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

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

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

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

其中

其中，Nj（ξl，ηk）為四點線性插值公式；｜J｜為坐標變換的雅戈比行列式。和的計算如下：

2 反問題

2.1 目標函數(shù)

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

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

2.2 復變量求導法

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

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

當h取極小值時，通?？梢院雎匀A以上的無窮小量，分別比較實部和虛部可得：

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

2.3 共軛梯度法

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

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

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

應(yīng)用復變量求導法的（28）式計算目標函數(shù)的梯度：

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

其中，▽fk為行向量，可由復變量求導法獲得。

2.4 反演問題的求解流程

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

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

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

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

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

3 結(jié)果與討論

3.1 算例

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

圖1 單元和節(jié)點

3.2 雙參數(shù)反演

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

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

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

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

圖4 導熱系數(shù)λ的計算結(jié)果

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

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

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

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

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

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

4 結(jié)束語

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

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

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

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

［4］ 程榮軍，程玉民.帶源參數(shù)的熱傳導反問題的無網(wǎng)格方法［J］.物理學報，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］ 周煥林，江 偉，胡 豪，等.二維彈性力學邊界條件反識別TSVD正則化方法［J］.合肥工業(yè)大學學報：自然科學版，2013，36（9）：1076－1081.

［9］ 卞歩喜，周煥林，程長征，等.二維位勢邊界條件反識別TSVD正則化法［J］.合肥工業(yè)大學學報：自然科學版，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.