基于順序函數(shù)法的圓管內(nèi)部流體溫度反演

熊 平，陸 祺，盧 濤，鄧 堅，劉 余，張 勇

(1.北京化工大學(xué) 機(jī)電工程學(xué)院，北京 100029；2.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實驗室，四川 成都 610041)

在反應(yīng)堆停運(yùn)階段，管道內(nèi)的流體處于靜止?fàn)顟B(tài)。在管道外部對流冷卻和流體自然對流作用下，內(nèi)部流體的溫度逐漸降低，并形成熱分層現(xiàn)象。再次啟動反應(yīng)堆時，由于管道內(nèi)部流體溫度未知，會對反應(yīng)堆的安全運(yùn)行造成不良影響，因此需監(jiān)測內(nèi)部流體溫度的變化。然而，由于核電管道完備性的要求，不允許隨意開孔安裝熱電偶，因此提出導(dǎo)熱反問題(IHCP)方法，通過監(jiān)測外壁溫度反演得到內(nèi)部流體溫度。

IHCP方法有梯度類算法和隨機(jī)類算法。梯度類算法包括共軛梯度法(CGM)[1-3]、Tikhonov正則化算法[4]和順序函數(shù)法(SFSM)[5-7]等。隨機(jī)類算法包括遺傳算法[8]、蟻群算法[9]等。盧濤等[1-2]利用CGM對管道內(nèi)壁溫度進(jìn)行了反演，得到了較好的結(jié)果，但并未對管道內(nèi)部流體進(jìn)行討論。錢煒祺等[5]應(yīng)用SFSM對表面熱流密度進(jìn)行了識別，研究了輻射及測量噪聲對反演結(jié)果的影響，并進(jìn)行了初步實驗驗證算法的有效性。孫金金等[10]應(yīng)用隨機(jī)算法中的人工神經(jīng)網(wǎng)絡(luò)算法對建筑墻體的傳熱系數(shù)進(jìn)行了反演識別，并進(jìn)行了4種不同的墻體實驗，驗證了此方法的精確性。

本文利用SFSM對二維圓管內(nèi)近壁面流體溫度和對流換熱系數(shù)進(jìn)行快速反演。通過CFD數(shù)值模擬計算二維圓管在自然對流冷卻作用下管道內(nèi)部流體溫度逐漸降低的過程。提取管道外壁面溫度作為IHCP輸入條件，內(nèi)壁面溫度、近壁面流體溫度及對流換熱系數(shù)作為IHCP的驗證條件，驗證IHCP程序的精確性。

1 物理模型

二維圓管物理模型如圖1所示，管道內(nèi)徑D0=50 mm，外徑D1=60 mm。在自然對流作用下，內(nèi)部流體自然冷卻，給定流體及管道初始溫度T0=200 ℃，外壁環(huán)境溫度為Tf，out=20 ℃。為加快流體溫度的冷卻速率，外壁面對流換熱系數(shù)hout=100 W/(m2·s)。圓管的導(dǎo)熱系數(shù)k=16.30 W/(m·K)，密度ρ=7 900 kg/m3，比熱容c=460 J/(kg·K)。

圖1 物理模型Fig.1 Physical model

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

求解IHCP包括求解導(dǎo)熱正問題和對解的優(yōu)化過程。IHCP根據(jù)測量已知參數(shù)，反演得到未知參數(shù)。為提高反演速率和降低初始參數(shù)對反演參數(shù)的影響，將內(nèi)壁面視為熱流密度邊界條件，利用IHCP反演計算出內(nèi)壁面熱流密度。再利用牛頓冷卻公式與對流換熱準(zhǔn)則關(guān)聯(lián)式，計算內(nèi)部流體溫度和對流換熱系數(shù)。因此，圓管的導(dǎo)熱偏微分方程、初始條件及邊界條件為：

S∈Ω

(1)

(2)

S∈wall_out

(3)

T=T0,t=0

(4)

