張經(jīng)豪 熊 平 郝睿智 盧 濤

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

引 言

在工業(yè)生產(chǎn)中，管道廣泛應(yīng)用于輸氣、供水和運(yùn)油等領(lǐng)域。由于管道長(zhǎng)期處于高溫高壓的工作狀態(tài)，在使用過程中容易造成內(nèi)壁的腐蝕、脫落，甚至產(chǎn)生污垢黏連附著在內(nèi)壁上導(dǎo)致出現(xiàn)內(nèi)壁增厚等現(xiàn)象，這些缺陷的存在成為安全生產(chǎn)和使用的潛在隱患[1]。因此，對(duì)管道內(nèi)壁缺陷的檢測(cè)具有重要意義[2-3]。

導(dǎo)熱反問題 (inverse heat conduction problem，IHCP)已在很多工程領(lǐng)域中得到廣泛應(yīng)用[4]。導(dǎo)熱反問題常用于計(jì)算反演邊界條件[5-9]、熱物性參數(shù)[10]、物體內(nèi)部熱源[11]及初始條件等未知項(xiàng)，也可用于物體內(nèi)部缺陷[12-15]和內(nèi)壁缺陷的識(shí)別。Chen等[16-17]采用共軛梯度法估計(jì)了管道外壁的未知霜層邊界形狀和管道內(nèi)壁上的未知污垢邊界形狀，但只單一討論了沿管道軸向變化的霜層或污垢分布情況；文獻(xiàn)[18-19]對(duì)雙層爐膛內(nèi)壁的邊界形狀進(jìn)行了反演，但沒有探討初始值及測(cè)點(diǎn)數(shù)目對(duì)反演結(jié)果精確度的影響；文獻(xiàn)[20-24]應(yīng)用共軛梯度法對(duì)管道內(nèi)壁面的邊界形狀進(jìn)行了穩(wěn)態(tài)或瞬態(tài)的識(shí)別，并分別對(duì)所建反演模型的可行性進(jìn)行了探究；張林等[25]在改進(jìn)的一維修正算法基礎(chǔ)上提出一維加權(quán)法，對(duì)二維圓管交界面形狀進(jìn)行了定量識(shí)別，提高了識(shí)別效率；Mohasseb 等[26]采用遺傳算法作為非適定導(dǎo)熱反問題的計(jì)算方法，求解了方形板的熱流密度邊界，但算法的穩(wěn)定性與收斂速度有待提高；文獻(xiàn)[3]、[27-30]均采用列文伯格- 馬夸爾特(Levenberg- Marquardt，L- M)算法對(duì)未知參數(shù)進(jìn)行了計(jì)算，算例驗(yàn)證了算法的有效性和優(yōu)越性。在眾多求解導(dǎo)熱反問題的數(shù)值計(jì)算方法中，屬于梯度類算法的L- M算法因具有穩(wěn)定性好、構(gòu)造思路簡(jiǎn)單、收斂快等優(yōu)勢(shì)而被廣泛應(yīng)用。

當(dāng)前的大多數(shù)研究都只對(duì)管道內(nèi)壁腐蝕減薄缺陷進(jìn)行單一的討論，涉及管道內(nèi)壁污垢增厚缺陷識(shí)別的相對(duì)較少，因此本文提出基于有限元法和L- M算法構(gòu)建反問題數(shù)學(xué)模型，利用COMSOL與MATLAB聯(lián)合仿真穩(wěn)定、快速的優(yōu)勢(shì)，將管道內(nèi)徑作為反演參數(shù)，分別對(duì)管道內(nèi)壁腐蝕減薄缺陷與污垢增厚缺陷進(jìn)行定量識(shí)別，并構(gòu)造數(shù)值實(shí)驗(yàn)對(duì)反演結(jié)果的精確性和穩(wěn)定性進(jìn)行分析。

1 二維圓管物理模型

二維圓管物理模型如圖1所示。管道外徑a=0.3 m，壁厚d=25 mm。管道的導(dǎo)熱系數(shù)λ1=17.60 W/(m·K)，污垢的導(dǎo)熱系數(shù)λ2=3.14 W/(m·K)，外壁面對(duì)流換熱系數(shù)hout=10 W/(m2·K)，環(huán)境溫度Ta=25 ℃，內(nèi)壁面對(duì)流換熱系數(shù)hin=1 000 W/(m2·K)，管內(nèi)流體溫度Tf=200 ℃，θ為極角，r為出現(xiàn)缺陷后管道的內(nèi)壁面極徑。

