CN109033591A - 基于逆瞬态模型的城市非金属管道泄漏定位方法 - Google Patents

Abstract

本发明提供一种基于逆瞬态模型的城市非金属管道泄漏定位方法，使用Fluent软件模拟非金属管道泄漏和无泄漏状态下的压力、流量、流速等参数分布，找出两种状态下各参数随时间和空间变化的规律；在实验室进行非金属埋地管道泄漏模拟试验，分别在泄漏和未泄漏状态下测量管道的流量、流速以及多处的压力参数数据；根据模拟软件得到的各参数变化规律，建立非金属管道泄漏状态的逆瞬态模型并对实验数据进行分析；应用全局和局部算法，确定收敛性最好的目标参数J_k和τ_k，得出瞬态求解器计算出的测压头；计算出的测压头与实验所得测压头代入逆瞬态模型，以确定泄漏大小和位置。与现有技术相比，该方法只需测量管道压力、流量参数即可，易实现。

Description

基于逆瞬态模型的城市非金属管道泄漏定位方法

技术领域

本发明涉及城市非金属管道泄漏检测技术领域，特别是涉及一种基于逆瞬态模型的城市非金属管道泄漏定位方法。

背景技术

非金属管材具有强度高、密度小、耐腐蚀性强、绝缘性优、使用寿命长等优点，在石化行业、建筑业、市政工程、农业等领域需求不断加大，促使非金属管道得到快速发展。但同时管道由于质量，第三方破坏等问题发生泄漏，造成人员伤亡和财产损失。管道的安全形势十分严峻，泄漏、火灾、爆炸等因素也日益为公众所关注。为了保障城市管道的安全运行，开发管道泄漏检测技术十分必要。

发明内容

本发明所要解决的技术问题是：为了克服现有技术中的不足，本发明提供一种基于逆瞬态模型的城市非金属管道泄漏定位方法，通过管道模拟和实地检测相结合，以实现管道泄漏精确定位。

本发明解决其技术问题所要采用的技术方案是：一种基于逆瞬态模型的城市非金属管道泄漏定位方法，该种方法主要用于检测复杂的存在弯折，众多阀门、仪表和管支的城市管道中较大流量压力下的泄漏，包括以下步骤：

S1：用AnsysIcem对管道进行建模，使用Fluent软件模拟非金属管道泄漏和无泄漏状态下的压力、流量、流速等参数分布，分别找出非金属管道泄漏和无泄漏状态下的压力、流量、流速等参数随时间和空间变化的规律；

本发明利用ICEM建模，用Fluent16.0模拟，在保证模拟精度达标且准确可靠的前提下，依据实际情况对管道泄漏的物理模型加以简化。模拟中假定埋地非金属管长24m，并且距离管道进口6m处存在一个泄漏孔，泄漏孔为圆形，大小为1.5mm，位于管道正上方，管道壁厚7.5mm，管道材质为聚乙烯(PE)。

通过ICEM建模对管道网格划分，对于整体管道网格，选择性加密了两侧管壁和泄漏孔，泄漏孔处加密，为保证网格质量，整体管道共有437002个非结构网格，且质量均在85.5％以上，网格数量和质量满足计算精度要求。

S2：在实验室进行非金属埋地管道泄漏模拟试验，分别在泄漏和未泄漏两种状态下测量管道流量、温度及多处的压力等参数数据，与步骤S1中获得的模拟结果对比和验证，突出泄漏状态下的规律，获得实验管道测量到的测压头；

S3：根据步骤S2中模拟软件模拟得到的压力、流量、流速等变化规律，建立非金属管道泄漏状态的逆瞬态模型；

逆瞬态模型，为实现计算出的测压头和实验测量到的测压头之间的加权均方误差最小化，具体如下

