非線性熱傳導(dǎo)逆問題的表面熱流辨識(shí)方法

錢煒祺， 周 宇， 何開鋒， 袁軍婭， 黃建棟

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621000；2.中國(guó)運(yùn)載火箭技術(shù)研究院，北京 100076）

0 引 言

熱傳導(dǎo)逆問題，或稱為熱傳導(dǎo)反問題（IHCP：Inverse Heat Conduction Problem），是利用實(shí)驗(yàn)手段測(cè)得材料內(nèi)部或邊界上某點(diǎn)或某些點(diǎn)上的溫度隨時(shí)間變化歷程，通過求解傳熱微分方程來反演邊界熱流、材料熱傳導(dǎo)系數(shù)或材料內(nèi)部熱源分布等參數(shù)。熱傳導(dǎo)逆問題在航天、核物理、冶金等工業(yè)研究領(lǐng)域中有著廣泛的應(yīng)用背景［1］，是傳熱學(xué)研究的主要領(lǐng)域之一?，F(xiàn)有大多數(shù)熱傳導(dǎo)逆問題研究針對(duì)的都是線性熱傳導(dǎo)方程，即材料的熱物性參數(shù)，包括熱傳導(dǎo)系數(shù)和比熱是常數(shù)的情況。而在實(shí)際工程應(yīng)用中，很多材料的熱物性參數(shù)是隨溫度變化的，此時(shí)的熱傳導(dǎo)方程是非線性的偏微分方程，對(duì)應(yīng)的熱傳導(dǎo)逆問題也稱為非線性熱傳導(dǎo)逆問題［2－3］。非線性熱傳導(dǎo)逆問題主要包括兩方面研究?jī)?nèi)容：一是邊界熱流已知，根據(jù)材料內(nèi)部或邊界溫度場(chǎng)的測(cè)量結(jié)果來辨識(shí)材料熱物性參數(shù)隨溫度的變化規(guī)律；二是材料熱物性參數(shù)隨溫度的變化規(guī)律已知，需要根據(jù)溫度的測(cè)量結(jié)果來辨識(shí)表面熱流。對(duì)于第一類問題，其求解方法在文獻(xiàn)［4－5］中有詳細(xì)介紹，在此不做討論；對(duì)于第二類問題，文獻(xiàn)［3］采用順序函數(shù)法來進(jìn)行了表面熱流辨識(shí)，文獻(xiàn)［7］給出了伴隨方程的形式，但未給出其具體推導(dǎo)。相比之下，國(guó)內(nèi)對(duì)這一問題的研究開展得相對(duì)較少，限制了表面熱流辨識(shí)方法在工程實(shí)際中的使用范圍。因此，本文工作將以一維非線性熱傳導(dǎo)逆問題為研究對(duì)象，系統(tǒng)地建立起兩類表面熱流辨識(shí)方法：順序函數(shù)法和共軛梯度法，并結(jié)合算例來對(duì)算法的有效性和魯棒性進(jìn)行分析。

1 非線性熱傳導(dǎo)逆問題的數(shù)學(xué)模型

典型的一維傳熱如圖1示，厚度為L(zhǎng)的材料上部受表面熱流Q（t），下部絕熱，溫度傳感器置于測(cè)點(diǎn)P處，此時(shí)的熱傳導(dǎo)方程為：

初始條件：t＝0：T＝T0。

觀測(cè)方程：

其中，ρ為材料密度，k（T）為材料熱傳導(dǎo)系數(shù)，Cp（T）為材料比熱，都是隨溫度T變化的函數(shù)，xm為測(cè)點(diǎn)位置，v（t）為測(cè)量噪聲。此時(shí)的熱傳導(dǎo)逆問題就是在k（T）和Cp（T）函數(shù)已知、表面熱流Q（t）未知的情況下，由觀測(cè)方程式（2）中的信息來反演辨識(shí)表面熱流Q（t）。

圖1 一維傳熱示意圖Fig.1 Sketch of 1Dheat conduction problem

該逆問題等價(jià)于求合適的Q（t）使如下目標(biāo)函數(shù)達(dá)極小的優(yōu)化問題：

