引言

在工业生产中，管道广泛应用于输气、供水和运油等领域。由于管道长期处于高温高压的工作状态，在使用过程中容易造成内壁的腐蚀、脱落，甚至产生污垢黏连附着在内壁上导致出现内壁增厚等现象，这些缺陷的存在成为安全生产和使用的潜在隐患^[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。管道的导热系数λ₁=17.60 W/(m·K)，污垢的导热系数λ₂=3.14 W/(m·K)，外壁面对流换热系数h_out=10 W/(m²·K)，环境温度T_a=25 ℃，内壁面对流换热系数h_in=1 000 W/(m²·K)，管内流体温度T_f=200 ℃，θ为极角，r为出现缺陷后管道的内壁面极径。

2 IHCP数学模型 2.1 导热正问题(DHCP)

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

2.1.1 管道增厚缺陷

圆管的导热偏微分方程为

$ \frac{\partial^{2} T}{\partial x^{2}}+\frac{\partial^{2} T}{\partial y^{2}}=0 $

(1)

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

$ -\lambda_{1}\left(\frac{\partial T_{1}}{\partial n_{\text {out }}}\right)=h_{\text {out }}\left(T_{\text {out }}-T_{\mathrm{a}}\right) $

(2)

$ -\lambda_{2}\left(\frac{\partial T_{2}}{\partial n_{\text {in }}}\right)=h_{\text {in }}\left(T_{\text {in }}-T_{\mathrm{f}}\right) $

(3)

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

$ T_{1}=T_{2} $

(4)

$ -\lambda_{1}\left(\frac{\partial T_{1}}{\partial n_{\mathrm{int}, 1}}\right)=-\lambda_{2}\left(\frac{\partial T_{2}}{\partial n_{\mathrm{int}, 2}}\right) $

(5)

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

2.1.2 管道减薄缺陷

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

$ -\lambda_{1}\left(\frac{\partial T_{1}}{\partial n_{\mathrm{in}}}\right)=h_{\mathrm{in}}\left(T_{\mathrm{in}}-T_{\mathrm{f}}\right) $

(6)

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

2.2 导热反问题

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

2.2.1 L-M算法

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

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

$ \boldsymbol{R}(\theta)=\left[r_{1}, r_{2}, \cdots, r_{n}, \cdots, r_{N}\right]^{\mathrm{T}} $

(7)

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

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

$ \boldsymbol{Y}=\left[Y_{1}, Y_{2}, \cdots, Y_{m}, \cdots, Y_{M}\right]^{\mathrm{T}} $

(8)

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

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

$ S(\boldsymbol{R})=[\boldsymbol{Y}-\boldsymbol{T}(\boldsymbol{R})]^{\mathrm{T}}[\boldsymbol{Y}-\boldsymbol{T}(\boldsymbol{R})] $

(9)

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

$ \boldsymbol{T}(\boldsymbol{R})=\left[T_{1}, T_{2}, \cdots, T_{m}, \cdots, T_{M}\right]^{\mathrm{T}} $

(10)

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

$ \nabla S(\boldsymbol{R})=-2 \boldsymbol{J}^{\mathrm{T}}(\boldsymbol{R})[\boldsymbol{Y}-\boldsymbol{T}(\boldsymbol{R})]=0 $

(11)

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

$ \boldsymbol{J}(\boldsymbol{R})=\left[\frac{\partial \boldsymbol{T}^{\mathrm{T}}(\boldsymbol{R})}{\partial \boldsymbol{R}}\right]^{\mathrm{T}} $

(12)

对式(11)的T(R)进行泰勒级数展开舍去高阶项，并加入衰减参数μ_LM^k来调节迭代速度，得到迭代式

$ \begin{array}{l} \ \ \ \ \ \ \ \ \boldsymbol{R}^{k+1}=\boldsymbol{R}^{k}+\left[\left(\boldsymbol{J}^{k}\right)^{\mathrm{T}} \boldsymbol{J}^{k}+\mu_{\mathrm{LM}}^{k} \boldsymbol{\varOmega}^{k}\right]^{-1}\left(\boldsymbol{J}^{k}\right)^{\mathrm{T}}[\boldsymbol{Y}- \\ \left.\boldsymbol{T}\left(\boldsymbol{R}^{k}\right)\right] \end{array} $

(13)

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

$ \boldsymbol{\varOmega}^{k}=\operatorname{diag}\left[\left(\boldsymbol{J}^{k}\right)^{\mathrm{T}} \boldsymbol{J}^{k}\right] $

(14)

迭代开始采用最速下降法以较大的μ_LM^k进行计算，后来降低μ_LM^k通过高斯牛顿法得到最终反演值。

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

$ S(\boldsymbol{R})<\varepsilon $

(15)

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

2.2.2 数值计算过程

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

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

2.2.3 反演识别误差

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

$ T_{\text {mea }}=T_{\text {exact }}+\sigma \omega $

(16)

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

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

$ \boldsymbol{\varepsilon}=\sum\limits_{m=1}^{M} \sigma_{m}^{2}=M \boldsymbol{\sigma}^{2} $

(17)

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

$ E_{\mathrm{RR}}=\frac{1}{N} \sum\limits_{n=1}^{N}\left|\frac{r_{n, \text { exact }}-r_{n, \text { est }}}{r_{n, \text { exact }}}\right| \times 100 \% $

(18)

式中，r_{n, exact}为内径精确值；r_{n, est}为内径反演值。

3 计算结果与分析

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

3.1 管道减薄缺陷

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

3.1.1 阶跃式函数

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

$ R(\theta)= \begin{cases}0.275, & \theta \in[0,5 {\rm{ \mathsf{ π} }} / 6) \cup(7 {\rm{ \mathsf{ π} }} / 6,2 {\rm{ \mathsf{ π} }}] \\ 0.287\ 5, & \theta \in[5 {\rm{ \mathsf{ π} }} / 6,7 {\rm{ \mathsf{ π} }} / 6]\end{cases} $

(19)

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

3.1.2 锲形函数

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

$ R(\theta)= \begin{cases}0.275, & \theta \in[0, {\rm{ \mathsf{ π} }}) \cup({\rm{ \mathsf{ π} }}, 2 {\rm{ \mathsf{ π} }}] \\ 0.287\ 5, & \theta={\rm{ \mathsf{ π} }}\end{cases} $

(20)

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

3.1.3 正弦函数

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

$ R(\theta)=0.285+0.01 \sin \theta $

(21)

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

3.2 管道增厚缺陷

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

3.2.1 均匀增厚

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

$ R(\theta)=0.265 $

(22)

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

3.2.2 三角形函数

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

$ \begin{array}{l} \ \ \ \ \ \ \ \ \ R(\theta)= \\ \begin{cases}0.005 \theta+0.25, & \theta \in[0, {\rm{ \mathsf{ π} }})\\ -0.005(\theta-{\rm{ \mathsf{ π} }})+(0.005 {\rm{ \mathsf{ π} }}+0.25), & \theta \in[{\rm{ \mathsf{ π} }}, 2 {\rm{ \mathsf{ π} }}]\end{cases} \end{array} $

(23)

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

3.2.3 椭圆函数

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

$ R(\theta)=\sqrt{(0.27 \cos \theta)^{2}+(0.26 \sin \theta)^{2}} $

(24)

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

4 结论

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

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

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

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