式中：T為溫度；t為時間；q為圓管內(nèi)壁面?zhèn)鳠釤崃髅芏龋籘w,out為外壁面溫度；r為圓管半徑位置；φ為圓管角度位置；S為圓管壁面；Ω為圓管內(nèi)部；n為壁面法向方向。

分別對式(1)～(4)所構(gòu)成的二維圓管傳熱在時間上和空間上進(jìn)行離散，時間項采用無條件穩(wěn)定的全隱格式進(jìn)行離散。通過網(wǎng)格獨(dú)立性驗證，選取周向網(wǎng)格節(jié)點(diǎn)數(shù)為36，軸向網(wǎng)格節(jié)點(diǎn)數(shù)為11。對以上瞬態(tài)導(dǎo)熱方程采用控制容積積分法進(jìn)行求解，得到圓管各區(qū)域溫度分布。

為提高反演速率，實現(xiàn)內(nèi)部流體溫度的準(zhǔn)在線監(jiān)測和對解的優(yōu)化過程，本文采用Beck等[7]提出的SFSM，并選取未來時間步長R=3，以降低測量噪聲對反演精度的影響。在反演熱流密度時，還須建立未來時間步長下熱流密度的關(guān)系，本文假設(shè)熱流密度在R-1個未來時刻相等，即：

qn,m=qn+1,m=…=qn+R-1,m

(5)

式中：n為反演時間；m為反演參數(shù)位置。

其最優(yōu)化的目標(biāo)函數(shù)為：

(6)

式中：M為測點(diǎn)數(shù)；Ti，j，cal(q)為測點(diǎn)處的計算溫度；Ti，j，mea為測量溫度。

迭代優(yōu)化計算公式為：

(7)

(8)

(9)

S∈wall_out

(10)

(11)

順序函數(shù)法的收斂目標(biāo)為：

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

(12)

式中，μ為很小的正數(shù)，本文選取μ≤10-6。

通過圓管IHCP反演得到圓管內(nèi)壁面局部傳熱熱流密度后，由內(nèi)壁面牛頓冷卻公式及自然對流換熱準(zhǔn)則關(guān)聯(lián)式，計算內(nèi)部近壁面流體溫度和對流換熱系數(shù)。

牛頓冷卻公式為：

q=hinΔT=hin(Tw-Tf)

(13)

式中：hin為流體與壁面間的局部對流換熱系數(shù)；Tw為圓管內(nèi)壁面溫度；Tf為局部近壁面流體溫度。

根據(jù)努賽爾準(zhǔn)則：

(14)

自然對流換熱準(zhǔn)則關(guān)聯(lián)式[11]為：

(15)

式中：Gr為格拉曉夫數(shù)；Pr為普朗特數(shù)；g為重力加速度；α為流體體積膨脹系數(shù)；ΔT為Tw和Tf之差；υ為運(yùn)動黏度；C和n為實驗確定的常數(shù)。

關(guān)聯(lián)式中的流體定性溫度為反演流體溫度Tf和對應(yīng)內(nèi)壁面溫度Tw的平均溫度。本研究中為了加快管道內(nèi)流體冷卻速率，提高了外壁面對流換熱系數(shù)，這必然對內(nèi)部流體自然對流換熱造成影響，將直接影響選取的實驗常數(shù)的大小。通過數(shù)值實驗，選擇C=0.40，n=1/4。

將式(14)、(15)聯(lián)立求解，hin表示為ΔT的函數(shù)，代入到式(13)中。由于Tw和q已由導(dǎo)熱反問題求出，故式(14)則轉(zhuǎn)變?yōu)殛P(guān)于Tf的方程，從而求得管內(nèi)近壁面流體溫度Tf，再根據(jù)式(15)可求得對流換熱系數(shù)hin。基于SFSM求解管內(nèi)近壁面流體溫度和對流換熱系數(shù)的具體流程如圖2所示。

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

通過數(shù)值仿真實驗驗證反問題模型的準(zhǔn)確性，并探討測量誤差對反演結(jié)果的影響。在導(dǎo)熱反演算法中，計算時間步長Δt=1 s，計算總時間步長N=200 s。