其中，t＝［0，tf］表示溫度測(cè)量的時(shí)間段，T（xm，t，Q）表示表面熱流為Q時(shí)測(cè)點(diǎn)溫度歷程的計(jì)算值；為相應(yīng)的實(shí)測(cè)值。由于Q（t）是時(shí)間的函數(shù)，因此式（3）實(shí)際上是一個(gè)泛函表達(dá)式，考慮到數(shù)值求解式（1）時(shí)對(duì)時(shí)間的離散，此時(shí)的待優(yōu)化變量為各時(shí)刻的熱流值Qi，即

其中Nt為數(shù)值計(jì)算中時(shí)間方向的離散層數(shù)。則式（3）的優(yōu)化變成了一參數(shù)優(yōu)化問題，下面介紹兩種優(yōu)化策略：順序函數(shù)法（SFM：Sequential Function specification Method）和共軛梯度法（CGM：Conjugate Gradient Method）。

2 順序函數(shù)法

針對(duì)式（1）描述的表面熱流辨識(shí)問題，當(dāng)用順序函數(shù)法來辨識(shí)未知熱流Q（t）時(shí)，問題可描述為：已知tM時(shí)刻的前M－1個(gè)時(shí)刻處的熱流1，…，－1和后r個(gè)時(shí)刻的溫度值TM，TM＋1，…，TM＋r－1，要求辨識(shí)tM時(shí)刻的熱流QM。例如，當(dāng)M＝1時(shí)，即需已知初始熱流Q0及后r個(gè)溫度測(cè)值T1，T2，…，Tr，再由Q0，T1，T2，…，Tr的信息來確定Q1。由于熱傳導(dǎo)的延遲性，QM決定tM時(shí)刻之后時(shí)間里的溫度值，即TM＋r－1不 僅 取 決 于QM，還 取 決 于QM＋1，…，QM＋r－1，因此，在辨識(shí)QM時(shí)，還必須建立QM＋i（i＝1，r）與QM的關(guān)系。若設(shè)QM－1至QM＋r－1是線性變化的，當(dāng)時(shí)間間隔是常數(shù)時(shí)，則有：

因此，辨識(shí)QM的狀態(tài)方程為：

觀測(cè)方程為（設(shè)采樣為等時(shí)間間隔）：

辨識(shí)的目標(biāo)函數(shù)取為：

視QM為參數(shù)，用靈敏度法來辨識(shí)，則其迭代的牛頓－拉夫遜（Newton－Raphson）算式為：

式中Q的上標(biāo)k、k＋1表示迭代層次，?T／?QM為靈敏度，其滿足的方程可近似寫為：

綜上所述，用順序函數(shù)法辨識(shí)熱流時(shí)，是在時(shí)間方向上按順序向后辨識(shí)。當(dāng)M＝1時(shí)，Q0已知，方程（6）中的T′（x）即為初始時(shí)刻已知的全場(chǎng)溫度分布，由初估的Q1及式（6）可求出后r個(gè)時(shí)刻的溫度測(cè)值T（xm，ti）（i＝1，…，r），再由式（10）求出靈敏度?T（xm，ti）／?QM（i＝1，…，r）后，用式（9）來對(duì)QM優(yōu)化，優(yōu)化得QM的最優(yōu)估值則取為Q1的辨識(shí)值。然后，再沿時(shí)間方向推進(jìn)，用及T2，…，Tr＋1按照同樣的方法來辨識(shí)Q2，依此類推。

在上述辨識(shí)過程中，涉及到狀態(tài)方程式（6）和靈敏度方程式（10）的求解，本文采用有限控制體積法［1，6］來對(duì)其進(jìn)行數(shù)值求解。

3 共軛梯度法

共軛梯度法的基本思想是將傳熱模型的偏微分方程看作為對(duì)未知參數(shù)的約束，利用拉格朗日（Lagrange）乘數(shù)法將辨識(shí)問題轉(zhuǎn)化為使如下目標(biāo)函數(shù)達(dá)極小的無約束優(yōu)化問題［7］：

式中的λ稱為伴隨變量。為了獲得伴隨變量滿足的方程，首先需推導(dǎo)靈敏度方程。