圖1 管道兩類缺陷的物理模型Fig.1 Physical model of two types of defects in pipelines

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

2.1 導(dǎo)熱正問題(DHCP)

由于所研究?jī)?nèi)容為二維圓管內(nèi)壁面缺陷的穩(wěn)態(tài)識(shí)別，基于此作如下合理假設(shè)：1)管道忽略軸向?qū)岬挠绊懀?)管道導(dǎo)熱系數(shù)和管內(nèi)外的對(duì)流換熱系數(shù)均為常數(shù)。

2.1.1管道增厚缺陷

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

(1)

假設(shè)給定內(nèi)壁面和外壁面均為第三類邊界條件，則邊界條件為

(2)

(3)

污垢增厚缺陷存在著交界面Γ，在此處

T1=T2

(4)

(5)

式中，下標(biāo)1代表管壁區(qū)域Ω1，下標(biāo)2代表污垢區(qū)域Ω2；下標(biāo)in、out分別為圓管內(nèi)壁面和外壁面；n代表法線方向；下標(biāo)int代表圓管與污垢交界面邊界。

2.1.2管道減薄缺陷

管道減薄缺陷只存在于管壁區(qū)域Ω1，將式(3)進(jìn)行修改可得式(6)。

(6)

式(6)與式(1)、式(2)共同構(gòu)成了管道減薄缺陷導(dǎo)熱正問題的數(shù)學(xué)描寫。

2.2 導(dǎo)熱反問題

導(dǎo)熱正問題是一個(gè)定解問題，而導(dǎo)熱反問題則是一個(gè)最優(yōu)化問題。

2.2.1L-M算法

L- M算法是非線性回歸中回歸參數(shù)最小二乘估計(jì)的一種方法，L- M算法將最速下降法和線性化方法相結(jié)合，從而較快地找到最優(yōu)值。

對(duì)于兩類內(nèi)壁面缺陷的反演，內(nèi)壁面幾何邊界都是未知的，可由極坐標(biāo)向量R(θ)描述

R(θ)=[r1,r2,…,rn,…,rN]T

(7)

式中，r為極半徑；N為內(nèi)壁面離散節(jié)點(diǎn)個(gè)數(shù)。

數(shù)學(xué)偏微分方程組(式(1)～(5)以及式(1)、(2)、(6))需要獲得外壁面測(cè)點(diǎn)溫度以滿足其封閉性，離散后的外壁面溫度測(cè)量值記為

Y=[Y1,Y2,…,Ym,…,YM]T

(8)

式中，M為外壁面測(cè)點(diǎn)數(shù)。

進(jìn)一步構(gòu)建該最優(yōu)化問題的目標(biāo)函數(shù)為

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

(9)

式中，T(R)為根據(jù)R(θ)所求得的外壁面溫度的計(jì)算值，表示為

T(R)=[T1,T2,…,Tm,…,TM]T

(10)

通過式(9)對(duì)未知參量R(θ)求偏導(dǎo)得到目標(biāo)函數(shù)梯度的最小化條件為

(11)

式中，J(R)為敏度系數(shù)矩陣，表示為

(12)

(13)

式中，k為迭代次數(shù)；Ωk為一對(duì)角矩陣，表示為

Ωk=diag[(Jk)TJk]

(14)

根據(jù)偏差原理給出收斂條件為

S(R)<ε

(15)

式中，ε在不考慮測(cè)量誤差時(shí)為一小正數(shù)。

2.2.2數(shù)值計(jì)算過程

MATLAB使用mphload函數(shù)加載COMSOL的mph文件，通過mphinterp函數(shù)便可讀取數(shù)據(jù)，無需使用輸入/輸出文件。COMSOL在被調(diào)用時(shí)可以自動(dòng)構(gòu)建幾何模型、更新網(wǎng)格和求解計(jì)算，有效提高了反演效率。

利用L- M算法求解導(dǎo)熱反問題的流程圖如圖2所示。

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

2.2.3反演識(shí)別誤差

