充分发展湍流管道内壁边界的红外定量识别*

张 林，杨 立，范春利

(海军工程大学 动力工程学院， 湖北 武汉 430033)

许多热力设备的管道内壁长期受管内环境的影响产生金属损失、腐蚀等缺陷，对设备的安全运行带来隐患[1]。因此，对管道内壁边界的检测和监测具有重要意义[2-4]。

红外无损检测技术是指通过红外热像仪得到设备表面温度的变化情况来判断设备内部信息的技术，由于其非接触、无损伤、检测速度快、扫描面积大等独特的优势已成为一种非常有效的无损检测技术[5-7]，被广泛应用于石化、冶金、电力等领域。

传热反问题理论是开展内部缺陷红外定量检测与识别的重要基础。基于导热反问题的红外检测研究已实现了对热物性参数、内热源强度、边界条件和边界形状等参数的精确检测和识别[8-11]。基于对流换热反问题的研究，由于流动的复杂性，不适定性严重，尚处于起步阶段。湍流由于质点相互掺混、碰撞的流动，不适定性更严重，解的存在性与唯一性更难确定，湍流对流换热反问题的研究比层流对流换热反问题的研究发展更慢[12]。

基于对流换热反问题的管道或平板通道的定量识别已由一些学者进行了研究和讨论[13-14]。 Chen等[15-16]使用共轭梯度法估计二维轴对称层流管道外壁的未知霜层边界形状和管道内壁上的未知污垢边界形状；Rouizi等[17]使用红外热像仪测量得到1 mm厚的平板微通道的外表面温度分布，基于测得的外表面温度分布对微通道内流体的温度分布进行了估计。

在工业生产中，大多数管道内流体为湍流流动状态。随着工业生产对设备安全性和可靠性的要求越来越高，发展基于对流换热反问题的湍流管道内壁边界的红外定量检测和识别具有非常重要的意义，但是关于湍流管道内壁边界的识别研究尚未见相关报道。

在实际中，内壁缺陷的存在往往会破坏管道几何的对称性，但是直接开展三维管道内壁边界识别研究，管道内流体变化复杂，管壁轴向和周向的传热互相影响，反问题研究的不适定性强，内壁边界识别规律难以分析归纳。参考Chen等[15-16]、肖丹等[18]的研究，首先开展二维轴对称管道的研究工作，为三维管道内部缺陷定量识别提供理论依据。

本文提出基于管道外壁温度分布，使用多物理场软件COMSOL建立充分发展的湍流管道模型，并与传热反问题算法——Levenberg-Marquardt算法相结合的识别方法，对二维轴对称充分发展的湍流管道内壁边界的识别问题进行研究。

1 对流换热正问题

构建如图1所示的分段充分发展的湍流管道模型。管道内部为远高于环境温度的热流体，管道置于自然环境中，通过管道外壁的自然对流和辐射换热将热量传递给周围环境。

对入口段和充分发展段分别进行建模，充分发展段会出现入出口效应，而且管壁轴向导热在分段处中断，因此设置虚拟管道前置入口段L1用以保证流体形成充分发展流，设置虚拟管道前置段L2和虚拟管道后置出口段L4来避免入出口效应和轴向导热的影响。本文研究的管道检测段为图中管道黑色区域L3。R1和R2分别为管道内壁和外壁的半径。入口段L1入口的流体均匀速度和均匀温度分别为uin和Tin，检测段L3入口的流体速度和温度分别为ufu,in和Tfu,in。在工作一段时间后，假定在管道检测段L3内壁上生成了一层腐蚀缺陷，缺陷处的不规则内壁边界形状为R(x)。

(a) 入口段(a) Entrance region

(b) 充分发展段(b) Fully developed region图1 管道的物理模型Fig.1 Schematic of pipeline

对管道的流动和传热特性进行数值模拟，作如下合理假设[19]：

1)管内流体稳态且不可压缩；