其中，OF为目标函数(用于实现泄漏定位)；DT为观测到的水力瞬变事件的时间步数(优化中使用的第一个时间实例，t＝1，表示瞬态事件的开始)；H为通过瞬态求解器计算出的测压头；H^*为实验设施中测量到的测压头；n^H为系统中压力观测点的数量；H_t,j为不同测量节点和泄漏时间下的通过瞬态求解器计算出的测压头，H*t,j为不同测量节点和泄漏时间下的实验设施中测量到的测压头；j为测量节点；t为泄漏后的时间，取t＝1,2,3…10。

当前适用于非金属管道可行的泄漏检测方法中，内部传感器法检测和定位效果较好，超声、视觉检测等方法检测可行性并不确定，涡流，漏磁等方法并不能使用，在此基础上，提出逆瞬态法实现泄漏检测与定位；

步骤4的处理是对计算出的测压头H进行求解，具体为应用全局和局部算法，得出Kelvin-Voigt元素的数目以及收敛性最好的目标函数值J_k和τ_k，为确定逆瞬态模型中通过瞬态求解器计算出的测压头做准备，计算出的测压头H的具体求解步骤包括；

S4.1：在matlab环境下，使用有效的全局搜索方法PSO(粒子群算法)，适应度线性定标，随机余数采样(SRS)作为选择方案，每个PSO评估使用九十代，找到Kelvin-Voigt元素的数目；

与当前主流的遗传算法相比较，PSO算法的优势在于简单，容易实现，无需梯度信息，参数少，其天然的实数编码特点特别适合于处理实优化问题；Kelvin-Voigt元素，凯尔文-沃伊特元素，可以描述应力和应变的弹性和粘性成分。

S4.2：分别投放1，10，20，50，100，200，500个随机种子，观察目标函数OF的函数值收敛情况，得出全局搜索下收敛性最佳时的J_k和τ_k值，其中Z＝{J_k，τ_k}，k＝1到N_KV，其中，N_KV为Kelvin-Voigt元素的数目；

S4.3：使用局部搜索(Levenberg-Marquardt，简称LM)算法来估计和优化系数J_k(Z＝{J_k，τ_k}；k＝1到N_KV)，并考虑由PSO确定的区域中τ_k的值，计算出三种常用的K值和使用LM方法的目标函数值，并绘折线图，挑选收敛性最好的目标函数值J_k和τ_k；

全局算法得到Kelvin-Voigt元素的数目以及目标函数值J_k，τ_k；局部算法进一步优化参数J_k，τ_k。

LM算法是一种迭代技术，它定位一个函数的局部最小值，该最小值表示为几个非线性实值函数的平方和，一旦初始值足够精确，LM会展现出快速收敛。

S4.4：将Kelvin-Voigt元素的数目和收敛性最好的目标函数值J_k和τ_k值代入蠕变函数方程及MOC特征方程，从而得出H值；

其中，蠕变函数方程可以通过不同的方法获得，

(1)蠕变函数方程为：

J(t)＝e_r(t)/σ₀ (2)

所述蠕变函数方程用以描述材料的延迟应变e_r(t)与载荷σ₀之间的比率，为解决特征方程中无法直接得到的延迟应变时间导数而引入。

(2)蠕变函数方程又可以通过参数模拟反向计算确定，

其中，J₀是管道的弹性模量的倒数，J_k是k元素弹簧的弹性模量的倒数，τ_k＝μ_k/E_k，μ_k是非金属管道的粘度，N_KV是Kelvin-Voigt元素的数量，参数J_k，τ_k和Kelvin-Voigt元素数目通过S4确定。

而MOC特征方程由联列封闭管道中一维瞬态流动的动量和连续性方程得到，具体为：

其中，t为时间；H为测压头；Q为流量；g为重力引起的加速度；A为管的横截面积；e_r为延迟应变分量；h_f为单位长度的损失(湍流条件下的其中f为达西-威斯巴赫摩擦系数和D为管内径)；a为水击波速，通常情况下远大于流体流速，故近似看作±a＝Δx/Δt，使用矩形计算网格并忽略对流项(与弹性波速度相比，流体速度可忽略不计)。