為驗(yàn)證反演結(jié)果對(duì)引入測(cè)量誤差的敏感性，在由導(dǎo)熱正問題計(jì)算得到的外壁面精確值的基礎(chǔ)上引入隨機(jī)誤差，來模擬實(shí)際的測(cè)量值。

Tmea=Texact+σω

(16)

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

當(dāng)考慮實(shí)際測(cè)溫誤差時(shí)，收斂標(biāo)準(zhǔn)ε寫成

(17)

為了表征反演值與精確值的偏離程度，定義平均相對(duì)誤差為

(18)

式中，rn,exact為內(nèi)徑精確值；rn,est為內(nèi)徑反演值。

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

通過運(yùn)行程序進(jìn)行數(shù)值計(jì)算，以驗(yàn)證管道減薄缺陷與增厚缺陷導(dǎo)熱反問題反演的精確性。首先設(shè)定內(nèi)壁面缺陷形狀，由導(dǎo)熱正問題得出外壁面溫度模擬的實(shí)際測(cè)量值，為導(dǎo)熱反問題提供輸入數(shù)據(jù)；所設(shè)定的內(nèi)壁面缺陷形狀參數(shù)作為校驗(yàn)數(shù)據(jù)。為驗(yàn)證程序的穩(wěn)定性和抗噪性，分別探討初值、測(cè)溫點(diǎn)數(shù)目及測(cè)量誤差對(duì)反演結(jié)果的影響。

3.1 管道減薄缺陷

對(duì)于管道減薄缺陷，在外壁面上均勻設(shè)置36個(gè)測(cè)點(diǎn)，內(nèi)壁面由36個(gè)均勻分布的離散節(jié)點(diǎn)擬合得到。設(shè)定沿周向變化的內(nèi)壁面幾何形狀為階躍式函數(shù)、鍥形函數(shù)和正弦函數(shù)，分別探討這3種不同工況下的導(dǎo)熱反問題。

3.1.1階躍式函數(shù)

設(shè)定內(nèi)壁面邊界形狀按階躍式變化，即Case 1。

(19)

圖3 不同初值下階躍式幾何邊界的反演值與精確值Fig.3 Inversion and exact values for a stepped geometrical boundary with different initial values

根據(jù)式(19)驗(yàn)證內(nèi)壁面幾何邊界呈階躍式規(guī)律變化下的反演結(jié)果，并探討不同初值R0對(duì)反演結(jié)果的影響。由圖3可以看出，當(dāng)R0取不同值時(shí)，反演曲線與精確曲線的走向基本趨于一致，靠近階躍邊界處的識(shí)別效果相對(duì)較差，但程序仍可較好地識(shí)別出管道內(nèi)壁的邊界形狀。由表1可以進(jìn)一步得出，3種初始邊界假設(shè)下，反演結(jié)果的平均相對(duì)誤差最大僅相差0.015%，說明初始邊界假設(shè)對(duì)內(nèi)壁面幾何邊界的反演幾乎沒有影響；此外，3種初值下得出計(jì)算結(jié)果的平均相對(duì)誤差均不超過0.1%，驗(yàn)證了程序反演的精確性。

表1 Case1～Case3平均相對(duì)誤差值

3.1.2鍥形函數(shù)

設(shè)定內(nèi)壁面幾何邊界形狀以鍥形函數(shù)變化，即Case 2。

(20)

在式(20)所示的內(nèi)壁邊界形狀以鍥形函數(shù)變化的條件下，研究不同外壁面溫度測(cè)點(diǎn)數(shù)目對(duì)識(shí)別結(jié)果的影響。由圖4可以看出，當(dāng)M=24時(shí)，程序?qū)?nèi)壁邊界的識(shí)別效果相對(duì)較差，尤其是在階躍點(diǎn)附近更為明顯；隨著外壁面溫度測(cè)點(diǎn)數(shù)目的增加，階躍點(diǎn)處的識(shí)別效果逐漸變好，當(dāng)M=72時(shí)，反演曲線與精確曲線高度吻合。結(jié)合表1可進(jìn)一步得出，外壁面測(cè)溫點(diǎn)數(shù)越多，識(shí)別效果就越好。M=36時(shí)的平均相對(duì)誤差為0.058%，M=72時(shí)的平均相對(duì)誤差為0.011%，均可較好地識(shí)別出內(nèi)壁邊界形狀。

