張經(jīng)豪 盧 濤 熊 平 郝睿智 周照春 田 源 陳柏宇 張琪琪

(北京化工大學(xué) 機(jī)電工程學(xué)院，北京 100029)

引 言

導(dǎo)熱正問題(direct heat conduction problem，DHCP)是根據(jù)邊界條件、初始條件等去計(jì)算導(dǎo)熱物體的整個(gè)溫度場，而導(dǎo)熱反問題(inverse heat conduction problem，IHCP)則是根據(jù)某一部分的溫度場來識(shí)別物體的幾何缺陷[1]、算子[2]、物體內(nèi)部熱源[3]、邊界條件[4]等未知參數(shù)。在核級(jí)管道系統(tǒng)中，很有必要對一些具有特殊安全性的管道進(jìn)行熱疲勞分析。熱疲勞分析需要得到內(nèi)壁面的溫度波動(dòng)，而管道的結(jié)構(gòu)完整性決定了直接測量存在一定困難，由于導(dǎo)熱反問題可用于間接、無損、快速、準(zhǔn)確地獲得內(nèi)壁面的溫度波動(dòng)，對于保證管道的安全運(yùn)行具有重要意義[5]。

導(dǎo)熱反問題的不適定性對此類問題的解決增加了一定的難度，國內(nèi)外學(xué)者所采用的算法也各不相同。Xie等[6]采用控制容積積分法與Levenberg-Marquardt優(yōu)化方法對絕熱材料的熱物性參數(shù)進(jìn)行了瞬態(tài)識(shí)別；張林等[7]利用Levenberg-Marquardt算法識(shí)別了充分發(fā)展的湍流管道內(nèi)壁的幾何邊界；Mohasseb等[8]應(yīng)用遺傳算法反演了方形板邊界熱流密度；錢煒祺等[9]采用順序函數(shù)法與有限控制體積法對二維圓環(huán)域外邊界熱流密度反演問題進(jìn)行了計(jì)算；Chen等[10]運(yùn)用共軛梯度法估算了雙層壁爐內(nèi)壁面的幾何形狀；盧濤等[11]采用共軛梯度法反演了三維T型管內(nèi)壁面的瞬態(tài)溫度。作為一種梯度類算法，共軛梯度法具有抗不適定性的優(yōu)點(diǎn)，因而被廣泛地應(yīng)用于導(dǎo)熱反問題中[12]。熊平等[13]將高斯消元法與共軛梯度法結(jié)合，可較快地反演出一維圓管的內(nèi)壁面溫度；Han等[14]基于共軛梯度法對對流邊界條件進(jìn)行了反演，同時(shí)得出了二維與三維的計(jì)算時(shí)間，但反演所需時(shí)間較長；Xiong等[15]運(yùn)用Gauss-Seidel點(diǎn)迭代法結(jié)合共軛梯度法精確地得出了內(nèi)壁面的瞬態(tài)溫度，但識(shí)別速度較低；后來他們進(jìn)一步改進(jìn)了優(yōu)化方法[16]，提出了一種結(jié)合順序函數(shù)法和共軛梯度法優(yōu)點(diǎn)的序列共軛梯度法來反演待定熱流，有效地提升了反演速率，但反演精度有所降低。

目前大多數(shù)文獻(xiàn)通過研究不同導(dǎo)熱反問題的優(yōu)化方法來提高反演速率，但很少有研究通過改進(jìn)反演算法中的導(dǎo)熱計(jì)算來提高其反演速率。本文通過改進(jìn)正問題的求解方法，將求解一維非穩(wěn)態(tài)導(dǎo)熱問題的托馬斯算法(tridiagonal matrix algorithm，TDMA)與線迭代相結(jié)合，提出TDMA線迭代法，將二維問題轉(zhuǎn)化為多個(gè)一維問題，以提高二維導(dǎo)熱反問題的計(jì)算速率；分別將TDMA線迭代法和常用的Gauss-Seidel點(diǎn)迭代法與適用性較強(qiáng)的共軛梯度法相結(jié)合，通過構(gòu)造數(shù)值試驗(yàn)分析兩種計(jì)算方法反演的精確性和抗噪性，并探討兩種導(dǎo)熱計(jì)算方法對反問題反演速率的影響。

1 物理模型的建立