S5：将步骤S4中计算出的测压头与步骤S2中实验测量到的测压头代入逆瞬态模型，得到目标函数值OF，用于实现泄漏定位；

S6：在实验所得多处的OF值中，挑选最小值的位置节点，并在最小值位置节点前后布置等距的且距离稍小于前一次(可选择前一次的一半距离)，个数与其相同的多个测量节点，重复步骤计算各节点H值，并得到OF，当各节点距离小于管道检测长度的4％时终止重复步骤；

2006年，大连理工大学李俊花在其博士论文中应用逆瞬态法和压力梯度法迭代求解长输管道泄漏点位置，利用遗传算法优选泄漏点前后的摩阻系数，与压力梯度法结合求解泄漏点的位置，通过输水非金属管道的泄漏实验得到这种方法漏点定位的平均误差达4％左右。因此本发明力求误差在4％以内。

本发明的有益效果是：本发明提供的一种基于逆瞬态模型的城市非金属管道泄漏定位方法，与现有技术相比，该方法只需测量管道压力、流量参数即可，易实现。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的流程示意图。

图2是本发明管道系统的结构示意图。

图3是实验管道示意图。

图4是具有特定时间间隔的特征网格示意图。

图5是0.18MPa压力泄漏时管道上下游压力变化图。

图6是0.10MPa压力泄漏时管道上下游压力变化图。

图中：1、管道系统，2、流量计，3、压力传感器，4、管道内物质流向。

具体实施方式

现在结合示例对本发明作详细的说明。此图为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。

如图1所示，一种基于逆瞬态模型的城市非金属管道泄漏定位方法，该种方法主要用于检测复杂的存在弯折，众多阀门、仪表和支管的城市管道中较大流量压力下的泄漏，包括以下步骤：

S1：用AnsysIcem对管道系统1进行建模，使用Fluent软件模拟。

如图2所示，实验管道系统1采用多段管道连接形成的具有分支和弯折的管道结构，在管道系统1上安装有多个压力传感器3，并且在管道系统1上游部分和下游部分安装流量计2，用于测量稳态流量，流量计2可以采用电磁流量计或体积流量计，图中箭头表示管道内物质流向4。

利用ICEM建模，用Fluent16.0模拟，在保证模拟精度达标且准确可靠的前提下，依据实际情况对管道泄漏的物理模型加以简化，简化后管道模型参考图3所示。

通过ICEM建模对管道网格划分，对于整体管道网格，选择性加密了两侧管壁和泄漏孔，泄漏孔处加密，为保证网格质量，整体管道共有437002个非结构网格，且质量均在85.5％以上，网格数量和质量满足计算精度要求。

Fluent模拟边界条件设置的流场进出口边界条件下如下：管道进口条件类型Velocity-inlet，取值5m/s；管道泄漏点条件类型Pressure-outlet，取值为0；管道出口条件类型Pressure-vent，分别取值0.10Mpa，0.18Mpa，0.26Mpa。

S2：在实验室进行非金属埋地管道泄漏模拟试验，分别在泄漏和未泄漏两种状态下测量管道流量、温度及多处的压力等参数数据，与步骤S1中获得的模拟结果对比、验证，突出泄漏状态下的规律；

由模拟与实验得出以下规律：

当液体管道距离较短时，未泄漏状态下，首末端压力差值很小；出现泄漏之后，管道各处压力都会呈近似线性梯度逐渐减小，并趋于一个稳定值；线性梯度在管道压力较小时(＜0.3Mpa)，影响因素主要是泄漏孔大小；管道沿程损失对管道质量流率的影响非常小；存在泄漏时，相较于流量变化，压力变化更加明显，压力变化是管道各段都会逐渐减小，流量变化是泄漏孔处流率骤减。

S3：根据步骤S1中模拟软件模拟得到的压力、流量、流速等变化规律，建立非金属管道泄漏状态的逆瞬态模型；

逆瞬态模型，为实现计算出的测压头和实验测量到的测压头之间的加权均方误差最小化，具体如下