圖4 鍥形幾何邊界不同測(cè)溫點(diǎn)數(shù)時(shí)的反演值與精確值Fig.4 Inversion and exact values for a wedge-shaped geometric boundary with different numbers of measurement points

3.1.3正弦函數(shù)

設(shè)定內(nèi)壁邊界形狀隨正弦函數(shù)變化，即Case 3。

R(θ)=0.285+0.01sinθ

(21)

內(nèi)壁面邊界形狀按式(21)以正弦規(guī)律變化，引入標(biāo)準(zhǔn)偏差，分析不同偏差對(duì)反演結(jié)果的影響。如圖5所示，可以看出當(dāng)σ=0 ℃時(shí)，反問題得到的反演曲線與正問題得出的精確曲線幾乎重合；當(dāng)σ=0.2 ℃和σ=0.5 ℃時(shí)，兩條反演曲線均在精確值附近上下小范圍波動(dòng)，且在σ=0.5 ℃時(shí)變化得相對(duì)更為明顯。根據(jù)表1不難進(jìn)一步得出，隨著標(biāo)準(zhǔn)偏差的增大，反演誤差也隨之增大，但3種標(biāo)準(zhǔn)偏差下的平均相對(duì)誤差均不超過1.5%，說明采用此方法反演邊界形狀對(duì)測(cè)量誤差的變化不敏感，具有一定的抗噪性。圖6是當(dāng)邊界形狀為正弦函數(shù)時(shí)的迭代收斂歷程圖，由圖可看出，迭代步長(zhǎng)先大后小，且只需要迭代3次即可滿足所設(shè)定的收斂條件。對(duì)比圖6(c)與圖6(d)可知，運(yùn)用L- M算法可在低迭代步數(shù)下準(zhǔn)確地識(shí)別出缺陷。

圖5 正弦?guī)缀芜吔绮煌瑴y(cè)量誤差時(shí)的反演值與精確值Fig.5 Inversion and exact values of a sine geometric boundary with different measurement errors

圖6 正弦函數(shù)邊界缺陷反演時(shí)的收斂歷程Fig.6 Convergence history of sine function boundary defect inversion

3.2 管道增厚缺陷

同樣地，對(duì)于管道增厚缺陷，在外邊界上也均勻設(shè)置36個(gè)測(cè)點(diǎn)，含有污垢的內(nèi)壁面亦由36個(gè)均勻分布的離散節(jié)點(diǎn)擬合得到，可通過反演管道內(nèi)徑進(jìn)而得出污垢的厚度。設(shè)定的內(nèi)壁面污垢幾何形狀包括均勻增厚缺陷、沿周向變化的三角形函數(shù)以及橢圓函數(shù)，分別探討這3種不同工況下的導(dǎo)熱反問題。

3.2.1均勻增厚

設(shè)定內(nèi)壁面有均勻厚度的污垢，即Case 4。

R(θ)=0.265

(22)

如式(22)所示，假設(shè)周向含有均勻污垢的圓管內(nèi)徑為0.265 m，即污垢的厚度為0.01 m，分別設(shè)置不同的初值R0，所得污垢厚度Δd反演結(jié)果如圖7所示。由圖7可以看出，不同初值下的反演值與精確值相差甚小，且反演曲線呈小范圍波動(dòng)。由表2可得，3種初始值下反演的平均相對(duì)誤差均不足0.1%，說明初值對(duì)均勻污垢反演結(jié)果的影響較小，程序?qū)τ趦?nèi)壁面均勻污垢有較好的識(shí)別能力。

圖7 不同初值下均勻污垢厚度的反演值與精確值Fig.7 Inversion and exact values of uniform fouling with different initial values

表2 Case4～Case6平均相對(duì)誤差值

3.2.2三角形函數(shù)

設(shè)定內(nèi)壁幾何邊界形狀按三角形函數(shù)變化，即Case 5。

R(θ)=

(23)

