充分發(fā)展湍流管道內(nèi)壁邊界的紅外定量識別*

張 林，楊 立，范春利

(海軍工程大學(xué) 動力工程學(xué)院， 湖北 武漢 430033)

許多熱力設(shè)備的管道內(nèi)壁長期受管內(nèi)環(huán)境的影響產(chǎn)生金屬損失、腐蝕等缺陷，對設(shè)備的安全運行帶來隱患[1]。因此，對管道內(nèi)壁邊界的檢測和監(jiān)測具有重要意義[2-4]。

紅外無損檢測技術(shù)是指通過紅外熱像儀得到設(shè)備表面溫度的變化情況來判斷設(shè)備內(nèi)部信息的技術(shù)，由于其非接觸、無損傷、檢測速度快、掃描面積大等獨特的優(yōu)勢已成為一種非常有效的無損檢測技術(shù)[5-7]，被廣泛應(yīng)用于石化、冶金、電力等領(lǐng)域。

傳熱反問題理論是開展內(nèi)部缺陷紅外定量檢測與識別的重要基礎(chǔ)?；趯?dǎo)熱反問題的紅外檢測研究已實現(xiàn)了對熱物性參數(shù)、內(nèi)熱源強度、邊界條件和邊界形狀等參數(shù)的精確檢測和識別[8-11]?；趯α鲹Q熱反問題的研究，由于流動的復(fù)雜性，不適定性嚴(yán)重，尚處于起步階段。湍流由于質(zhì)點相互摻混、碰撞的流動，不適定性更嚴(yán)重，解的存在性與唯一性更難確定，湍流對流換熱反問題的研究比層流對流換熱反問題的研究發(fā)展更慢[12]。

基于對流換熱反問題的管道或平板通道的定量識別已由一些學(xué)者進(jìn)行了研究和討論[13-14]。 Chen等[15-16]使用共軛梯度法估計二維軸對稱層流管道外壁的未知霜層邊界形狀和管道內(nèi)壁上的未知污垢邊界形狀；Rouizi等[17]使用紅外熱像儀測量得到1 mm厚的平板微通道的外表面溫度分布，基于測得的外表面溫度分布對微通道內(nèi)流體的溫度分布進(jìn)行了估計。

在工業(yè)生產(chǎn)中，大多數(shù)管道內(nèi)流體為湍流流動狀態(tài)。隨著工業(yè)生產(chǎn)對設(shè)備安全性和可靠性的要求越來越高，發(fā)展基于對流換熱反問題的湍流管道內(nèi)壁邊界的紅外定量檢測和識別具有非常重要的意義，但是關(guān)于湍流管道內(nèi)壁邊界的識別研究尚未見相關(guān)報道。

在實際中，內(nèi)壁缺陷的存在往往會破壞管道幾何的對稱性，但是直接開展三維管道內(nèi)壁邊界識別研究，管道內(nèi)流體變化復(fù)雜，管壁軸向和周向的傳熱互相影響，反問題研究的不適定性強，內(nèi)壁邊界識別規(guī)律難以分析歸納。參考Chen等[15-16]、肖丹等[18]的研究，首先開展二維軸對稱管道的研究工作，為三維管道內(nèi)部缺陷定量識別提供理論依據(jù)。

本文提出基于管道外壁溫度分布，使用多物理場軟件COMSOL建立充分發(fā)展的湍流管道模型，并與傳熱反問題算法——Levenberg-Marquardt算法相結(jié)合的識別方法，對二維軸對稱充分發(fā)展的湍流管道內(nèi)壁邊界的識別問題進(jìn)行研究。

1 對流換熱正問題

構(gòu)建如圖1所示的分段充分發(fā)展的湍流管道模型。管道內(nèi)部為遠(yuǎn)高于環(huán)境溫度的熱流體，管道置于自然環(huán)境中，通過管道外壁的自然對流和輻射換熱將熱量傳遞給周圍環(huán)境。