其中，OF＝目标函数(用于实现泄漏定位)；DT＝观测到的水力瞬变事件的时间步数(优化中使用的第一个时间实例，t＝1，表示瞬态事件的开始)；H＝通过瞬态求解器计算出的测压头；H^*＝实验设施中测量到的测压头；n^H系统中压力观测点的数量；H_t,j为不同测量节点和泄漏时间下的通过瞬态求解器计算出的测压头，H^* _t,j为不同测量节点和泄漏时间下的实验设施中测量到的测压头；j为测量节点(实施例中共上游、下游两个节点则j＝1,2)；t为泄漏后的时间，取t＝1,2,3…10。

步骤4的处理是对计算出的测压头H进行求解，具体为应用全局和局部算法，得出Kelvin-Voigt元素的数目以及收敛性最好的目标函数值J_k和τ_k，为确定逆瞬态模型中通过瞬态求解器计算出的测压头做准备，计算出的测压头H的具体求解步骤包括；

S4.1：在matlab环境下，使用有效的全局搜索方法PSO(粒子群算法)，适应度线性定标，随机余数采样(SRS)作为选择方案，每个PSO评估使用九十代，找到Kelvin-Voigt元素的数目；

与当前主流的遗传算法相比较，PSO算法的优势在于简单，容易实现，无需梯度信息，参数少，其天然的实数编码特点特别适合于处理实优化问题；Kelvin-Voigt元素，凯尔文-沃伊特元素，可以描述应力和应变的弹性和粘性成分。

S4.2：分别投放1，10，20，50，100，200，500个随机种子，观察目标函数OF的函数值收敛情况，得出全局搜索下收敛性最佳时的J_k和τ_k值，其中Z＝{J_k，τ_k}，k＝1到N_KV；

S4.3：使用局部搜索(LM)算法来估计和优化系数J_k(Z＝{J_k，τ_k}；k＝1到N_KV)，并考虑由PSO确定的区域中τ_k的值，计算出三种常用的K值和使用LM方法的目标函数值，并绘折线图，挑选收敛性最好的目标函数值J_k和τ_k；

LM(Levenberg-Marquardt)算法是一种迭代技术，它定位一个函数的局部最小值，该最小值表示为几个非线性实值函数的平方和，一旦初始值足够精确，LM会展现出快速收敛。

S4.4：将Kelvin-Voigt元素的数目和收敛性最好的目标函数值J_k和τ_k值代入蠕变函数方程及MOC特征方程，得出H值；

其中，蠕变函数方程可以通过不同的方法获得，

(1)蠕变函数方程为：

J(t)＝e_r(t)/σ₀ (2)

所述蠕变函数方程用以描述材料的延迟应变e_r(t)与载荷σ₀之间的比率，为解决特征方程中无法直接得到的延迟应变时间导数而引入。

(2)蠕变函数方程又可以通过参数模拟反向计算确定，

其中，J₀是管道的弹性模量的倒数，J_k是k元素弹簧的弹性模量的倒数，τ_k＝μ_k/E_k，μ_k是非金属管道的粘度，N_KV是Kelvin-Voigt元素的数量，参数J_k，τ_k和Kelvin-Voigt元素数目通过S4确定；

而MOC特征方程由联连列封闭管道中一维瞬态流动的动量和连续性方程得到，具体为：

其中，t为时间；H为测压头；Q为流量；g为重力引起的加速度；A为管的横截面积；e_r为延迟应变分量；h_f为单位长度的损失(湍流条件下的其中f为达西-威斯巴赫摩擦系数和D为管内径)；a为水击波速，通常情况下远大于流体流速，故近似看作±a＝Δx/Δt，使用矩形计算网格并忽略对流项(与弹性波速度相比，流体速度可忽略不计)；

S5：将步骤S4中计算出的测压头与步骤S2中实验测量到的测压头代入逆瞬态模型，得到目标函数值OF，用于实现泄漏定位；