二維圓管物理模型見圖1。管道外徑R1=30 mm，壁厚d=5 mm。管道導(dǎo)熱系數(shù)k=17.60 W/(m·K)，外壁面對流換熱系數(shù)hf=100 W/(m2·K)，環(huán)境溫度Tf=20 ℃，密度ρ=7 850 kg/m3,比熱容c=465 J/(kg·K)，時(shí)間步長Δt=1 s，總計(jì)算時(shí)間N=30 s。對二維圓管截面在空間上采用均分網(wǎng)格進(jìn)行離散，其中周向均分為36份，徑向均分為10份，網(wǎng)格量為360，網(wǎng)格節(jié)點(diǎn)數(shù)為396個(gè)。

圖1 二維圓管的物理模型和網(wǎng)格數(shù)

2 IHCP數(shù)學(xué)模型

2.1 導(dǎo)熱正問題

圓管的導(dǎo)熱偏微分方程為

(1)

假設(shè)導(dǎo)熱系數(shù)和對流換熱系數(shù)均為常數(shù)。給定內(nèi)壁面為第一類邊界條件，外壁面為第三類邊界條件，則邊界條件為

TS=T(t),S∈wall_in

(2)

(3)

初始條件為

T=T0,t=0

(4)

式中，T為圓管溫度；t為時(shí)間；r為圓管極徑位置；φ為圓管極角位置；Tw為圓管外壁面溫度；S為圓管壁面；Ω為圓管內(nèi)部。

2.1.1離散方法

分別對式(1)～(4)在時(shí)間上和空間上離散，本文采用控制容積積分法進(jìn)行數(shù)學(xué)模型離散，極坐標(biāo)下的網(wǎng)格節(jié)點(diǎn)如圖2所示。在時(shí)間間隔[t,t+Δt]內(nèi)，將式(1)對圖2(a)所示的控制容積作積分可得式(5)；假定非穩(wěn)態(tài)項(xiàng)中的T隨x呈階梯形分布，擴(kuò)散項(xiàng)中T隨t作隱式階梯變化，則可得式(6)；簡化式(6)，省略上標(biāo)t+Δt來表示所求時(shí)刻的節(jié)點(diǎn)溫度，用上標(biāo)0替換上標(biāo)t表示初始時(shí)刻的已知節(jié)點(diǎn)溫度，最終得到內(nèi)節(jié)點(diǎn)離散方程式(7)。同理，對圓管內(nèi)外邊界條件采用相同的方法進(jìn)行數(shù)值離散，可推得內(nèi)邊界離散方程(8)和外邊界離散方程(9)。在給定的定解條件下，即可對上述瞬態(tài)導(dǎo)熱正問題進(jìn)行數(shù)值求解。

圖2 極坐標(biāo)下的網(wǎng)格節(jié)點(diǎn)

(5)

(6)

(7)

TP=T(t)

(8)

(9)

式中，P、N、S、W、E分別表示所研究節(jié)點(diǎn)及其相鄰的4個(gè)節(jié)點(diǎn)；n、e、w、s表示相應(yīng)的界面；Δφ表示相鄰兩節(jié)點(diǎn)間的周向角度；Δr表示相鄰兩節(jié)點(diǎn)間的徑向長度。

2.1.2數(shù)值求解方法

數(shù)值求解方法分為直接解法和迭代解法，本文分別采用Gauss-Seidel點(diǎn)迭代法和TDMA線迭代法進(jìn)行求解。

離散方程可寫成以下通用形式

aPTP=aETE+aWTW+aNTN+aSTS+b

(10)

1)Gauss-Seidel點(diǎn)迭代法

Gauss-Seidel點(diǎn)迭代法的單步迭代通過取鄰點(diǎn)的最新值來代入計(jì)算。假設(shè)先沿著徑向由內(nèi)向外進(jìn)行掃描，再沿著周向從0°位置開始順時(shí)針方向進(jìn)行掃描，Gauss-Seidel迭代式為

(11)

(12)

(13)

式中，n表示本次迭代，n-1表示上一次迭代。

2)TDMA線迭代法

TDMA線迭代法構(gòu)造了一種采用直接法與迭代法相結(jié)合的多維問題代數(shù)方程求解方法，即沿徑向?qū)⒍S問題轉(zhuǎn)化為一維隱式格式，在每個(gè)一維隱式格式的求解中均采用高斯消元法進(jìn)行計(jì)算，將五對角矩陣問題轉(zhuǎn)換為三對角矩陣問題并采用TDMA算法進(jìn)行求解，最終通過沿周向從0°位置開始順時(shí)針方向掃描進(jìn)行線迭代求解。

根據(jù)式(10)，假設(shè)某一極角下的徑向節(jié)點(diǎn)均為待求量，而周向節(jié)點(diǎn)為已知量，設(shè)此時(shí)的常數(shù)項(xiàng)為b′，則離散表達(dá)式為