對入口段和充分發(fā)展段分別進(jìn)行建模，充分發(fā)展段會出現(xiàn)入出口效應(yīng)，而且管壁軸向?qū)嵩诜侄翁幹袛啵虼嗽O(shè)置虛擬管道前置入口段L1用以保證流體形成充分發(fā)展流，設(shè)置虛擬管道前置段L2和虛擬管道后置出口段L4來避免入出口效應(yīng)和軸向?qū)岬挠绊?。本文研究的管道檢測段為圖中管道黑色區(qū)域L3。R1和R2分別為管道內(nèi)壁和外壁的半徑。入口段L1入口的流體均勻速度和均勻溫度分別為uin和Tin，檢測段L3入口的流體速度和溫度分別為ufu,in和Tfu,in。在工作一段時間后，假定在管道檢測段L3內(nèi)壁上生成了一層腐蝕缺陷，缺陷處的不規(guī)則內(nèi)壁邊界形狀為R(x)。

(a) 入口段(a) Entrance region

(b) 充分發(fā)展段(b) Fully developed region圖1 管道的物理模型Fig.1 Schematic of pipeline

對管道的流動和傳熱特性進(jìn)行數(shù)值模擬，作如下合理假設(shè)[19]：

1)管內(nèi)流體穩(wěn)態(tài)且不可壓縮；

2)流體在邊界上為無滑移流動；

3)忽略重力影響。

則管道模型內(nèi)部流體流動的質(zhì)量守恒方程、動量守恒方程和能量守恒方程為：

ρ·(u)=0

(1)

ρ(u·)u=·{-pI+μ[u+(u)T]}

(2)

ρcpuT=·(kT)

(3)

其中，u是速度矢量，I是單位矩陣，ρ是流體密度，p是壓力，μ是動力黏度，cp是定壓比熱，T是溫度，k是導(dǎo)熱系數(shù)。

管壁穩(wěn)態(tài)導(dǎo)熱的能量守恒方程為：

2T=0

(4)

對流換熱正問題是在管道檢測段L3的內(nèi)壁邊界形狀R(x)和其他所有熱物性參數(shù)、邊界條件均為已知時，求解管道對流換熱問題得到管道檢測段L3外壁溫度分布。

對于充分發(fā)展的湍流管道內(nèi)壁形狀識別問題，管道內(nèi)壁在每次迭代過程中都會發(fā)生變化，因此管道內(nèi)流場不能被指定，需要反復(fù)求解對流換熱問題，計算過程中的模型自動重建、網(wǎng)格自動劃分和求解計算是必需的。

多物理場軟件COMSOL可以在模型計算過程中自動重建幾何模型、更新網(wǎng)格和求解計算，同時COMSOL和MATLAB之間不需要使用輸入/輸出文件就可以進(jìn)行數(shù)據(jù)傳遞[20-21]，可以用來計算許多實際而復(fù)雜的對流換熱問題。因此，本文使用COMSOL對對流換熱正問題進(jìn)行求解，在MATLAB中使用讀取函數(shù)mphinterp直接讀取正問題的外壁溫度分布。

2 反問題識別

對于反問題，不規(guī)則內(nèi)壁邊界R(x)被認(rèn)為是未知的，而式(1)～(4)中其他所有熱物性參數(shù)、邊界條件均為已知的。此外，在管外壁上通過紅外熱像儀采集得到的溫度分布被認(rèn)為是可用的，將測量得到的溫度用Y表示，Y=[Y1,Y2，…，Yi，…，Ym]T，i=1…m，其中m表示測溫點數(shù)目。反問題的求解方法是以式(5)最小化的方式獲得：

S(R)=[Y-T(R)]T[Y-T(R)]

(5)

式中：T(R)為在迭代計算中，根據(jù)識別參數(shù)R求解對流換熱正問題得到的管道檢測段L3外壁溫度分布。

2.1 Levenberg-Marquardt算法

為使式(5)取最小值，通過求導(dǎo)可得[22-23]：

(6)

(7)

式中，

(8)

Ωn=diag[(Jn)TJn]

(9)

n為迭代次數(shù)。

2.2 收斂條件

通過Levenberg-Marquardt算法迭代識別管道檢測段L3內(nèi)壁邊界形狀的收斂條件為：

S(R)<ε

(10)

式中，ε為一個很小的正數(shù)[24]。

2.3 MATLAB調(diào)用COMSOL

MATLAB調(diào)用COMSOL的語句為：

model=mphload(′E:pipeline_heat_thransfer′);

其中，“E:pipeline_heat_thransfer”為求解外壁溫度分布的mph文件的所在路徑，mphload表示MATLAB加載COMSOL的mph文件。

2.4 數(shù)值計算過程