S6：在实验所得多处的OF值中，挑选最小值的位置节点，并在最小值位置节点前后布置等距的且距离稍小于前一次(可选择前一次的一半距离)，个数与其相同的多个测量节点，重复步骤计算各节点H值，并得到OF，当各节点距离小于管道检测长度的4％时终止重复步骤；

2006年，大连理工大学李俊花在博士论文中应用逆瞬态法和压力梯度法迭代求解长输管道泄漏点位置，利用遗传算法优选泄漏点前后的摩阻系数，与压力梯度法结合求解泄漏点的位置。通过输水非金属管道的泄漏实验得到这种方法漏点定位的平均误差达4％左右。因此本发明力求误差在4％以内。

根据上述方法进行实验，实验管道总长34m，其中PE管长24m，管道类型为埋地管；传输介质为水，传感器分别布置在被测管道上下游，并且距离PE管道进口6m处存在一个泄漏孔，泄漏孔为圆形，泄漏孔径为1.5mm，位于管道正上方，管道壁厚7.5mm，实验管道布置如图3所示。

实验所需流量、压力、温度等数据经管道上各传感器测量后由WS-5921/U系列的USB接口将数据传输至电脑中。

使用Vib’SYS振动信号采集、处理和分析软件中自带绘图工具进行绘图，得到图5和图6，其中横坐标为时间(30s)，纵坐标为测量到的电信号(vlot)，为清晰显示压力变化图，选择起始纵坐标为1.1，结束纵坐标为1.5；其中，图5为0.18Mpa压力下，6m处1.5mm的泄漏孔泄漏时，压力随时间变化曲线，其中，A为管道上游压力变化曲线，B为管道下游压力变化曲线；图6为0.10Mpa压力下，6m处1.5mm的泄漏孔泄漏时，压力随时间变化曲线，其中，A为管道上游压力变化曲线，B为管道下游压力变化曲线。

从测量仪软件内导出数据至matlab2017，计算上下游泄漏后30s内的参数均值，经传感器量程标定后得到表1和表2。

表1管道上游各参数值

表2管道下游各参数值

数据代入模型之前，需要确定计算得测压头。具体如下：经PSO算法模拟，得到最佳Kelvin-Voigt元素的数目为7。进一步数值模拟得到J_K，τ_k见表3所示。

表3模拟所得J_K和τ_k值

可以看出随机种子数量增加时，参数J_K，τ_k的值逐渐减小，种子数量达200以上时，变化梯度趋于平缓，故，计算测压管水头时，取J_K＝0.03000，τ_k＝0.0530。

此后，应用局部搜索(LM)算法优化系数J_k(Z＝{J_k，τ_k}；k＝1到7)，计算得到收敛性最好的目标函数值J_k和τ_k，最后代入S4.4得到计算得测压头H，通过与实验测量得测压头的计算得到最终结果，见表4。

表4上下游计算结果

由表4数据结果可知，管道内压力较大时，OF值较小，说明该方法能更好地模拟流量较大的管道内流体的流动情况；此外，相同压力下，上游OF值小于下游OF值，可知，相较于下游传感器，泄漏孔的位置更加靠近上游传感器。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关的工作人员完全可以在不偏离本发明的范围内，进行多样的变更以及修改。本项发明的技术范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

Claims (4)

1.一种基于逆瞬态模型的城市非金属管道泄漏定位方法，其特征在于：包括如下步骤：

S1：用AnsysIcem对管道进行建模，使用Fluent软件模拟非金属管道泄漏和无泄漏状态下的压力、流量和流速参数分布，分别找出非金属管道泄漏和无泄漏状态下的压力、流量和流速参数随时间和空间变化的规律；

S2：在实验室进行非金属埋地管道泄漏模拟试验，分别在泄漏和未泄漏两种状态下测量管道流量、温度及多处的压力参数数据，与步骤S1中获得的模拟结果对比和验证，突出泄漏状态下的规律，获得实验管道测量到的测压头；