aPTP=aNTN+aSTS+b′

(14)

式中，b′=aETE+aWTW+b。

將式(14)轉(zhuǎn)換為如下一般方程組形式

AiTi=BiTi+1+CiTi-1+Di

(15)

對于邊界節(jié)點(diǎn)，代數(shù)方程為二元方程；而對于中間節(jié)點(diǎn)，代數(shù)方程均為三元方程，可經(jīng)過逐步消元將中間節(jié)點(diǎn)三元方程轉(zhuǎn)化為二元方程，具體形式如下

Ti-1=Pi-1Ti+Qi-1

(16)

當(dāng)從前向后消元進(jìn)行到最后的邊界方程時(shí)，該二元方程即化為一元方程，可直接求得該邊界節(jié)點(diǎn)的溫度值。聯(lián)立式(15)與式(16)可得出Pi、Qi的遞歸通式，再結(jié)合式(16)從后向前即可依次求出不同徑向位置處的節(jié)點(diǎn)溫度Ti。

將上述求解過程程序化，并在給定定解條件下利用該通用計(jì)算程序?qū)ΧS圓管溫度場進(jìn)行數(shù)值求解，即可獲得管道各節(jié)點(diǎn)的時(shí)域溫度變化。

2.2 導(dǎo)熱反問題

導(dǎo)熱正問題是一個(gè)定解問題，而導(dǎo)熱反問題是是一個(gè)最優(yōu)化問題。因此,本文采用共軛梯度法作為最優(yōu)化問題的求解方法。

2.2.1共軛梯度法

為了獲得精確的內(nèi)壁面溫度，構(gòu)建該最優(yōu)化問題的目標(biāo)函數(shù)為

(17)

式中，M為外壁面測點(diǎn)數(shù)；TI,j,t,cal為反問題求得的外壁面溫度計(jì)算值；TI,j,t,mea為外壁面溫度測量值。

求解敏度問題。設(shè)內(nèi)壁面節(jié)點(diǎn)數(shù)為K，將導(dǎo)熱正問題的控制方程式分別對T0,k,n求偏導(dǎo)，得到敏度系數(shù)求解方程組為

(18)

(19)

(20)

(21)

(22)

內(nèi)壁面溫度迭代式為

(Tk,n)b+1=(Tk,n)b-βb(dk,n)b

(23)

式中，b為迭代步數(shù)；(Tk,n)b+1為迭代計(jì)算得到的新的內(nèi)壁面節(jié)點(diǎn)溫度；βb為迭代步長；(dk,n)b表示迭代搜索方向。βb和(dk,n)b可由式(24)～(26)求得。

(24)

(25)

式中，γb為共軛系數(shù)，且當(dāng)b=0時(shí)，設(shè)定γ0=0。

迭代步長為

βb=

(26)

共軛梯度法的收斂目標(biāo)為

|J(T)b+1-J(T)b|≤μ

(27)

式中，μ為一小正數(shù)。

2.2.2反演流程圖

利用共軛梯度法求解導(dǎo)熱反問題的流程圖如圖3所示。首先，計(jì)算敏度系數(shù)方程，求解0 s穩(wěn)態(tài)導(dǎo)熱反問題，得到初始條件；其次，設(shè)定迭代初始場及迭代步數(shù)b=0，滿足定解條件后進(jìn)行導(dǎo)熱正問題的計(jì)算，得到外壁面溫度值，并求出目標(biāo)函數(shù)；最后根據(jù)是否達(dá)到收斂目標(biāo)來決定返回重新計(jì)算或輸出最終反演值。

圖3 反演流程圖

3 計(jì)算結(jié)果與分析

3.1 導(dǎo)熱反問題算例結(jié)果

通過第2節(jié)所構(gòu)建的二維瞬態(tài)導(dǎo)熱反問題數(shù)學(xué)模型進(jìn)行C語言通用計(jì)算程序的編寫，分別驗(yàn)證基于兩種求解方法(Gauss-Seidel點(diǎn)迭代與TDMA線迭代)導(dǎo)熱反問題的精確性。首先設(shè)定內(nèi)壁面溫度值，并利用導(dǎo)熱正問題程序計(jì)算所得出的外壁面溫度模擬真實(shí)測量值，為導(dǎo)熱反問題分析提供輸入數(shù)據(jù)，所設(shè)定的內(nèi)壁面溫度值作為校驗(yàn)數(shù)據(jù)。分別探討3種不同工況下的導(dǎo)熱反問題，這3種工況設(shè)定的內(nèi)壁面節(jié)點(diǎn)溫度包括隨時(shí)間變化的階躍式函數(shù)(Case 1)、三角形函數(shù)(Case 2)以及同時(shí)隨時(shí)間和空間變化的正弦函數(shù)(Case 3)。