2)流体在边界上为无滑移流动；

3)忽略重力影响。

则管道模型内部流体流动的质量守恒方程、动量守恒方程和能量守恒方程为：

ρ·(u)=0

(1)

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

(2)

ρcpuT=·(kT)

(3)

其中，u是速度矢量，I是单位矩阵，ρ是流体密度，p是压力，μ是动力黏度，cp是定压比热，T是温度，k是导热系数。

管壁稳态导热的能量守恒方程为：

2T=0

(4)

对流换热正问题是在管道检测段L3的内壁边界形状R(x)和其他所有热物性参数、边界条件均为已知时，求解管道对流换热问题得到管道检测段L3外壁温度分布。

对于充分发展的湍流管道内壁形状识别问题，管道内壁在每次迭代过程中都会发生变化，因此管道内流场不能被指定，需要反复求解对流换热问题，计算过程中的模型自动重建、网格自动划分和求解计算是必需的。

多物理场软件COMSOL可以在模型计算过程中自动重建几何模型、更新网格和求解计算，同时COMSOL和MATLAB之间不需要使用输入/输出文件就可以进行数据传递[20-21]，可以用来计算许多实际而复杂的对流换热问题。因此，本文使用COMSOL对对流换热正问题进行求解，在MATLAB中使用读取函数mphinterp直接读取正问题的外壁温度分布。

2 反问题识别

对于反问题，不规则内壁边界R(x)被认为是未知的，而式(1)～(4)中其他所有热物性参数、边界条件均为已知的。此外，在管外壁上通过红外热像仪采集得到的温度分布被认为是可用的，将测量得到的温度用Y表示，Y=[Y1,Y2，…，Yi，…，Ym]T，i=1…m，其中m表示测温点数目。反问题的求解方法是以式(5)最小化的方式获得：

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

(5)

式中：T(R)为在迭代计算中，根据识别参数R求解对流换热正问题得到的管道检测段L3外壁温度分布。

2.1 Levenberg-Marquardt算法

为使式(5)取最小值，通过求导可得[22-23]：

(6)

(7)

式中，

(8)

Ωn=diag[(Jn)TJn]

(9)

n为迭代次数。

2.2 收敛条件

通过Levenberg-Marquardt算法迭代识别管道检测段L3内壁边界形状的收敛条件为：

S(R)<ε

(10)

式中，ε为一个很小的正数[24]。

2.3 MATLAB调用COMSOL

MATLAB调用COMSOL的语句为：

model=mphload(′E:pipeline_heat_thransfer′);

其中，“E:pipeline_heat_thransfer”为求解外壁温度分布的mph文件的所在路径，mphload表示MATLAB加载COMSOL的mph文件。

2.4 数值计算过程

通过Levenberg-Marquardt算法根据管道外壁温度分布定量识别管道检测段L3内壁边界形状R(x)的迭代求解步骤：

1)给出管道内壁边界形状R(x)的初始假设R0作为输入，开始迭代计算；

2)将估计值Rn传递到COMSOL中，MATLAB调用COMSOL，求解得到管道检测段L3的外壁温度分布T(Rn)；

3)根据式(10)判断是否收敛，如果收敛则输出结果，否则继续；

4)采用迭代式(7)修正识别参数Rn，求得Rn+1返回到步骤2。

3 仿真结果与讨论

数值实验中，建立二维轴对称管道对流换热模型，管道内半径R1为0.047 m，壁厚0.010 m。依据管道湍流充分发展经验公式L/D≥60，管道模型数据取为：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。

管内流质为空气，物性参数取自COMSOL材料库。管壁导热系数为17.6 W·m-1·K-1；外界环境温度为25 ℃，管道外壁发射率为0.95。使用COMSOL 5.2，采用标准k-ε湍流模型，进行正问题模型的求解计算。湍流模型收敛条件设置为1×10-6。