對(duì)式（1）中的Q（t）做擾動(dòng)ΔQ，設(shè)相應(yīng)的溫度場(chǎng)變化為ΔT，熱傳導(dǎo)系數(shù)和比熱變?yōu)椋?/p>

顯然，擾動(dòng)后的物理量也滿足方程式（1）的形式，與未擾動(dòng)時(shí)的方程相減，即可得ΔT滿足的靈敏度方程：

接下來對(duì)式（11）中的目標(biāo)函數(shù)做變分，并引入上面靈敏度方程的邊界條件后，有：

由ΔT的任意性，可得伴隨變量滿足的伴隨方程為：

邊界條件：x＝0：－k（T）

時(shí)域條件：t＝tf，λ＝0。

同時(shí)可求出目標(biāo)函數(shù)對(duì)Q的梯度為：

將此梯度值代入如下優(yōu)化算法，即為共軛梯度法：

其中，γn

U是由式（12）求得的ΔQ＝Pn引起的溫度場(chǎng)變化值；該優(yōu)化算法的計(jì)算停止準(zhǔn)則為：

這一辨識(shí)方法中涉及到的狀態(tài)方程式（6）、靈敏度方程式（12）和伴隨方程（14），同樣都采用有限控制體積法［1，6］來進(jìn)行數(shù)值求解。

4 算例分析

下面給出表面熱流辨識(shí)的算例，對(duì)于如圖1所示的傳熱模型，設(shè)材料密度為ρ＝1650kg／m3，厚度為L(zhǎng)＝5mm，材料熱物性參數(shù)隨溫度變化規(guī)律如圖2示，表面x＝0處的給定熱流如圖3中的“Exact”示，溫度測(cè)點(diǎn)位于x＝L處，利用此給定的熱流可以仿真計(jì)算出測(cè)點(diǎn)溫度歷程如圖4。先不考慮測(cè)量噪聲，將此溫度歷程作為實(shí)測(cè)值，采用順序函數(shù)法和共軛梯度法對(duì)表面熱流進(jìn)行辨識(shí)，順序函數(shù)法中的r取為30，對(duì)應(yīng)的將來時(shí)間長(zhǎng)度tM＋r－1－tM－1＝6s。圖3中示出了采用順序函數(shù)法（圖中“SFM”）和共軛梯度法（圖中“CGM”）辨識(shí)出的表面熱流與給定值的比較，從中可以看到：兩種方法的辨識(shí)結(jié)果與給定值符合很好，兩種辨識(shí)方法都是可行的。

圖2 材料熱物性參數(shù)隨溫度變化函數(shù)Fig.2 Temperature dependent thermal property of the material

圖3 表面熱流辨識(shí)結(jié)果比較Fig.3 Comparison of estimated heat flux

就此問題而言，如果在辨識(shí)過程中不考慮材料熱物性參數(shù)隨溫度的變化，以往的做法是將材料熱物性參數(shù)用近似的常數(shù)值代替后再用線性熱傳導(dǎo)逆問題的表面熱流辨識(shí)算法［8］來處理，這里也給出這一處理方法的結(jié)果。從圖4可以看到，測(cè)點(diǎn)溫度區(qū)間的中值約為T＝550K，可將材料的熱物性參數(shù)取為T＝550K時(shí)對(duì)應(yīng)的常數(shù)值：k＝2.4W／（m·K），Cp＝1800J／（kg·K）；同樣以圖4中的測(cè)點(diǎn)溫度計(jì)算結(jié)果做為測(cè)量值，采用順序函數(shù)法來對(duì)表面熱流進(jìn)行辨識(shí)，得到的表面熱流如圖5中的“k＝2.4，Cp＝1800”所示，可以看到，由于溫度測(cè)量值是在考慮了材料熱物性參數(shù)隨溫度變化的情況下計(jì)算出來的，但此時(shí)辨識(shí)過程中材料熱物性參數(shù)采用了常數(shù)而未考慮隨溫度的變化，因而此時(shí)的辨識(shí)結(jié)果與圖3中的“SFM”辨識(shí)結(jié)果有相當(dāng)大的差異。由此可見，對(duì)于當(dāng)前工況，材料熱物性參數(shù)用近似常數(shù)值代替的辨識(shí)方法得到的表面熱流會(huì)有較大誤差，在辨識(shí)過程中考慮材料熱物性參數(shù)隨溫度的變化是必要的。