通過Levenberg-Marquardt算法根據(jù)管道外壁溫度分布定量識別管道檢測段L3內(nèi)壁邊界形狀R(x)的迭代求解步驟：

1)給出管道內(nèi)壁邊界形狀R(x)的初始假設(shè)R0作為輸入，開始迭代計算；

2)將估計值Rn傳遞到COMSOL中，MATLAB調(diào)用COMSOL，求解得到管道檢測段L3的外壁溫度分布T(Rn)；

3)根據(jù)式(10)判斷是否收斂，如果收斂則輸出結(jié)果，否則繼續(xù)；

4)采用迭代式(7)修正識別參數(shù)Rn，求得Rn+1返回到步驟2。

3 仿真結(jié)果與討論

數(shù)值實驗中，建立二維軸對稱管道對流換熱模型，管道內(nèi)半徑R1為0.047 m，壁厚0.010 m。依據(jù)管道湍流充分發(fā)展經(jīng)驗公式L/D≥60，管道模型數(shù)據(jù)取為：0～5.64 m為虛擬管道前置入口段L1，0～0.5 m為虛擬管道前置段L2，0.5～1.5 m為管道檢測段L3，1.5～1.8 m為虛擬管道后置出口段L4，管道總長為7.44 m。

管內(nèi)流質(zhì)為空氣，物性參數(shù)取自COMSOL材料庫。管壁導(dǎo)熱系數(shù)為17.6 W·m-1·K-1；外界環(huán)境溫度為25 ℃，管道外壁發(fā)射率為0.95。使用COMSOL 5.2，采用標(biāo)準(zhǔn)k-ε湍流模型，進(jìn)行正問題模型的求解計算。湍流模型收斂條件設(shè)置為1×10-6。

管道內(nèi)壁缺陷只存在于管道檢測段L3，虛擬管道段L1、L2和L4內(nèi)壁均不存在缺陷。在管道檢測段L3外壁表面上均勻布置的測溫點數(shù)為11個，內(nèi)壁邊界均勻分布的離散點亦為11個。算例中，含內(nèi)壁缺陷的管道內(nèi)壁邊界形狀R(x)為：

R(x)=0.047+0.007sin[(x-0.5)π]，0.5≤x≤1.5

(11)

3.1 不規(guī)則內(nèi)壁邊界的紅外特征規(guī)律

對模型的流體域和管道壁固體域均進(jìn)行結(jié)構(gòu)化網(wǎng)格剖分。為保證模型的計算精度，對模型網(wǎng)格進(jìn)行獨立性檢驗。對徑向(流體域+固體域)×(虛擬入口段L1、虛擬前置段L2、檢測段L3和虛擬出口段L4)的網(wǎng)格劃分分別為(12+5)×(141+12+100+7)、(25+5)×(282+25+200+15)和(50+5)×(564+50+400+30)的3套網(wǎng)格進(jìn)行試算。在檢測段入口速度ufu,in=35 m·s-1、入口溫度Tfu,in=573.15 K的情況下，得到管道檢測段L3外壁面中心處的溫度分別為487.15 K，491.68 K，491.91 K，相對誤差分別為0.93%和0.05%。為同時保證計算精度和計算效率，選用第二套網(wǎng)格劃分方案。

使用第二套網(wǎng)格劃分方案，對管道內(nèi)壁形狀為式(11)的管道檢測段L3進(jìn)行網(wǎng)格劃分，如圖2所示?？梢钥闯?，管道檢測段L3管壁中間部分的厚度小于兩端部分的厚度，且從中間部分到兩端部分內(nèi)壁面呈圓滑過渡。

圖2 管道檢測段L3的網(wǎng)格劃分Fig.2 Grid-meshing of pipeline inspection section L3

圖3為管道檢測段L3檢測表面溫度分布和管道檢測段L3入口速度ufu,in、入口溫度Tfu,in的關(guān)系。圖中，檢測表面最大溫差(Maximum Temperature Difference，MTD)是指管道檢測段L3

(a) ufu,in=15 m·s-1, Tfu,in=373.15 K

(b) ufu,in=35 m·s-1, Tfu,in=373.15 K

(c) ufu,in=35 m·s-1, Tfu,in=573.15 K圖3 內(nèi)壁發(fā)生變化前后，檢測表面的溫度分布和 入口速度ufu,in、入口溫度Tfu,in的關(guān)系Fig.3 Relationship between temperature distribution of outer pipeline surface, the inlet velocity ufu,in and the inlet temperature Tfu,in before and after corrosion