為驗證測量噪聲對反演精度的影響，引入服從正態(tài)分布的測量誤差。

Tmea=Texact+σω

(16)

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

圖2 反演流程圖Fig.2 Chart of inversion progress

為驗證該反演模型的反演精度，定義平均相對誤差ξ為：

(17)

式中：Ti,exact為流體溫度精確值，即數(shù)值模擬值；Ti,est為流體溫度反演值。ξ越小，則說明反演值與精確值偏離程度越小，反演值越接近于精確值。

3.1 CFD數(shù)值模擬結(jié)果

通過CFD數(shù)值模擬，為導(dǎo)熱反問題分析提供輸入數(shù)據(jù)(外壁面溫度)和校驗數(shù)據(jù)(近壁面流體溫度和對流換熱系數(shù))。圖3示出不同時刻流體及管道數(shù)值模擬溫度云圖。在流體自然對流和浮升力作用下，流體溫度逐漸降低，并呈現(xiàn)上熱下冷的熱分層現(xiàn)象。在近壁面附近，流體存在明顯的熱邊界層。提取管道外壁面模擬溫度作為反演程序的測量溫度。

圖4示出靠近壁面1、2和3 mm位置處的流體溫度隨時間的變化。由圖4可看出，除180°位置，其余位置處距壁面2 mm和3 mm的流體溫度曲線幾乎重合。然而，在1 mm位置處溫度存在較大的階躍變化，這主要是因為在靠近壁面邊界處，在自然對流作用下流體往下流動，并出現(xiàn)湍動現(xiàn)象，如圖3所示。在圓管底部(180°位置)位置，流體存在明顯的熱分層現(xiàn)象。因此在180°位置處，將距壁面1 mm位置的流體溫度作為近壁面流體溫度；在其他位置處，將距壁面2 mm位置的流體溫度作為近壁面流體溫度。

a——t=50 s；b——t=100 s；c——t=150 s；d——t=200 s圖3 管道及流體溫度分布云圖Fig.3 Temperature distribution cloud of pipe and fluid

位置：a——0°；b——90°；c——180°；d——270°圖4 不同位置流體溫度分布Fig.4 Fluid temperature distribution at different locations

3.2 導(dǎo)熱反演結(jié)果

圖5示出在無測量噪聲條件下，反演溫度與數(shù)值模擬溫度比較。由圖5可見，內(nèi)壁面溫度反演值與模擬值具有很好的吻合度，而流體溫度反演值與模擬值的吻合度要稍差些。這主要由兩方面造成：一是數(shù)值模擬中提取的流體溫度為靠近內(nèi)壁面不同位置處，用該位置表征近壁面流體溫度，造成反演值與模擬值的差別；二是選擇的實驗常數(shù)C和n造成的反演誤差。在無測量誤差(σ=0)時，各位置處近壁面流體溫度的反演平均相對誤差均小于1%，如表1所列。說明該導(dǎo)熱反問題對內(nèi)壁面流體溫度的反演具有較高精度。

圖6示出在無測量噪聲條件下，不同位置內(nèi)部流體對流換熱系數(shù)反演值與模擬值隨時間的變化。由于初始溫度的影響，開始階段對流換熱系數(shù)有失真現(xiàn)象。隨反演時間的推移，反演對流換熱系數(shù)在模擬值的平均值附近趨于平穩(wěn)，這主要是因為在數(shù)值模擬中，計算時間間隔為0.005 s，計算得到的對流換熱系數(shù)每隔1 s(200個時間步)保存1次，近似為一瞬時對流換熱系數(shù)。而導(dǎo)熱反問題中，計算時間間隔為1 s，每個時間步保存對流換熱系數(shù)，反演得到的對流換熱系數(shù)為該時間間隔內(nèi)的平均對流換熱系數(shù)。因此，反演的對流換熱系數(shù)相對平穩(wěn)，并趨于數(shù)值模擬的平均值附近。在圓管底部(180°)位置，該位置具有明顯的熱分層，自然對流換熱波動已不再明顯。因此，數(shù)值模擬計算得到的對流換熱系數(shù)相對平穩(wěn)。在圓管底部(180°)位置，對流換熱系數(shù)反演值與模擬值差別較大，這主要是由于數(shù)值模擬得到的對流換熱系數(shù)是通過提取的內(nèi)壁面局部熱流密度和近壁面溫度計算得到，造成反演值與模擬值有一定差別。不論是反演結(jié)果還是數(shù)值模擬結(jié)果，圓管底部位置的對流換熱系數(shù)均低于其他位置。這是因為在浮升力的作用下，底部流體溫度較低，與壁面換熱強(qiáng)度降低，使得底部對流換熱系數(shù)低于其他位置。