圖4 測(cè)點(diǎn)溫度計(jì)算結(jié)果Fig.4 Calculated temperature history at measurement point

圖5 物性參數(shù)選取對(duì)熱流辨識(shí)結(jié)果影響Fig.5 Influence of thermal property on estimated heat flux

接下來在非線性熱傳導(dǎo)逆問題的表面熱流辨識(shí)中考慮測(cè)量噪聲的影響，在圖4的測(cè)點(diǎn)溫度歷程計(jì)算結(jié)果上疊加標(biāo)準(zhǔn)差σ＝2.5K的白噪聲（對(duì)應(yīng)測(cè)量的相對(duì)誤差約為0.5%）來作為模擬的測(cè)量結(jié)果，如圖6中的“Exp.”所示。同樣采用順序函數(shù)法和共軛梯度法來對(duì)表面熱流進(jìn)行辨識(shí)，順序函數(shù)法中的r取為40，對(duì)應(yīng)的將來時(shí)間長(zhǎng)度tM＋r－1－tM－1＝8s。圖7中示出了這兩種方法的辨識(shí)結(jié)果與給定值的比較，圖6中的“Fitted”則示出了采用順序函數(shù)法辨識(shí)結(jié)果計(jì)算出的測(cè)點(diǎn)溫度擬合曲線，可以看到：采用熱流辨識(shí)結(jié)果計(jì)算出的測(cè)點(diǎn)溫度與測(cè)量值擬合較好，表面熱流的辨識(shí)結(jié)果也與給定值符合。

圖6 測(cè)點(diǎn)溫度比較（σ＝2.5K）Fig.6 Comparison of temperature histories at measurement point（σ＝2.5K）

圖7 表面熱流辨識(shí)結(jié)果比較（σ＝2.5K）Fig.7 Comparison of estimated heat flux（σ＝2.5K）

進(jìn)一步增大測(cè)量噪聲的標(biāo)準(zhǔn)差到σ＝10K（對(duì)應(yīng)測(cè)量的相對(duì)誤差約為2%），并將該測(cè)量噪聲疊加到測(cè)點(diǎn)溫度歷程計(jì)算結(jié)果上來模擬測(cè)量結(jié)果，如圖8中的“Exp.”示。圖9中示出了此時(shí)采用順序函數(shù)法和共軛梯度法的辨識(shí)結(jié)果與給定值的比較，順序函數(shù)法中的r取為45，對(duì)應(yīng)的將來時(shí)間長(zhǎng)度tM＋r－1－tM－1＝9s。圖8中的“Fitted”則示出了采用共軛梯度法辨識(shí)結(jié)果計(jì)算出的測(cè)點(diǎn)溫度擬合曲線，可以看到：在相對(duì)誤差約為2%的情況下，采用熱流辨識(shí)結(jié)果計(jì)算出的測(cè)點(diǎn)溫度與測(cè)量值擬合較好，表面熱流的辨識(shí)結(jié)果仍與給定值符合較好，但是由于測(cè)量噪聲的增大，辨識(shí)結(jié)果與給定值之間的偏差有所增大。

圖8 測(cè)點(diǎn)溫度比較（σ＝10K）Fig.8 Comparison of temperature histories at measurement point（σ＝10K）

圖9 表面熱流辨識(shí)結(jié)果比較（σ＝10K）Fig.9 Comparison of estimated heat flux（σ＝10K）

從上述辨識(shí)結(jié)果可以看到：對(duì)于非線性熱傳導(dǎo)逆問題的表面熱流辨識(shí)，順序函數(shù)法和共軛梯度法都能給出較好的辨識(shí)結(jié)果，并且算法都具有較好的抗噪性能。相比較而言：

（1）順序函數(shù)法是對(duì)熱流按時(shí)間順序估計(jì)，算法的推導(dǎo)和實(shí)現(xiàn)較為簡(jiǎn)便；而共軛梯度法是在全時(shí)間域內(nèi)對(duì)熱流一次進(jìn)行估計(jì)，算法的推導(dǎo)和實(shí)現(xiàn)相對(duì)較為復(fù)雜。

