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

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

引 言

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

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

當(dāng)前的大多數(shù)研究都只對管道內(nèi)壁腐蝕減薄缺陷進(jìn)行單一的討論，涉及管道內(nèi)壁污垢增厚缺陷識別的相對較少，因此本文提出基于有限元法和L- M算法構(gòu)建反問題數(shù)學(xué)模型，利用COMSOL與MATLAB聯(lián)合仿真穩(wěn)定、快速的優(yōu)勢，將管道內(nèi)徑作為反演參數(shù)，分別對管道內(nèi)壁腐蝕減薄缺陷與污垢增厚缺陷進(jìn)行定量識別，并構(gòu)造數(shù)值實驗對反演結(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)，外壁面對流換熱系數(shù)hout=10 W/(m2·K)，環(huán)境溫度Ta=25 ℃，內(nèi)壁面對流換熱系數(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)

由于所研究內(nèi)容為二維圓管內(nèi)壁面缺陷的穩(wěn)態(tài)識別，基于此作如下合理假設(shè)：1)管道忽略軸向?qū)岬挠绊懀?)管道導(dǎo)熱系數(shù)和管內(nèi)外的對流換熱系數(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)熱正問題是一個定解問題，而導(dǎo)熱反問題則是一個最優(yōu)化問題。

2.2.1L-M算法

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

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

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

(7)

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

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

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

(8)

式中，M為外壁面測點(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(θ)所求得的外壁面溫度的計算值，表示為

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

(10)

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

(11)

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

(12)

(13)

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

Ωk=diag[(Jk)TJk]

(14)

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

S(R)<ε

(15)

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

2.2.2數(shù)值計算過程

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

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

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

2.2.3反演識別誤差

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

Tmea=Texact+σω

(16)

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

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

(17)

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

(18)

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

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

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

3.1 管道減薄缺陷

對于管道減薄缺陷，在外壁面上均勻設(shè)置36個測點(diǎn)，內(nèi)壁面由36個均勻分布的離散節(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)驗證內(nèi)壁面幾何邊界呈階躍式規(guī)律變化下的反演結(jié)果，并探討不同初值R0對反演結(jié)果的影響。由圖3可以看出，當(dāng)R0取不同值時，反演曲線與精確曲線的走向基本趨于一致，靠近階躍邊界處的識別效果相對較差，但程序仍可較好地識別出管道內(nèi)壁的邊界形狀。由表1可以進(jìn)一步得出，3種初始邊界假設(shè)下，反演結(jié)果的平均相對誤差最大僅相差0.015%，說明初始邊界假設(shè)對內(nèi)壁面幾何邊界的反演幾乎沒有影響；此外，3種初值下得出計算結(jié)果的平均相對誤差均不超過0.1%，驗證了程序反演的精確性。

表1 Case1～Case3平均相對誤差值

3.1.2鍥形函數(shù)

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

(20)

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

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

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

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

3.2 管道增厚缺陷

同樣地，對于管道增厚缺陷，在外邊界上也均勻設(shè)置36個測點(diǎn)，含有污垢的內(nèi)壁面亦由36個均勻分布的離散節(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可以看出，不同初值下的反演值與精確值相差甚小，且反演曲線呈小范圍波動。由表2可得，3種初始值下反演的平均相對誤差均不足0.1%，說明初值對均勻污垢反演結(jié)果的影響較小，程序?qū)τ趦?nèi)壁面均勻污垢有較好的識別能力。

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

表2 Case4～Case6平均相對誤差值

3.2.2三角形函數(shù)

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

R(θ)=

(23)

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

圖8 不同測溫點(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)偏差時，反演值圍繞精確值附近上下波動，且隨著標(biāo)準(zhǔn)偏差的增大，波動的幅度和頻率也愈為明顯。根據(jù)表2可知，在標(biāo)準(zhǔn)偏差增大到0.5 ℃時，波動值仍可保持在1%以內(nèi)，說明反演模型在識別內(nèi)壁污垢時穩(wěn)定性較好。圖10為當(dāng)污垢以橢圓函數(shù)變化時的管道迭代收斂歷程圖。由圖可以看出，對橢圓污垢的反演與對正弦減薄缺陷的反演一樣，僅是經(jīng)歷3次迭代就達(dá)到了最終收斂，迭代步數(shù)較少。

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

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

4 結(jié)論

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

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

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

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