位置：a——0°；b——90°；c——180°；d——270°圖5 溫度反演值與模擬值的比較Fig.5 Comparison of temperature between inversion value and simulation value

表1 流體溫度反演平均相對誤差Table 1 Mean relative error of fluid temperature inversion value

圖7示出不同測量誤差條件下內(nèi)部近壁面流體溫度的反演結(jié)果。由圖7可見，反演結(jié)果在精確值之間來回波動，隨測量誤差的增大，其波動幅度增大。當(dāng)標(biāo)準(zhǔn)偏差σ=0.1時，各位置處近壁面流體溫度的反演值與模擬值之間仍符合得較好，波動的平均相對誤差均為1%左右。說明該反演程序在存在測量誤差條件下仍具有較好的反演精度。

圖8示出存在測量噪聲條件下內(nèi)部流體的對流換熱系數(shù)的反演結(jié)果。與反演流體溫度結(jié)果類似，反演值隨測量誤差的增大，波動幅度增大，尤其在圓管底部位置。表2列出對流換熱系數(shù)反演的平均相對誤差。由表2可見，在有無測量誤差條件下，其對流換熱系數(shù)的反演平均相對誤差變化較小，測量誤差對對流換熱系數(shù)的反演結(jié)果影響較小。說明該反演模型具有較好的抗噪性。

基于未來時間步長R=3的SFSM計算二維圓管導(dǎo)熱模型，計算總時間N=200 s，計算時間間隔Δt=1 s，總耗時約10 s，極大地提高了計算效率，能實現(xiàn)對內(nèi)部流體溫度和對流換熱系數(shù)的快速反演。

位置：a——0°；b——90°；c——180°；d——270°圖6 對流換熱系數(shù)反演值與模擬值的比較Fig.6 Comparison of convective heat transfer coefficient between inversion value and simulation value

位置：a——0°；b——90°；c——180°；d——270°圖7 測量噪聲對反演近壁面流體溫度的影響Fig.7 Influence of measurement noise on inversion of near-wall fluid temperature

位置：a——0°；b——90°；c——180°；d——270°圖8 測量噪聲對反演對流換熱系數(shù)的影響Fig.8 Influence of measurement noise on inversion of convective heat transfer coefficient

表2 對流換熱系數(shù)反演平均相對誤差Table 2 Mean relative error of convective heat transfer coefficient inversion value

4 結(jié)論

本文對二維圓管內(nèi)流體自然對流冷卻過程的近壁面流體溫度和對流換熱系數(shù)同時反演。建立了基于SFSM的快速反演方法，采用控制容積積分法求解導(dǎo)熱正問題。通過數(shù)值仿真實驗驗證導(dǎo)熱反問題程序的精確性。數(shù)值實驗結(jié)果表明，通過牛頓冷卻公式和自然對流關(guān)系式，利用SFSM能較好地反演得到圓管近壁面流體溫度和對流換熱系數(shù)。通過引入高斯分布的隨機(jī)測量誤差，探討了該反演模型對測量誤差的敏感性。在有、無測量誤差條件下，近壁面流體溫度的反演平均相對誤差均在1%左右，而對流換熱系數(shù)的反演精度略差一些，其平均相對誤差均小于20%。從反演結(jié)果可看出，該反演模型具有較好的抗噪性。