S3：根据步骤S1得到的压力、流量和流速的变化规律，建立非金属管道泄漏状态的逆瞬态模型；

为实现计算出的测压头和实验测量到的测压头之间的加权均方误差最小化，逆瞬态模型具体如下：

其中，OF为目标函数；DT为观测到的水力瞬变事件的时间步数；H为通过瞬态求解器计算出的测压头；H^*为实验管道测量到的测压头；n^H为系统中压力观测点的数量；H_t,j为不同测量节点和泄漏时间下的通过瞬态求解器计算出的测压头，H^* _t,j不同测量节点和泄漏时间下的实验设施中测量到的测压头；j为测量节点；t为泄漏后的时间；

S4：应用全局和局部算法，得出Kelvin-Voigt元素的数目以及收敛性最好的目标函数值J_k和τ_k，确定逆瞬态模型中通过瞬态求解器计算出的测压头H，具体步骤包括：

S4.1：在matlab环境下，使用有效的全局搜索方法PSO(粒子群算法)，适应度线性定标，随机余数采样(SRS)作为选择方案，每个PSO评估使用九十代，找到Kelvin-Voigt元素的数目；

S4.2：分别投放1，10，20，50，100，200，500个随机种子，观察目标函数OF的函数值收敛情况，得出全局搜索下收敛性最佳时的J_k和τ_k值，其中Z＝{J_k，τ_k}，k＝1到N_KV，其中，N_KV为Kelvin-Voigt元素的数目；

S4.3：使用局部搜索(LM)算法来估计和优化系数J_k(Z＝{J_k，τ_k}；k＝1到N_KV)，并考虑由PSO确定的区域中τ_k的值，计算出三种常用的K值和使用LM方法的目标函数值，并绘折线图，挑选收敛性最好的目标函数值J_k和τ_k；

S4.4：将Kelvin-Voigt元素的数目和收敛性最好的目标函数值J_k和τ_k值代入蠕变函数方程及MOC特征方程，从而得出H值；

S5：将实验测量到的测压头与计算出的测压头共同代入逆瞬态模型，得到目标函数值OF；

S6：在实验所得多处的OF值中，挑选最小值的位置节点，并在最小值位置节点前后布置等距的且距离小于前一次，个数与其相同的多个测量节点，重复步骤计算各节点H值，并得到OF，当各节点距离小于管道检测长度的4％时终止重复步骤。

2.如权利要求1所述的基于逆瞬态模型的城市非金属管道泄漏定位方法，其特征在于：步骤4.4所述的蠕变函数方程具体为：

J(t)＝e_r(t)/σ₀

所述蠕变函数方程用以描述材料的延迟应变e_r(t)与载荷σ₀之间的比率，为解决特征方程中无法直接得到的延迟应变时间导数而引入。

3.如权利要求1所述的基于逆瞬态模型的城市非金属管道泄漏定位方法，其特征在于：步骤4.4所述的蠕变函数方程通过参数模拟反向计算确定，具体为：

其中，J₀是管道的弹性模量的倒数，J_k是k元素弹簧的弹性模量的倒数，τ_k＝μ_k/E_k，μ_k是非金属管道的粘度，N_KV＝Kelvin-Voigt元素的数量，参数J_k，τ_k和Kelvin-Voigt元素数目通过步骤S4确定。

4.如权利要求1所述的基于逆瞬态模型的城市非金属管道泄漏定位方法，其特征在于：步骤4.4所述的MOC特征方程由联列封闭管道中一维瞬态流动的动量和连续性方程得到，具体为：

其中，t为时间；H为测压头；Q为流量；g为重力引起的加速度；A＝管的横截面积；e_r为延迟应变分量；h_f为单位长度的损失(湍流条件下的其中f为达西-威斯巴赫摩擦系数和D为管内径)；a为水击波速，通常情况下远大于流体流速，故近似看作±a＝Δx/Δt，使用矩形计算网格并忽略对流项。