在實(shí)際測溫過程中，測量溫度必定會(huì)存在測量誤差。為驗(yàn)證反演結(jié)果對引入測量誤差的敏感性，本文在導(dǎo)熱正問題計(jì)算得到的外壁面溫度精確值的基礎(chǔ)上引入隨機(jī)誤差，模擬實(shí)際的測量值。

TI,j,t,mea=TI,j,t,exact+σω

(28)

式中，TI,j,t,mea為外壁面溫度測量值；TI,j,t,exact為由正問題計(jì)算得到的外壁面溫度精確值；σ為標(biāo)準(zhǔn)偏差；ω為[-2.576,2.576]區(qū)間內(nèi)服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)，該區(qū)間能達(dá)到99%的測量可靠度。

為了驗(yàn)證內(nèi)壁面溫度反演值與精確值之間的偏離程度，定義平均相對誤差

(29)

式中，T0,j,t,exact為內(nèi)壁面溫度精確值；T0,j,t,est為內(nèi)壁面溫度反演值；ξ越小，則說明反演值與精確值之間的偏離程度越小，反演值越接近于精確值。

3.1.1隨時(shí)間變化的階躍式溫度

設(shè)定內(nèi)壁面階躍式溫度變化為Case 1

(30)

式(30)表示內(nèi)壁面溫度按照階躍式規(guī)律變化。圖4表示當(dāng)內(nèi)壁面呈階躍式溫度變化時(shí)，在不同的標(biāo)準(zhǔn)偏差下，基于兩種求解方法所得的反演溫度值與精確溫度值。由圖4可以看出，當(dāng)標(biāo)準(zhǔn)偏差σ=0和0.1時(shí)，反問題反演溫度曲線與正問題精確溫度曲線幾乎重合。當(dāng)σ=0.2時(shí)，兩種方法的反演溫度均在精確溫度附近上下小范圍波動(dòng)，且Gauss-Seidel點(diǎn)迭代反演溫度的波動(dòng)程度要略明顯于TDMA線迭代。表1為不同標(biāo)準(zhǔn)偏差下的平均相對誤差。由表1可以看出，當(dāng)標(biāo)準(zhǔn)偏差σ=0.2時(shí)，Gauss-Seidel點(diǎn)迭代與TDMA線迭代反演平均相對誤差分別為0.914 01%和0.750 64%，說明二者的反演精度較好且相差較小。當(dāng)標(biāo)準(zhǔn)偏差達(dá)到σ=0.5時(shí)，Gauss-Seidel點(diǎn)迭代與TDMA線迭代反演的平均相對誤差分別為1.793 02%和1.800 00%，說明測量偏差對兩種方法的反演精度都有一定的影響。

圖4 階躍式溫度變化的反演值與精確值

表1 Case 1和Case 2平均相對誤差值

3.1.2隨時(shí)間變化的三角形溫度

設(shè)定內(nèi)壁面三角形溫度變化為Case 2

(31)

根據(jù)式(31)驗(yàn)證內(nèi)壁面溫度呈三角形規(guī)律變化情況下的反演結(jié)果。由圖5可以看出，隨著標(biāo)準(zhǔn)偏差的增大，反演溫度誤差也隨之增大，且在精確值附近上下波動(dòng)。由圖5和表1可得出，當(dāng)標(biāo)準(zhǔn)偏差σ=0時(shí)，兩種方法的反演平均相對誤差在0.000 02%左右，說明在無測量誤差下，兩種方法的反演值非常接近精確值，正問題和反問題所獲得的溫度分布高度相似。當(dāng)σ=0.2時(shí)，應(yīng)用兩種方法得到的反演溫度波動(dòng)不明顯，平均相對誤差較小。當(dāng)標(biāo)準(zhǔn)偏差σ=0.5時(shí)，Gauss-Seidel點(diǎn)迭代法與TDMA線迭代法得到的反演溫度的波動(dòng)幅度相對較大，平均相對誤差分別為1.512 72%和1.611 71%。由上述結(jié)果可以看出平均相對誤差均隨著標(biāo)準(zhǔn)偏差的增大而增大，且在同一標(biāo)準(zhǔn)偏差下，兩種方法的反演誤差較為接近。

圖5 三角形溫度變化的反演值與精確值

3.1.3隨時(shí)間和空間變化的正弦溫度

設(shè)定內(nèi)壁面正弦溫度變化為Case 3

(32)

