基于导热反问题圆管内壁腐蚀减薄及污垢增厚的识别

张经豪 熊 平 郝睿智 卢 涛

(北京化工大学 机电工程学院， 北京 100029)

引 言

在工业生产中，管道广泛应用于输气、供水和运油等领域。由于管道长期处于高温高压的工作状态，在使用过程中容易造成内壁的腐蚀、脱落，甚至产生污垢黏连附着在内壁上导致出现内壁增厚等现象，这些缺陷的存在成为安全生产和使用的潜在隐患[1]。因此，对管道内壁缺陷的检测具有重要意义[2-3]。

导热反问题 (inverse heat conduction problem，IHCP)已在很多工程领域中得到广泛应用[4]。导热反问题常用于计算反演边界条件[5-9]、热物性参数[10]、物体内部热源[11]及初始条件等未知项，也可用于物体内部缺陷[12-15]和内壁缺陷的识别。Chen等[16-17]采用共轭梯度法估计了管道外壁的未知霜层边界形状和管道内壁上的未知污垢边界形状，但只单一讨论了沿管道轴向变化的霜层或污垢分布情况；文献[18-19]对双层炉膛内壁的边界形状进行了反演，但没有探讨初始值及测点数目对反演结果精确度的影响；文献[20-24]应用共轭梯度法对管道内壁面的边界形状进行了稳态或瞬态的识别，并分别对所建反演模型的可行性进行了探究；张林等[25]在改进的一维修正算法基础上提出一维加权法，对二维圆管交界面形状进行了定量识别，提高了识别效率；Mohasseb 等[26]采用遗传算法作为非适定导热反问题的计算方法，求解了方形板的热流密度边界，但算法的稳定性与收敛速度有待提高；文献[3]、[27-30]均采用列文伯格- 马夸尔特(Levenberg- Marquardt，L- M)算法对未知参数进行了计算，算例验证了算法的有效性和优越性。在众多求解导热反问题的数值计算方法中，属于梯度类算法的L- M算法因具有稳定性好、构造思路简单、收敛快等优势而被广泛应用。

当前的大多数研究都只对管道内壁腐蚀减薄缺陷进行单一的讨论，涉及管道内壁污垢增厚缺陷识别的相对较少，因此本文提出基于有限元法和L- M算法构建反问题数学模型，利用COMSOL与MATLAB联合仿真稳定、快速的优势，将管道内径作为反演参数，分别对管道内壁腐蚀减薄缺陷与污垢增厚缺陷进行定量识别，并构造数值实验对反演结果的精确性和稳定性进行分析。

1 二维圆管物理模型

二维圆管物理模型如图1所示。管道外径a=0.3 m，壁厚d=25 mm。管道的导热系数λ1=17.60 W/(m·K)，污垢的导热系数λ2=3.14 W/(m·K)，外壁面对流换热系数hout=10 W/(m2·K)，环境温度Ta=25 ℃，内壁面对流换热系数hin=1 000 W/(m2·K)，管内流体温度Tf=200 ℃，θ为极角，r为出现缺陷后管道的内壁面极径。

图1 管道两类缺陷的物理模型Fig.1 Physical model of two types of defects in pipelines

2 IHCP数学模型

2.1 导热正问题(DHCP)

由于所研究内容为二维圆管内壁面缺陷的稳态识别，基于此作如下合理假设：1)管道忽略轴向导热的影响；2)管道导热系数和管内外的对流换热系数均为常数。

2.1.1管道增厚缺陷

圆管的导热偏微分方程为

(1)

假设给定内壁面和外壁面均为第三类边界条件，则边界条件为

(2)

(3)

污垢增厚缺陷存在着交界面Γ，在此处

T1=T2

(4)

(5)

式中，下标1代表管壁区域Ω1，下标2代表污垢区域Ω2；下标in、out分别为圆管内壁面和外壁面；n代表法线方向；下标int代表圆管与污垢交界面边界。

2.1.2管道减薄缺陷

管道减薄缺陷只存在于管壁区域Ω1，将式(3)进行修改可得式(6)。

(6)

式(6)与式(1)、式(2)共同构成了管道减薄缺陷导热正问题的数学描写。

2.2 导热反问题

导热正问题是一个定解问题，而导热反问题则是一个最优化问题。

2.2.1L-M算法

L- M算法是非线性回归中回归参数最小二乘估计的一种方法，L- M算法将最速下降法和线性化方法相结合，从而较快地找到最优值。

对于两类内壁面缺陷的反演，内壁面几何边界都是未知的，可由极坐标向量R(θ)描述

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

(7)

式中，r为极半径；N为内壁面离散节点个数。

数学偏微分方程组(式(1)～(5)以及式(1)、(2)、(6))需要获得外壁面测点温度以满足其封闭性，离散后的外壁面温度测量值记为

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

(8)

式中，M为外壁面测点数。