外壁溫度分布中最高溫度值和最低溫度值之差；而內(nèi)壁變化后管道檢測段L3外壁溫度分布和內(nèi)壁沒有發(fā)生變化時的溫度分布之差為絕對溫差(Absolute Temperature Difference, ATD)。

由圖3可得，在內(nèi)壁無缺陷時，管壁溫度分布呈線性降低，且入口速度ufu,in越慢，入口溫度Tfu,in越高，換熱越充分，管壁溫度分布最大溫差越大；存在內(nèi)壁缺陷時，管壁溫度分布隨內(nèi)壁缺陷變化，最大溫差略有變化，且入口速度ufu,in越慢，入口溫度Tfu,in越高，管壁溫度分布最大溫差越大。綜上，缺陷的存在與否，不改變外壁溫度分布最大溫差隨入口速度ufu,in、入口溫度Tfu,in的變化規(guī)律。

圖4和圖5分別為管道檢測段L3檢測表面溫度分布絕對溫差和管道檢測段L3入口速度ufu,in、入口溫度Tfu,in的關(guān)系。

(a) Tfu,in=373.15 K

(b) Tfu,in=473.15 K

(c) Tfu,in=573.15 K圖4 檢測表面絕對溫差和入口速度ufu,in的關(guān)系Fig.4 Relationship between absolute temperature difference of outer pipeline surface and inlet velocity ufu,in

(a) ufu,in=5 m·s-1

(b) ufu,in=15 m·s-1

(c) ufu,in=25 m·s-1

(d) ufu,in=35 m·s-1圖5 檢測表面絕對溫差和入口溫度Tfu,in的關(guān)系Fig.5 Relationship between absolute temperature difference of outer pipeline surface and inlet temperature Tfu,in

從圖4中可得：由缺陷引起的絕對溫差在缺陷初始處基本為零，即有缺陷時與無缺陷時的外壁溫度在缺陷初始處差別很小。入口速度ufu,in越快，絕對溫差正增長的幅值越大。在入口速度ufu,in較慢時，由于檢測表面溫度受管壁軸向?qū)岬挠绊?，絕對溫差在前半缺陷處呈下降趨勢，為負(fù)增長；隨著入口速度ufu,in的增大，絕對溫差在前半缺陷處負(fù)增長的幅值逐漸變小直至負(fù)增長消失。綜合絕對溫差中正增長和負(fù)增長的變化規(guī)律，入口速度ufu,in越快，外壁溫度的絕對溫差越大。

從圖5中可以看出：入口溫度Tfu,in越高，絕對溫差正增長的幅值越大；同時，由于檢測表面溫度受管壁軸向?qū)岬挠绊?，隨著入口溫度Tfu,in的升高，絕對溫差在前半缺陷處逐漸形成負(fù)增長，且負(fù)增長的幅值越來越大。綜合絕對溫差中正增長和負(fù)增長的變化規(guī)律，入口溫度Tfu,in越高，外壁溫度絕對溫差越大。

因此，入口速度ufu,in越快，入口溫度Tfu,in越高，由缺陷引起的外壁溫度絕對溫差越大。結(jié)合圖3中的分析可得，在含內(nèi)壁缺陷的高溫管道中，最大溫差和由缺陷引起的絕對溫差并不是同步增加的。

3.2 內(nèi)壁邊界識別

通過在不同入口速度ufu,in和入口溫度Tfu,in下，求解對流換熱正問題得到的管道檢測段L3的外壁溫度分布，來模擬在實際檢測工作中測量得到的外壁溫度分布Y，對式(11)所示的內(nèi)壁缺陷進(jìn)行識別研究。

要說明的是，在實際檢測過程中，熱像儀對管道外壁進(jìn)行檢測工作，得到的是二維熱圖，可以使用熱像儀自帶的軟件進(jìn)行溫度提取工作，使用“線提取”得到管道外壁上沿軸向變化的溫度分布。

圖6和圖7分別為在不同入口速度ufu,in和入口溫度Tfu,in時，不同初始邊界假設(shè)(初始邊界假設(shè)分別為0.047 m和0.055 m)的內(nèi)壁邊界識別結(jié)果(ε=0.000 1)。