管道內(nèi)壁面溫度隨時(shí)間和空間呈正弦規(guī)律變化，如圖6所示。根據(jù)式(32)對反演結(jié)果進(jìn)行驗(yàn)證。圖7為兩種計(jì)算方法在給定4個(gè)不同的標(biāo)準(zhǔn)偏差下，對應(yīng)4個(gè)不同位置(圖1所示)處的內(nèi)壁面時(shí)域溫度變化?？梢钥闯鰞煞N方法的反演結(jié)果雖然圍繞精確值波動(dòng)的方向不一致，但波動(dòng)頻率和幅度較為接近。表2為在不同標(biāo)準(zhǔn)偏差下，沿周向均勻分布的4個(gè)角度處的平均相對誤差值。通過對比可進(jìn)一步發(fā)現(xiàn)，在同一標(biāo)準(zhǔn)偏差下兩種方法各個(gè)角度處的平均相對誤差相差較小，并且隨著標(biāo)準(zhǔn)偏差的增大，兩種求解方法在各個(gè)角度處的平均相對誤差也相應(yīng)增大。由表1和表2可得出，當(dāng)標(biāo)準(zhǔn)偏差增大到σ=0.5時(shí)ξ波動(dòng)值仍可以保持在2%左右，說明基于兩種方法的反演精度對測量誤差的變化不敏感，具有一定的抗噪性。

圖6 隨時(shí)間和空間正弦變化的內(nèi)壁面溫度

圖7 隨時(shí)間和空間正弦溫度變化的反演值與精確值

表2 Case 3平均相對誤差值

3.2 程序運(yùn)行時(shí)間對比

表3為3臺(tái)不同PC計(jì)算機(jī)的性能參數(shù)，以同一算例為基礎(chǔ)，采用C語言編譯器VC++6.0，分別在上述3臺(tái)PC計(jì)算機(jī)上進(jìn)行計(jì)算。對基于兩種求解方法的通用計(jì)算程序，采用計(jì)時(shí)函數(shù)clock()精準(zhǔn)計(jì)時(shí)，并將程序運(yùn)行時(shí)間進(jìn)行對比，結(jié)果如表4所示。可以看出，Gauss-Seidel點(diǎn)迭代法程序運(yùn)行所用時(shí)間是TDMA線迭代法程序運(yùn)行時(shí)間的3倍左右，表明基于TDMA線迭代法程序的運(yùn)行速度要明顯優(yōu)于基于Gauss-Seidel點(diǎn)迭代法的程序，具有更高的計(jì)算效率。

表3 各PC計(jì)算機(jī)性能參數(shù)

表4 程序運(yùn)行時(shí)間對比

4 結(jié)論

(1)在不考慮外壁測溫誤差條件下，通過數(shù)值模擬驗(yàn)證了基于控制容積積分法的導(dǎo)熱正問題以及基于共軛梯度法的優(yōu)化算法所構(gòu)建的二維瞬態(tài)導(dǎo)熱反問題數(shù)學(xué)模型的有效性與精確性。研究結(jié)果顯示，Gauss-Seidel點(diǎn)迭代法與TDMA線迭代法反演的內(nèi)壁面溫度均具有較高的識(shí)別精度。

(2)當(dāng)存在外壁測溫誤差時(shí)，探討了反演結(jié)果對測量誤差的敏感性。計(jì)算結(jié)果表明，反演誤差隨著標(biāo)準(zhǔn)偏差的增大而增大，但Gauss-Seidel點(diǎn)迭代法與TDMA線迭代法仍能夠較準(zhǔn)確地識(shí)別出內(nèi)壁面溫度波動(dòng)，即使標(biāo)準(zhǔn)偏差增大到σ=0.5時(shí)，所反演內(nèi)壁面溫度的平均相對誤差仍可保持在2%左右，具有一定的抗噪性。

(3)將基于兩種求解方法的通用程序分別在不同PC計(jì)算機(jī)上運(yùn)行，結(jié)果顯示，相較于常用的Gauss-Seidel點(diǎn)迭代法，TDMA線迭代法在保持求解精度的基礎(chǔ)上，求解速度有明顯的提升，能夠較快地反演得到內(nèi)壁面溫度波動(dòng)，具有更高的計(jì)算效率。

從目前的研究可以看出，三維IHCP與二維IHCP的計(jì)算效率差異很大，三維IHCP所耗費(fèi)的時(shí)間要遠(yuǎn)遠(yuǎn)高于二維IHCP。因此，今后的研究重點(diǎn)是進(jìn)一步探討如何通過改進(jìn)正問題的求解方法來提高三維IHCP的計(jì)算效率。