进一步构建该最优化问题的目标函数为

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

(9)

式中，T(R)为根据R(θ)所求得的外壁面温度的计算值，表示为

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

(10)

通过式(9)对未知参量R(θ)求偏导得到目标函数梯度的最小化条件为

(11)

式中，J(R)为敏度系数矩阵，表示为

(12)

(13)

式中，k为迭代次数；Ωk为一对角矩阵，表示为

Ωk=diag[(Jk)TJk]

(14)

根据偏差原理给出收敛条件为

S(R)<ε

(15)

式中，ε在不考虑测量误差时为一小正数。

2.2.2数值计算过程

MATLAB使用mphload函数加载COMSOL的mph文件，通过mphinterp函数便可读取数据，无需使用输入/输出文件。COMSOL在被调用时可以自动构建几何模型、更新网格和求解计算，有效提高了反演效率。

利用L- M算法求解导热反问题的流程图如图2所示。

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

2.2.3反演识别误差

为验证反演结果对引入测量误差的敏感性，在由导热正问题计算得到的外壁面精确值的基础上引入随机误差，来模拟实际的测量值。

Tmea=Texact+σω

(16)

式中，Tmea为外壁面温度测量值；Texact为外壁面温度精确值；σ为标准偏差；ω为区间[-2.576,2.576]内服从标准正态分布的随机数，该区间能达到99%的测量可靠度。

当考虑实际测温误差时，收敛标准ε写成

(17)

为了表征反演值与精确值的偏离程度，定义平均相对误差为

(18)

式中，rn,exact为内径精确值；rn,est为内径反演值。

3 计算结果与分析

通过运行程序进行数值计算，以验证管道减薄缺陷与增厚缺陷导热反问题反演的精确性。首先设定内壁面缺陷形状，由导热正问题得出外壁面温度模拟的实际测量值，为导热反问题提供输入数据；所设定的内壁面缺陷形状参数作为校验数据。为验证程序的稳定性和抗噪性，分别探讨初值、测温点数目及测量误差对反演结果的影响。

3.1 管道减薄缺陷

对于管道减薄缺陷，在外壁面上均匀设置36个测点，内壁面由36个均匀分布的离散节点拟合得到。设定沿周向变化的内壁面几何形状为阶跃式函数、锲形函数和正弦函数，分别探讨这3种不同工况下的导热反问题。

3.1.1阶跃式函数

设定内壁面边界形状按阶跃式变化，即Case 1。

(19)

图3 不同初值下阶跃式几何边界的反演值与精确值Fig.3 Inversion and exact values for a stepped geometrical boundary with different initial values

根据式(19)验证内壁面几何边界呈阶跃式规律变化下的反演结果，并探讨不同初值R0对反演结果的影响。由图3可以看出，当R0取不同值时，反演曲线与精确曲线的走向基本趋于一致，靠近阶跃边界处的识别效果相对较差，但程序仍可较好地识别出管道内壁的边界形状。由表1可以进一步得出，3种初始边界假设下，反演结果的平均相对误差最大仅相差0.015%，说明初始边界假设对内壁面几何边界的反演几乎没有影响；此外，3种初值下得出计算结果的平均相对误差均不超过0.1%，验证了程序反演的精确性。

表1 Case1～Case3平均相对误差值

3.1.2锲形函数

设定内壁面几何边界形状以锲形函数变化，即Case 2。

(20)

在式(20)所示的内壁边界形状以锲形函数变化的条件下，研究不同外壁面温度测点数目对识别结果的影响。由图4可以看出，当M=24时，程序对内壁边界的识别效果相对较差，尤其是在阶跃点附近更为明显；随着外壁面温度测点数目的增加，阶跃点处的识别效果逐渐变好，当M=72时，反演曲线与精确曲线高度吻合。结合表1可进一步得出，外壁面测温点数越多，识别效果就越好。M=36时的平均相对误差为0.058%，M=72时的平均相对误差为0.011%，均可较好地识别出内壁边界形状。

图4 锲形几何边界不同测温点数时的反演值与精确值Fig.4 Inversion and exact values for a wedge-shaped geometric boundary with different numbers of measurement points

3.1.3正弦函数

设定内壁边界形状随正弦函数变化，即Case 3。

R(θ)=0.285+0.01sinθ

(21)

内壁面边界形状按式(21)以正弦规律变化，引入标准偏差，分析不同偏差对反演结果的影响。如图5所示，可以看出当σ=0 ℃时，反问题得到的反演曲线与正问题得出的精确曲线几乎重合；当σ=0.2 ℃和σ=0.5 ℃时，两条反演曲线均在精确值附近上下小范围波动，且在σ=0.5 ℃时变化得相对更为明显。根据表1不难进一步得出，随着标准偏差的增大，反演误差也随之增大，但3种标准偏差下的平均相对误差均不超过1.5%，说明采用此方法反演边界形状对测量误差的变化不敏感，具有一定的抗噪性。图6是当边界形状为正弦函数时的迭代收敛历程图，由图可看出，迭代步长先大后小，且只需要迭代3次即可满足所设定的收敛条件。对比图6(c)与图6(d)可知，运用L- M算法可在低迭代步数下准确地识别出缺陷。