如圖8所示，含有污垢的管道內(nèi)壁面半徑(即內(nèi)壁污垢形狀)隨周向呈三角形規(guī)律變化。根據(jù)式(23)對(duì)反演結(jié)果進(jìn)行驗(yàn)證，同時(shí)探討改變外壁面測(cè)溫點(diǎn)個(gè)數(shù)對(duì)反演結(jié)果的影響。由表2可知，增加外壁面測(cè)溫點(diǎn)個(gè)數(shù)可以更加準(zhǔn)確地識(shí)別出內(nèi)壁面的污垢形狀，但提高的程度有限，當(dāng)測(cè)溫點(diǎn)的個(gè)數(shù)達(dá)到一定值時(shí)，增加節(jié)點(diǎn)個(gè)數(shù)對(duì)識(shí)別精度的提高影響并不大。另外，由表2結(jié)果也可看出，當(dāng)測(cè)溫點(diǎn)個(gè)數(shù)較少時(shí)識(shí)別精度仍然較高，M=24時(shí)的平均相對(duì)誤差為0.033%，表明仍能較準(zhǔn)確地識(shí)別出內(nèi)壁面的污垢形狀。

圖8 不同測(cè)溫點(diǎn)數(shù)下三角形污垢厚度的 反演值與精確值Fig.8 Inversion and exact values of triangle fouling with different numbers of temperature measurement points

3.2.3橢圓函數(shù)

設(shè)定內(nèi)壁面幾何邊界形狀以橢圓函數(shù)變化，即Case 6。

(24)

式(24)為內(nèi)壁邊界形狀隨橢圓函數(shù)變化的情況，即內(nèi)壁污垢形狀亦呈橢圓函數(shù)變化。圖9為不同標(biāo)準(zhǔn)偏差下污垢尺寸的反演值與精確值。由圖可看出，無偏差下反演值與精確值二者的曲線走向基本同步；當(dāng)存在標(biāo)準(zhǔn)偏差時(shí)，反演值圍繞精確值附近上下波動(dòng)，且隨著標(biāo)準(zhǔn)偏差的增大，波動(dòng)的幅度和頻率也愈為明顯。根據(jù)表2可知，在標(biāo)準(zhǔn)偏差增大到0.5 ℃時(shí)，波動(dòng)值仍可保持在1%以內(nèi)，說明反演模型在識(shí)別內(nèi)壁污垢時(shí)穩(wěn)定性較好。圖10為當(dāng)污垢以橢圓函數(shù)變化時(shí)的管道迭代收斂歷程圖。由圖可以看出，對(duì)橢圓污垢的反演與對(duì)正弦減薄缺陷的反演一樣，僅是經(jīng)歷3次迭代就達(dá)到了最終收斂，迭代步數(shù)較少。

圖9 不同測(cè)量誤差下橢圓污垢厚度的反演值與精確值Fig.9 Inversion and exact values of elliptical fouling for different measurement errors

圖10 橢圓函數(shù)邊界缺陷反演時(shí)的收斂歷程Fig.10 Convergence history of elliptic function boundary defect inversion

4 結(jié)論

(1)通過設(shè)定多種工況進(jìn)行數(shù)值模擬，驗(yàn)證了二維穩(wěn)態(tài)導(dǎo)熱反問題數(shù)學(xué)模型的有效性與精確性。研究結(jié)果顯示，該反演模型能夠?qū)崿F(xiàn)對(duì)內(nèi)壁面兩類缺陷的識(shí)別，且具有較高的識(shí)別精度。

(2)對(duì)圓管內(nèi)壁腐蝕減薄及污垢增厚兩類缺陷進(jìn)行了反演結(jié)果穩(wěn)定性的探究，計(jì)算結(jié)果表明，初值選取對(duì)反演結(jié)果無明顯影響，反演精度隨著外壁面測(cè)溫點(diǎn)數(shù)目的增加而提高，但精度提高有限。

(3)引入隨機(jī)測(cè)量誤差，探討了反演結(jié)果對(duì)測(cè)量誤差的敏感性。數(shù)值算例結(jié)果顯示，本文算法具有一定的抗噪性，當(dāng)存在σ=0.5 ℃的標(biāo)準(zhǔn)偏差時(shí)，仍能得到較為精確的結(jié)果。

目前，本文只研究了二維圓管內(nèi)壁面缺陷的穩(wěn)態(tài)識(shí)別，缺少基于具體實(shí)驗(yàn)的驗(yàn)證工作。下一步需要結(jié)合實(shí)際工況來優(yōu)化二維圓管穩(wěn)態(tài)導(dǎo)熱反演模型，并進(jìn)一步開發(fā)二維圓管瞬態(tài)導(dǎo)熱反演模型。