管道内壁缺陷只存在于管道检测段L3，虚拟管道段L1、L2和L4内壁均不存在缺陷。在管道检测段L3外壁表面上均匀布置的测温点数为11个，内壁边界均匀分布的离散点亦为11个。算例中，含内壁缺陷的管道内壁边界形状R(x)为：

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

(11)

3.1 不规则内壁边界的红外特征规律

对模型的流体域和管道壁固体域均进行结构化网格剖分。为保证模型的计算精度，对模型网格进行独立性检验。对径向(流体域+固体域)×(虚拟入口段L1、虚拟前置段L2、检测段L3和虚拟出口段L4)的网格划分分别为(12+5)×(141+12+100+7)、(25+5)×(282+25+200+15)和(50+5)×(564+50+400+30)的3套网格进行试算。在检测段入口速度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%。为同时保证计算精度和计算效率，选用第二套网格划分方案。

使用第二套网格划分方案，对管道内壁形状为式(11)的管道检测段L3进行网格划分，如图2所示。可以看出，管道检测段L3管壁中间部分的厚度小于两端部分的厚度，且从中间部分到两端部分内壁面呈圆滑过渡。

图2 管道检测段L3的网格划分Fig.2 Grid-meshing of pipeline inspection section L3

图3为管道检测段L3检测表面温度分布和管道检测段L3入口速度ufu,in、入口温度Tfu,in的关系。图中，检测表面最大温差(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 内壁发生变化前后，检测表面的温度分布和 入口速度ufu,in、入口温度Tfu,in的关系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

外壁温度分布中最高温度值和最低温度值之差；而内壁变化后管道检测段L3外壁温度分布和内壁没有发生变化时的温度分布之差为绝对温差(Absolute Temperature Difference, ATD)。

由图3可得，在内壁无缺陷时，管壁温度分布呈线性降低，且入口速度ufu,in越慢，入口温度Tfu,in越高，换热越充分，管壁温度分布最大温差越大；存在内壁缺陷时，管壁温度分布随内壁缺陷变化，最大温差略有变化，且入口速度ufu,in越慢，入口温度Tfu,in越高，管壁温度分布最大温差越大。综上，缺陷的存在与否，不改变外壁温度分布最大温差随入口速度ufu,in、入口温度Tfu,in的变化规律。

图4和图5分别为管道检测段L3检测表面温度分布绝对温差和管道检测段L3入口速度ufu,in、入口温度Tfu,in的关系。

(a) Tfu,in=373.15 K

(b) Tfu,in=473.15 K

(c) Tfu,in=573.15 K图4 检测表面绝对温差和入口速度ufu,in的关系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的关系Fig.5 Relationship between absolute temperature difference of outer pipeline surface and inlet temperature Tfu,in

从图4中可得：由缺陷引起的绝对温差在缺陷初始处基本为零，即有缺陷时与无缺陷时的外壁温度在缺陷初始处差别很小。入口速度ufu,in越快，绝对温差正增长的幅值越大。在入口速度ufu,in较慢时，由于检测表面温度受管壁轴向导热的影响，绝对温差在前半缺陷处呈下降趋势，为负增长；随着入口速度ufu,in的增大，绝对温差在前半缺陷处负增长的幅值逐渐变小直至负增长消失。综合绝对温差中正增长和负增长的变化规律，入口速度ufu,in越快，外壁温度的绝对温差越大。

从图5中可以看出：入口温度Tfu,in越高，绝对温差正增长的幅值越大；同时，由于检测表面温度受管壁轴向导热的影响，随着入口温度Tfu,in的升高，绝对温差在前半缺陷处逐渐形成负增长，且负增长的幅值越来越大。综合绝对温差中正增长和负增长的变化规律，入口温度Tfu,in越高，外壁温度绝对温差越大。

因此，入口速度ufu,in越快，入口温度Tfu,in越高，由缺陷引起的外壁温度绝对温差越大。结合图3中的分析可得，在含内壁缺陷的高温管道中，最大温差和由缺陷引起的绝对温差并不是同步增加的。

3.2 内壁边界识别

通过在不同入口速度ufu,in和入口温度Tfu,in下，求解对流换热正问题得到的管道检测段L3的外壁温度分布，来模拟在实际检测工作中测量得到的外壁温度分布Y，对式(11)所示的内壁缺陷进行识别研究。

要说明的是，在实际检测过程中，热像仪对管道外壁进行检测工作，得到的是二维热图，可以使用热像仪自带的软件进行温度提取工作，使用“线提取”得到管道外壁上沿轴向变化的温度分布。

图6和图7分别为在不同入口速度ufu,in和入口温度Tfu,in时，不同初始边界假设(初始边界假设分别为0.047 m和0.055 m)的内壁边界识别结果(ε=0.000 1)。

从图6可以看出，入口速度ufu,in越快，不同初值得到的识别结果曲线越贴合真实缺陷曲线，即识别结果越好，符合3.1节图4的分析，即入口速度ufu,in越快，温度分布绝对温差越来越大。由图7得，随着入口温度Tfu,in的升高，识别结果曲线先是贴近真实缺陷曲线，然后远离真实缺陷曲线，即识别结果先是变好，然后变差，结合3.1节图5(c)的分析，可得识别结果变好主要是由于温度分布绝对温差随着入口温度Tfu,in的升高而增大，识别结果变差主要是由于在入口温度Tfu,in继续升高的过程中温度分布绝对温差在前半缺陷处逐渐形成负增长。

因此，在对流换热反问题中，绝对温差大，识别结果未必好。在进行湍流管道的内壁边界识别时，应尽可能提高管内入口速度，提高绝对温差。同时，将入口温度尽可能地控制在一个合理的范围内：入口温度太低，绝对温差过小；入口温度过高，绝对温差前半部分出现的负增长会导致识别结果变差。

(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时内壁边界识别结果Fig.6 Identification results of inner boundary with inlet temperature Tfu,in=573.15 K

从图6和图7亦可以看出，在入口速度较快且入口温度没有导致绝对温差出现负增长时，不同的初值得到几乎相同的识别结果；而在绝对温差出现负增长时，不同的初值得到不同的识别结果，识别结果有好有坏，即初值对识别结果的影响大小取决于绝对温差中是否出现负增长。

同时还可得，无论哪种入口条件下的识别结果，前半缺陷的识别精度都很高，而后半缺陷的识别精度变低，在不规则内壁终点处识别精度略微变差。这是一个明显而且非常有用的经验。在不同的入口条件下识别时，尽管不能识别得到完整的内壁缺陷，但能得到前半缺陷的准确形状和缺陷的最大值等信息，这对于缺陷检测工作来说已经是很有用的了。

(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时内壁边界识别结果Fig.7 Identification results of inner boundary with inlet velocity ufu,in=25 m·s-1

4 结论

针对基于表面测温的充分发展湍流管道内壁边界形状的识别问题，使用COMSOL构建分段充分发展湍流管道模型，关联COMSOL和MATLAB，利用Levenberg-Marquardt算法对二维轴对称充分发展的湍流管道的内壁边界形状的识别进行了研究。算例验证了该方法识别充分发展湍流管道内壁边界形状的有效性，可以得出如下结论：

1)在含内壁缺陷的湍流管道中，入口速度越慢，且入口温度越高，外壁最大温差越大；入口速度越快，且入口温度越高，由缺陷引起的外壁绝对温差越大。因此，外壁最大温差和由缺陷引起的绝对温差并不是同步增加的。

2)对流换热反问题的识别规律和导热反问题的识别规律并不一致，在对流换热反问题中，由于绝对温差中负增长的出现，绝对温差大，识别结果未必好。

3)反问题方法的识别精度在含缺陷内壁的终点处略微变差，但能准确地预测得到内壁的其余轮廓，这对管道维护有很大的帮助。