從圖6可以看出，入口速度ufu,in越快，不同初值得到的識別結(jié)果曲線越貼合真實缺陷曲線，即識別結(jié)果越好，符合3.1節(jié)圖4的分析，即入口速度ufu,in越快，溫度分布絕對溫差越來越大。由圖7得，隨著入口溫度Tfu,in的升高，識別結(jié)果曲線先是貼近真實缺陷曲線，然后遠(yuǎn)離真實缺陷曲線，即識別結(jié)果先是變好，然后變差，結(jié)合3.1節(jié)圖5(c)的分析，可得識別結(jié)果變好主要是由于溫度分布絕對溫差隨著入口溫度Tfu,in的升高而增大，識別結(jié)果變差主要是由于在入口溫度Tfu,in繼續(xù)升高的過程中溫度分布絕對溫差在前半缺陷處逐漸形成負(fù)增長。

因此，在對流換熱反問題中，絕對溫差大，識別結(jié)果未必好。在進(jìn)行湍流管道的內(nèi)壁邊界識別時，應(yīng)盡可能提高管內(nèi)入口速度，提高絕對溫差。同時，將入口溫度盡可能地控制在一個合理的范圍內(nèi)：入口溫度太低，絕對溫差過小；入口溫度過高，絕對溫差前半部分出現(xiàn)的負(fù)增長會導(dǎo)致識別結(jié)果變差。

(a) ufu,in=15 m·s-1

(b) ufu,in=25 m·s-1

(c) ufu,in=35 m·s-1圖6 入口溫度Tfu,in=573.15 K時內(nèi)壁邊界識別結(jié)果Fig.6 Identification results of inner boundary with inlet temperature Tfu,in=573.15 K

從圖6和圖7亦可以看出，在入口速度較快且入口溫度沒有導(dǎo)致絕對溫差出現(xiàn)負(fù)增長時，不同的初值得到幾乎相同的識別結(jié)果；而在絕對溫差出現(xiàn)負(fù)增長時，不同的初值得到不同的識別結(jié)果，識別結(jié)果有好有壞，即初值對識別結(jié)果的影響大小取決于絕對溫差中是否出現(xiàn)負(fù)增長。

同時還可得，無論哪種入口條件下的識別結(jié)果，前半缺陷的識別精度都很高，而后半缺陷的識別精度變低，在不規(guī)則內(nèi)壁終點處識別精度略微變差。這是一個明顯而且非常有用的經(jīng)驗。在不同的入口條件下識別時，盡管不能識別得到完整的內(nèi)壁缺陷，但能得到前半缺陷的準(zhǔn)確形狀和缺陷的最大值等信息，這對于缺陷檢測工作來說已經(jīng)是很有用的了。

(a) Tfu,in=373.15 K

(b) Tfu,in=473.15 K

(c) Tfu,in=573.15 K圖7 入口速度ufu,in=25 m·s-1時內(nèi)壁邊界識別結(jié)果Fig.7 Identification results of inner boundary with inlet velocity ufu,in=25 m·s-1

4 結(jié)論

針對基于表面測溫的充分發(fā)展湍流管道內(nèi)壁邊界形狀的識別問題，使用COMSOL構(gòu)建分段充分發(fā)展湍流管道模型，關(guān)聯(lián)COMSOL和MATLAB，利用Levenberg-Marquardt算法對二維軸對稱充分發(fā)展的湍流管道的內(nèi)壁邊界形狀的識別進(jìn)行了研究。算例驗證了該方法識別充分發(fā)展湍流管道內(nèi)壁邊界形狀的有效性，可以得出如下結(jié)論：

1)在含內(nèi)壁缺陷的湍流管道中，入口速度越慢，且入口溫度越高，外壁最大溫差越大；入口速度越快，且入口溫度越高，由缺陷引起的外壁絕對溫差越大。因此，外壁最大溫差和由缺陷引起的絕對溫差并不是同步增加的。

2)對流換熱反問題的識別規(guī)律和導(dǎo)熱反問題的識別規(guī)律并不一致，在對流換熱反問題中，由于絕對溫差中負(fù)增長的出現(xiàn)，絕對溫差大，識別結(jié)果未必好。

3)反問題方法的識別精度在含缺陷內(nèi)壁的終點處略微變差，但能準(zhǔn)確地預(yù)測得到內(nèi)壁的其余輪廓，這對管道維護(hù)有很大的幫助。