图5 正弦几何边界不同测量误差时的反演值与精确值Fig.5 Inversion and exact values of a sine geometric boundary with different measurement errors

图6 正弦函数边界缺陷反演时的收敛历程Fig.6 Convergence history of sine function boundary defect inversion

3.2 管道增厚缺陷

同样地，对于管道增厚缺陷，在外边界上也均匀设置36个测点，含有污垢的内壁面亦由36个均匀分布的离散节点拟合得到，可通过反演管道内径进而得出污垢的厚度。设定的内壁面污垢几何形状包括均匀增厚缺陷、沿周向变化的三角形函数以及椭圆函数，分别探讨这3种不同工况下的导热反问题。

3.2.1均匀增厚

设定内壁面有均匀厚度的污垢，即Case 4。

R(θ)=0.265

(22)

如式(22)所示，假设周向含有均匀污垢的圆管内径为0.265 m，即污垢的厚度为0.01 m，分别设置不同的初值R0，所得污垢厚度Δd反演结果如图7所示。由图7可以看出，不同初值下的反演值与精确值相差甚小，且反演曲线呈小范围波动。由表2可得，3种初始值下反演的平均相对误差均不足0.1%，说明初值对均匀污垢反演结果的影响较小，程序对于内壁面均匀污垢有较好的识别能力。

图7 不同初值下均匀污垢厚度的反演值与精确值Fig.7 Inversion and exact values of uniform fouling with different initial values

表2 Case4～Case6平均相对误差值

3.2.2三角形函数

设定内壁几何边界形状按三角形函数变化，即Case 5。

R(θ)=

(23)

如图8所示，含有污垢的管道内壁面半径(即内壁污垢形状)随周向呈三角形规律变化。根据式(23)对反演结果进行验证，同时探讨改变外壁面测温点个数对反演结果的影响。由表2可知，增加外壁面测温点个数可以更加准确地识别出内壁面的污垢形状，但提高的程度有限，当测温点的个数达到一定值时，增加节点个数对识别精度的提高影响并不大。另外，由表2结果也可看出，当测温点个数较少时识别精度仍然较高，M=24时的平均相对误差为0.033%，表明仍能较准确地识别出内壁面的污垢形状。

图8 不同测温点数下三角形污垢厚度的 反演值与精确值Fig.8 Inversion and exact values of triangle fouling with different numbers of temperature measurement points

3.2.3椭圆函数

设定内壁面几何边界形状以椭圆函数变化，即Case 6。

(24)

式(24)为内壁边界形状随椭圆函数变化的情况，即内壁污垢形状亦呈椭圆函数变化。图9为不同标准偏差下污垢尺寸的反演值与精确值。由图可看出，无偏差下反演值与精确值二者的曲线走向基本同步；当存在标准偏差时，反演值围绕精确值附近上下波动，且随着标准偏差的增大，波动的幅度和频率也愈为明显。根据表2可知，在标准偏差增大到0.5 ℃时，波动值仍可保持在1%以内，说明反演模型在识别内壁污垢时稳定性较好。图10为当污垢以椭圆函数变化时的管道迭代收敛历程图。由图可以看出，对椭圆污垢的反演与对正弦减薄缺陷的反演一样，仅是经历3次迭代就达到了最终收敛，迭代步数较少。

图9 不同测量误差下椭圆污垢厚度的反演值与精确值Fig.9 Inversion and exact values of elliptical fouling for different measurement errors

图10 椭圆函数边界缺陷反演时的收敛历程Fig.10 Convergence history of elliptic function boundary defect inversion

4 结论

(1)通过设定多种工况进行数值模拟，验证了二维稳态导热反问题数学模型的有效性与精确性。研究结果显示，该反演模型能够实现对内壁面两类缺陷的识别，且具有较高的识别精度。

(2)对圆管内壁腐蚀减薄及污垢增厚两类缺陷进行了反演结果稳定性的探究，计算结果表明，初值选取对反演结果无明显影响，反演精度随着外壁面测温点数目的增加而提高，但精度提高有限。

(3)引入随机测量误差，探讨了反演结果对测量误差的敏感性。数值算例结果显示，本文算法具有一定的抗噪性，当存在σ=0.5 ℃的标准偏差时，仍能得到较为精确的结果。

目前，本文只研究了二维圆管内壁面缺陷的稳态识别，缺少基于具体实验的验证工作。下一步需要结合实际工况来优化二维圆管稳态导热反演模型，并进一步开发二维圆管瞬态导热反演模型。