（2）順序函數(shù)法中存在人為確定的參數(shù)r，該參數(shù)值的選取需遵循一定的原則，文獻(xiàn)［1，8］對(duì)此進(jìn)行了專門分析。共軛梯度法中雖然沒有自由參數(shù)，但需事先知道式（17）中的測(cè)量結(jié)果標(biāo)準(zhǔn)差σ，以對(duì)優(yōu)化計(jì)算的步數(shù)進(jìn)行控制。實(shí)際上，順序函數(shù)法中確定r的原則也需要用到測(cè)量結(jié)果標(biāo)準(zhǔn)差σ的信息［1］，從這一點(diǎn)來看，兩種辨識(shí)方法是類似的，但具體的要求形式不同。

（3）從計(jì)算結(jié)果的精度來看，兩種辨識(shí)方法的計(jì)算結(jié)果精度大致相當(dāng)，但在有測(cè)量噪聲的情況下，t＜5s段共軛梯度法的辨識(shí)結(jié)果精度明顯低于順序函數(shù)法，其原因尚有待進(jìn)一步深入分析。

（4）從計(jì)算效率看，兩種辨識(shí)方法的計(jì)算效率與溫度測(cè)量的時(shí)間段長(zhǎng)度tf大小有關(guān)，當(dāng)測(cè)量時(shí)間較短時(shí)，順序函數(shù)法的計(jì)算效率較高，共軛梯度法的計(jì)算效率較低；而當(dāng)測(cè)量時(shí)間較長(zhǎng)時(shí)，順序函數(shù)法的計(jì)算效率則相對(duì)較低，共軛梯度法的計(jì)算效率較高。

5 小 結(jié)

本文建立了非線性熱傳導(dǎo)逆問題的兩種表面熱流辨識(shí)方法：順序函數(shù)法和共軛梯度法，給出了這兩種辨識(shí)方法的基本思想和具體算法推導(dǎo)，并針對(duì)典型算例進(jìn)行了仿真辨識(shí)，結(jié)果表明：（1）對(duì)于算例中的工況而言，采用材料熱物性參數(shù)用近似常數(shù)值代替的辨識(shí)方法得到的表面熱流會(huì)有較大誤差，在辨識(shí)過程中考慮材料熱物性參數(shù)隨溫度的變化是必要的。（2）本文針對(duì)非線性熱傳導(dǎo)逆問題建立的兩種表面熱流辨識(shí)方法雖然在算法構(gòu)造、計(jì)算效率方面存在一定的差異，但都能給出較好的辨識(shí)結(jié)果，并且算法受測(cè)量噪聲的影響較小，具有較好的魯棒性。此外，本文工作雖然只討論了一維非線性熱傳導(dǎo)逆問題的表面熱流辨識(shí)算法，但相關(guān)的公式推導(dǎo)和求解方法可以較容易地推廣應(yīng)用于高維情況。

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

［2］OSMAN A M，DOWDING K J，BECK J V.Numerical solution of the general two－dimensional inverse heat conduction problem （IHCP）［J］.JournalofHeat Transfer，1997，119：38－45.

［3］YANG C Y.Estimation of boundary conditions in nonlinear inverse heat conduction problems［J］.Journalof ThermophysicsandHeatTransfer，2003，17（3）：389－395.

［4］HUANG C H，YAN J Y.An inverse problem in simultaneously measuring temperature－dependent thermal conductivity and heat capacity［J］.Int.J.HeatMass Transfer，1995，38（18）：3433－3441.

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

［6］錢煒祺，蔡金獅.用靈敏度法辨識(shí)熱傳導(dǎo)系數(shù)和熱流參數(shù)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1998，16（2）：226－231.

［7］ALIFANOV O M.Inverse heat transfer problems［M］.Springer－Verlag，Berlin，1994.

［8］錢煒祺，何開鋒，桂業(yè)偉，等.非穩(wěn)態(tài)表面熱流反演算法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28（2）：155－161.