CN110500511B - 一种城市非金属管道泄漏定位方法 - Google Patents

一种城市非金属管道泄漏定位方法 Download PDF

Info

Publication number
CN110500511B
CN110500511B CN201910742549.0A CN201910742549A CN110500511B CN 110500511 B CN110500511 B CN 110500511B CN 201910742549 A CN201910742549 A CN 201910742549A CN 110500511 B CN110500511 B CN 110500511B
Authority
CN
China
Prior art keywords
leakage
pipeline
flow
pressure
state
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910742549.0A
Other languages
English (en)
Other versions
CN110500511A (zh
Inventor
马逸飞
郝永梅
蒋军成
邢志祥
杨克
倪磊
吴洁
朱一龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changzhou University
Original Assignee
Changzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changzhou University filed Critical Changzhou University
Priority to CN201910742549.0A priority Critical patent/CN110500511B/zh
Publication of CN110500511A publication Critical patent/CN110500511A/zh
Priority to US17/254,285 priority patent/US11592351B2/en
Priority to PCT/CN2020/108372 priority patent/WO2021027803A1/zh
Application granted granted Critical
Publication of CN110500511B publication Critical patent/CN110500511B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F17—STORING OR DISTRIBUTING GASES OR LIQUIDS
    • F17D—PIPE-LINE SYSTEMS; PIPE-LINES
    • F17D5/00—Protection or supervision of installations
    • F17D5/02—Preventing, monitoring, or locating loss
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Energy or water supply
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M3/00—Investigating fluid-tightness of structures
    • G01M3/02—Investigating fluid-tightness of structures by using fluid or vacuum
    • G01M3/26—Investigating fluid-tightness of structures by using fluid or vacuum by measuring rate of loss or gain of fluid, e.g. by pressure-responsive devices, by flow detectors
    • G01M3/28—Investigating fluid-tightness of structures by using fluid or vacuum by measuring rate of loss or gain of fluid, e.g. by pressure-responsive devices, by flow detectors for pipes, cables or tubes; for pipe joints or seals; for valves ; for welds
    • G01M3/2807—Investigating fluid-tightness of structures by using fluid or vacuum by measuring rate of loss or gain of fluid, e.g. by pressure-responsive devices, by flow detectors for pipes, cables or tubes; for pipe joints or seals; for valves ; for welds for pipes
    • G01M3/2815—Investigating fluid-tightness of structures by using fluid or vacuum by measuring rate of loss or gain of fluid, e.g. by pressure-responsive devices, by flow detectors for pipes, cables or tubes; for pipe joints or seals; for valves ; for welds for pipes using pressure measurements
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Economics (AREA)
  • Evolutionary Computation (AREA)
  • Computing Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • General Business, Economics & Management (AREA)
  • Marketing (AREA)
  • Human Resources & Organizations (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Strategic Management (AREA)
  • Tourism & Hospitality (AREA)
  • Water Supply & Treatment (AREA)
  • Artificial Intelligence (AREA)
  • Operations Research (AREA)
  • Fluid Mechanics (AREA)
  • Mechanical Engineering (AREA)
  • Computational Linguistics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Examining Or Testing Airtightness (AREA)
  • Pipeline Systems (AREA)

Abstract

本发明提供一种城市非金属管道泄漏定位方法，首先通过数值模拟规律分析或基于马尔科夫链的流量分析法判断管道是否发生泄漏；对于泄漏管道建立非金属管道气体泄漏的逆瞬态控制方程，通过试验取得各个测量点不同时段的压力、流量数据代入控制方程，对实验数据进行分析；定义具有最小二乘准则的目标函数的非线性规划问题并应用序列二次规划法，将目标函数最小化，以确定泄漏大小和位置。与现有技术相比，该方法可以进行管道泄漏检测同时对有效泄漏面积进行检测，通过对测量点压力、流量的数据的分析可以准确地确定管道的泄漏参数及泄漏定位。

Description

一种城市非金属管道泄漏定位方法
技术领域
本发明涉及城市非金属管道泄漏检测技术领域，特别是涉及一种城市非金属管道泄漏定位方法。
背景技术
因非金属管材具有强度高、密度小、腐蚀性强、绝缘性优、使用寿命长等优点，其逐步取代金属管道用于能源、原材料等的管道运输。随着城市埋地管网越来越密集，伴随着这些管道数量的增多和运行时间的增长，非金属管道因接头渗漏、拉断，或遭第三方破坏等原因导致的管道泄漏隐患或事故日益增多。因此对管道故障的有效检测和及时处理显得尤为重要。近些年，国内外在管道泄漏检测技术以及泄漏定位等方面都做了大量的研究工作。其中基于逆瞬态法的管道泄漏检测方法研究发展尤为迅速。
逆向瞬态分析(ITA，inverse-transient analysis)，最初由国外的Pudar和Liggett在1992年提出，这种方法通常用于配水网络，是一种用于加压管道中摩擦系数泄漏检测和校准的强大方法，通过这种方法，启动瞬态流动并在系统某处测量压力，通过减小测量压力与实际压力的误差以达到校准泄漏参数，从而进行泄漏检测定位。为此目的，数值模型的建立以及各种优化技术的选择成为了ITA技术应用于管道泄漏检测的主要挑战。国外学者针对配水管网提出了液压瞬态求解器以提高模型精度减小误差，针对ITA中的优化求解采取了许多算法，2001年Vitkovsky采用的混合复杂演化(SCE)针对大孔泄漏(泄漏孔直径与管道直径的比值大于0.1)优化计算精度较好，却不适于小孔泄漏(泄漏孔直径与管道直径的比值等于小于0.1)，2010年B.S.Jung采用的粒子群算法(PSO)和遗传算法(GA)全局搜索能力较强，但计算量大，对目标函数的收敛性较差，对摩擦系数的校准能力较低。
发明内容
本发明所要解决的技术问题是：为了克服现有技术中的不足，本发明提供一种城市非金属管道泄漏定位方法。本发明通过数值模拟管道泄漏的特征或马尔科夫法对管道是否产生泄漏作出判断标准；针对输气管道采用了气体在管道中的连续性、运动方程从而描述气体在等温管道中瞬变流动，采用序列二次规划(SQP，sequence quadratic program)法来解决非线性约束问题从而实现最小化目标函数的目的，SQP方法具有很好的兼容性和快速收敛性，计算效率高，能够对小孔泄漏的泄漏参数精确校准，从而达到对输气管道泄漏的精确定位。
本发明解决其技术问题所要采用的技术方案是：一种城市非金属管道泄漏定位方法，该种方法主要用于检测复杂的存在弯折，众多阀门、仪表和管支的城市管道中较大流量压力下的泄漏，整个过程主要包括泄漏判断和泄漏定位两个步骤，其中，泄漏判断是指通过软件数值模拟规律分析或马尔科夫法两种方法判断管道泄漏；泄漏定位是指通过逆瞬态的分析方法将假定泄漏点的计算压力与测量压力差值最小化获得符合现实情况的泄漏参数来判别泄漏位置；只有准确判断出哪个管道存在泄漏，才能对泄漏位置进行定位，通过双重的筛选提高了检测效率。
具体包括以下步骤：
S1：判断管道是否发生泄漏；
在对管道泄漏状态进行筛选判断时，可以采用多种方法，本发明给出两种技术方案。
技术方案一：采用Fluent软件模拟规律分析进行泄漏判断，具体包括以下步骤：
用AnsysIcem对管道进行建模，使用Fluent软件模拟气体在非金属管道泄漏和无泄漏状态下的压力、流量参数分布，并通过实验进行验证，分别找出气体在非金属管道泄漏和无泄漏状态下的压力、流量参数随时间和空间变化的规律。对管道是否发生泄漏给出判断标准。
步骤A1，对实验管道进行建模，采用结构化的四面体网格划分区域，并对泄漏口附近进行适当的网格加密处理；
步骤A2，控制方程和算法选择：
数值模型选用标准k-ε双方程模型，在模型中，用湍能k反映特征速度，用湍能耗散率ε反映特征长度尺度，利用Boussinesq假定简化而形成了k方程和ε方程，联合上述两个方程一道组成封闭的方程组对流动进行模拟，具体为：
k方程：
Figure BDA0002164450830000031
ε方程：
Figure BDA0002164450830000032
其中：xi，xj，xk为坐标分量；u为速度矢量；ρ为流体密度；μ为粘性系数；μt为湍动粘度；Gk、Gb为速度梯度、浮力引发的湍动能k的产生项；各常数C1ε＝1.44，C2ε＝1.92，C3ε＝0.09，σk＝1.0，σt＝1.3；Sk、Sε为用户定义源项；YM为可压湍流中脉动扩张项；
模型的数值计算算法选用SIMPLE算法，对于给定的压力，求解离散形式的动量方程，得出速度场，然后对给定的压力加以修正，反复迭代，直到获得收敛的解；
步骤A3，边界条件设置，包括对进口和出口的条件设置、流体和管道材质的物性参数设置；
步骤A4，在给定的多组进口压力下，分别对管道泄漏和非泄漏状态的管道内压力及流量参数的变化进行模拟，找出压力和流量参数的变化规律，对管道是否发生泄漏给出判断标准，然后将实际测得的管道数据与该标准进行对比，确定泄漏管道和无泄漏管道；
步骤A4.1，对数值模拟泄漏和非泄漏状态的管道内压力、流量流场进行分析，找出泄漏前后整个管道内压力、流量的变化规律；
步骤A4.2，对实验管道泄漏和非泄漏状态下的压力、流量流场进行分析，找出泄漏前后整个管道内压力、流量的变化规律；
步骤A4.3，将数值模拟和管道实验所得到的泄漏前后管道内压力、流量的变化规律进行对比分析，总结得到判断管道产生泄漏的具体标准。
技术方案二：采用基于马尔科夫链的方法实现非金属管道泄漏的判断，具体包括以下步骤：
步骤B1，利用安装在管线上的流量传感器收集第一组X个流量变化率q1，q2，q3，q4，q5...qx，其中qi为流量的变化率＝ΔQn/Qn×100％，qi(i＝1，2，3，4，5...X)中，i为采集到的X个流量变化率的排序，ΔQn为传感器传输数据的时间点Tn的流量与时间点Tn-1的流量的变化量，Qn为时间段T内的实时流量，Tn为流量计收集传输流量数据的时间，n＝1，2，3，4，5…，然后设置四个变化状态，分别为第一状态：qi的值为0，表示管道为非泄漏状态；第二状态：qi的值为0-1(％)，表示管道为小型泄漏状态；第三状态：qi的值为1-3(％)，表示管道为泄漏扩大状态；第四状态：qi的值为3-100(％)，表示管道为大泄漏状态。
设定条件：管道的泄漏为突变过程，泄漏状态的变化为渐变构成，管道的工作状态由于其他各类因素的影响可能会直接从运行状态变为泄漏状态中某一程度的泄漏，且管道的泄漏状态不可逆，但是泄漏程度可逆，根据以上条件，得到管道泄漏的邻接矩阵A为：
Figure BDA0002164450830000051
其中，当Aij＝1表示第i状态与第j状态之间存在转移关系，当Aij＝0表示第i状态与第j状态之间不存在转移关系，(i＝1，2，3，4；j＝1，2，3，4)。
步骤B2，进行第一组流量变化率的状态分类，将各流量变化率分别归为这四个状态分类，以时间为序，按时间顺序统计在X个参数中，相邻两个参数选取的时间间隔，流量变化率的状态转移情况，形成状态转移概率矩阵P中的各个元素；
步骤B2.1，得到处于第一状态的流量变化率数量为n1，由第一状态至第一状态转移的数量为n11，由第一状态至第二状态转移的数量为n12，由第一状态至第三状态转移的数量为n13，由第一状态至第四状态转移的数量为n14，统计这n1个流量变化率在下一时刻的状态转移情况：
第一状态至第一状态的转移概率
Figure BDA0002164450830000052
第一状态至第二状态的转移概率
Figure BDA0002164450830000053
第一状态至第三状态的转移概率
Figure BDA0002164450830000054
第一状态至第四状态的转移概率
Figure BDA0002164450830000055
步骤B2.2，得到处于第二状态的流量变化率数量为n2，由第二状态至第二状态转移的数量为n22，由第二状态至第三状态转移的数量为n23，由第二状态至第四状态转移的数量为n24，统计这n2个流量变化率在下一时刻的状态转移情况：
第二状态至第一状态的转移概率p21＝0；
第二状态至第二状态的转移概率
Figure BDA0002164450830000056
第二状态至第三状态的转移概率
Figure BDA0002164450830000061
第二状态至第四状态的转移概率
Figure BDA0002164450830000062
步骤B2.3，得到处于第三状态的流量变化率数量为n3，由第三状态至第二状态转移的数量为n32，由第三状态至第三状态转移的数量为n33，由第三状态至第四状态转移的数量为n34，统计这n3个流量变化率在下一时刻的状态转移情况：
第三状态至第一状态的转移概率p31＝0；
第三状态至第二状态的转移概率
Figure BDA0002164450830000063
第三状态至第三状态的转移概率
Figure BDA0002164450830000064
第三状态至第四状态的转移概率
Figure BDA0002164450830000065
步骤B2.4，得到处于第四状态的流量变化率数量为n4，由第四状态至第二状态转移的数量为n42，由第四状态至第三状态转移的数量为n43，由第四状态至第四状态转移的数量为n44，统计这n4个流量变化率在下一时刻的状态转移情况：
第四状态至第一状态的转移概率p41＝0；
第四状态至第二状态的转移概率
Figure BDA0002164450830000066
第四状态至第三状态的转移概率
Figure BDA0002164450830000067
第四状态至第四状态的转移概率
Figure BDA0002164450830000068
则该马尔科夫链中的状态转移概率矩阵P为：
Figure BDA0002164450830000069
且满足：n1+n2+n3+n4＝X。
步骤B3，按照步骤B1的方法收集第二组X个流量变化率，并按照步骤B2的方法统计这X个流量变化率在每个变化状态内的数量，将此时实际流量变化率在每个区间的数量编为第一参数向量，设为
Figure BDA0002164450830000071
其中，α1表示第一参数向量，
Figure BDA0002164450830000072
表示在第一参数向量下，X个流量变化率中，在第一状态下的流量变化率的数量；
Figure BDA0002164450830000073
表示在第一参数向量下，X个流量变化率中，在第二状态下的流量变化率的数量；
Figure BDA0002164450830000074
表示在第一参数向量下，X个流量变化率中，在第三状态下的流量变化率的数量；
Figure BDA0002164450830000075
表示在第一参数向量下，X个流量变化率中，在第四状态下的流量变化率的数量；
由马尔科夫链计算，得到将来时间段的基于马尔可夫链的预测流量变化率的参数向量β1，其公式为：
Figure BDA0002164450830000076
由矩阵乘法：
Figure BDA0002164450830000077
Figure BDA0002164450830000078
Figure BDA0002164450830000079
Figure BDA00021644508300000710
其中，
Figure BDA00021644508300000711
表示在预测得出的在第一状态下的流量变化率的数量，
Figure BDA00021644508300000712
表示在预测得出的在第二状态下的流量变化率的数量，
Figure BDA00021644508300000713
表示在预测得出的在第三状态下的流量变化率的数量，
Figure BDA00021644508300000714
表示在预测得出的在第四状态下的流量变化率的数量；
在这X个流量参数中，以S表示为预测状态序列号，(S＝1，2，3，4)，选取含有流量参数呈现数量最多的状态，即选取
Figure BDA00021644508300000715
认作此时预测管道状态将呈现为第Smax状态。
步骤B4，按照步骤A1的方法获取收集第三组X个流量变化率的值，并按照步骤B2的方法统计其在各个变化状态的匹配数量情况，设为
Figure BDA0002164450830000081
Figure BDA0002164450830000082
在这X个流量参数中，以W表示为实际管道状态序列号(W＝1，2，3，4)，选取含有流量参数呈现数量最多的状态，即选取
Figure BDA0002164450830000083
认作管道实际处于第Wmax状态。
步骤B5，比较预测状态
Figure BDA0002164450830000084
和实际状态
Figure BDA0002164450830000085
中状态序列号S和W的值：
①若S＝W＝1，则判断管道处于正常运行状态；
②若S＜W，则管道处于开始泄漏前期阶段，且判断其为第W状态；
③若S＞W，则判断管道处于开始泄漏后期阶段，且判断其状态为第W状态；
④若S＝W≠1，则判断其状态为泄漏稳定阶段，且判断其状态为第W状态。
步骤B6，以步骤B3所得数据为第一组数据，步骤B4所得数据成为第二组数据，继续预测判定管道实时状态或泄漏情况，开始重新迭代数据。步骤B1-步骤B6的泄漏检测方法可作为利用次声波泄漏定位的辅助手段，当利用此方法检测到泄漏状态时。
筛选出存在泄漏的管道后，进入泄漏位置的定位步骤，以确定具体的泄漏位置，便于维护和维修。泄漏位置的定位具体包括以下步骤：
S2：管道泄漏实验：获取管道压力、流量数据
具体方式为通过逐渐关闭阀门产生瞬态管流，阀门关闭时间持续足够长，以最小化不稳定和摩擦系数不确定性的影响，得到各个测量点不同时段的压力、流量数据。
S3：建立管道泄漏控制方程判断泄漏位置
利用气体的连续性和运动方程求解出的非线性方程与泄漏孔泄漏量方程建立非金属管道气体泄漏状态的控制方程，将S2步骤测得监测点两端Δx处的压力PA、PB，流量MA、MB带入控制方程，得到关于监测点计算压力PL与流量M1、M2的方程，以此为依据，根据逆瞬态计算原理，将某一时间点监测点a计算压力PL与试验真实数据进行比较，并采用SQP算法对两者差值进行收敛达到最小值从而获得符合现实情况的泄漏参数；
其中，气体连续性和运动方程为：
Figure BDA0002164450830000091
Figure BDA0002164450830000092
其中：M为气体流速，kg·s-1；A为管道横截面积，m2；P为压力，Pa；ρ为气体密度，kg·m-3；D为管内径，m；λ为摩阻系数。
采用特征线法对公式(1)和(2)进行求解，得到两组特征线方程，用C+，C-进行表示，具体如：
Figure BDA0002164450830000093
Figure BDA0002164450830000094
Figure BDA0002164450830000095
Figure BDA0002164450830000096
其中，a为压力波速，m/s，气体管道等温流动时，a为定值。
公式(11)和(13)分别为C+，C-特征线AP，BP；公式(10)和(12)为满足各自特征线的相容性方程。将管道均分为N等份，步长为Δx，时间步长为Δt＝aΔx，分别沿着C+特征线AP和C-特征线BP对公式(10)和(12)进行积分，其中对相容性方程的第三项即摩阻项采用一阶近似，得到C+特征线和C-特征线上的两个非线性方程，具体如：
Figure BDA0002164450830000101
Figure BDA0002164450830000102
管道泄漏视为小孔泄漏，则燃气泄漏公式取决于燃气在泄漏口处的流动速度，且燃气在泄漏口处的流动速度一般为亚音速流动，具体如：
Figure BDA0002164450830000103
其中：ML为泄漏孔泄漏流量，kg·s-1；Ae＝CA0为有效泄漏面积，m2；C为孔口系数，跟泄漏孔的形状有关；A0为泄漏孔面积，m2；Pd为管道起点压力，Pa；k为气体的绝热系数，无量纲；R为气体常数，J·kg-1·K-1；T为气体的温度，K。
假设监测节点为泄漏点，建立泄漏控制方程对监测节点的压力、流量数据进行分析，以泄漏孔处为边界处，流入泄漏孔处之前Δx距离内的流动特性满足公式(14)，流出泄漏孔处之后Δx距离内的流动特性满足公式(15)。设流经泄漏孔处之前的流体参数下标为1，流经泄漏孔处之后的流体参数下标为2。则建立的泄漏控制方程具体如：
Figure BDA0002164450830000104
Figure BDA0002164450830000105
Figure BDA0002164450830000106
M1-M2＝ML (20)
P1＝P2＝PL (21)
其中，P1、M1为流入泄漏点之前的一段距离的流体压力，Pa、流量，kg·s-1；P2、M2为流出泄漏点之后的一段距离的流体压力，Pa、流量，kg·s-1；PL为泄漏点处的压力，Pa。
S4：定义目标函数
将监测点处的压力计算值与测量值之差的最小值作为控制目标，控制方程和边界条件作为限制条件，对确定泄漏参数进行反问题分析，定义具有最小二乘目标准则的目标函数，某个时间点距离监测点a两端Δx处测量点A、B的压力PA、PB和流量MA、MB由试验给出，单位时间步长Δt后监测点a两端的计算压力P1、P2(P1＝P2＝PL)和流量M1、M2是未知的，流量之差可以判断监测点a是否产生泄漏，根据S3泄漏控制方程将已知的数据带入并且表示出监测点a两端流量M1、M2关于压力PL的方程，通过将计算压力PL与真实压力的差值收敛到最小，得到符合真实情况的监测点流量数据，并将流量数据代入燃气泄漏量方程求出有效泄漏面积，从而确定监测点泄漏情况；
S4.1：定义最小二乘准则目标函数：
Figure BDA0002164450830000111
其中：E是目标函数；M为时间步数；Pi是计算出的压力，Pa；P′i是测得的压力，Pa。通过将目标函数E最小化为零，从而产生最符合实际情况的有效泄漏面积，并根据有效面积的值确定该节点是否泄漏。
S4.2：确定泄漏有效面积的限制范围：
0≤Aei≤Aemax (23)
其中：Aei为节点i处的有效泄漏面积；Aemax为泄漏面积的最大限制，确定为管道横截面积的合理比例。
S4.3：摩阻系数根据实验数据由摩阻计算公式计算给出，而输气管道的流动一般处于紊流区，具体为：
伯拉修斯公式：
Figure BDA0002164450830000121
其中，雷诺系数计算公式采用经验公式：
Figure BDA0002164450830000122
Figure BDA0002164450830000123
式中，入为摩阻系数，Re为雷诺系数，Qv为体积流量m3·s-1，D为管道内径，m，υ为运动粘度，m2·s-1，μ为动力粘度，Pa·s，ρ为气体密度，kg·m-3。
S5：算法优化
混合复杂演化(SCE)针对大孔泄漏优化计算精度较好，却不适于小孔泄漏，粒子群算法(PSO)和遗传算法(GA)全局搜索能力较强，但计算量大，对目标函数的收敛性较差，对摩擦系数的校准能力较低。SQP算法针对求解非线性优化问题具有收敛性好、计算效率高、边界搜索能力强等优势，且可以解决小孔泄漏的泄漏参数校准精确问题。为了解决目标函数的非线性规划问题，本专利采取SQP方法在S4设定监测点两端流量限制范围内根据S3所获得压力、流量数据将目标函数最小化确定有效泄漏面积。
S5.1：建立M文件，根据实验数据定义目标函数E(x)；
S5.2：建立M文件，定义约束条件中的非线性约束Ax≤b或Aeq·X＝Beq；
S5.3：确定迭代的初始值X0；
S5.4：确定变量的上下限VLB，VUB；
S5.5：建立主程序，非线性规划求解的函数为fmincon，运行求解。
S6：泄漏定位
S6.1：在S5最小化结束时，如果发现节点的泄漏面积为非零，则将其视为泄漏；
S6.2：由于逆瞬态法的定位方法主要是判断监测点处是否泄漏，所以定位精度与监测点的设置数量有密切关系，并不能给出传统方法的定位公式，所以如果泄漏点并不在监测点处，即两个相邻的监测点的泄漏面积都不为零，则认为泄漏点在两个相邻且泄漏面积都不为零的监测点之间，并在两个监测点之间布置等距的多个测量节点，将测得的数据和计算的数据代入目标函数重复S5，直至各节点距离小于管道检测长度的4％时终止重复步骤。
本发明的有益效果是：本发明提供的一种城市非金属管道泄漏定位方法，与现有技术相比，该方法在对非金属管道定位的同时可以确定有效泄漏面积，通过定义最小二乘准则目标函数，将反问题应用于泄漏检测，采用序列二次规划(SQP)方法确定泄漏参数，保证了目标函数的收敛速度和计算精度，从而对每个节点的泄漏有效面积的检测更加准确，快速。
附图说明
下面结合附图和实施例对本发明作进一步说明。
图1是本发明的流程示意图。
图2是本发明管道系统的结构示意图。
图3是具有特定时间间隔的特征网格示意图。
图4是泄漏网格划分示意图。
图5是两种状态下模拟管道压力分布图。
图6是两种状态下模拟管道速度矢量分布图。
图7是管道实验数据采集图。
图中：1-上游流量传感器，2-上有压力传感器，3-上游次声波传感器，4-泄漏阀，5-下游次声波传感器，6-下游压力传感器，7-下游流量传感器。
具体实施方式
现在结合附图对本发明作详细的说明。此图为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。
如图2所示，本发明的实验管段采用三段直管段依次连接形成U形结构，其中，X1-X2管段长3.2m，X2-X3长1.6m，X3-X4长3.2m，泄漏孔在X5处，距离管道进口2m处，泄漏孔径为1mm。其中，X1位于管道进口处，X4位于管道出口处，X6距离进口处2.8m处，且在管道进口端依次设置上游流量传感器1、上有压力传感器2和上游次声波传感器3，在管道出口端依次设置下游流量传感器7、下有压力传感器6和下游次声波传感器5，泄漏阀4设置在X3处。
实施例一：本实施例中采用Fluent软件模拟规律分析进行泄漏判断。
如图1所示，本发明的一种城市非金属管道泄漏定位方法，包括使用Fluent软件模拟非金属管道泄漏和非泄漏状态下的压力、流量参数分布，研究气体在非金属管道发生泄漏的特征，给出判断管道是否发生泄漏的判断标准，并将实验数据与标准相比较判断出泄漏管道和非泄漏管道，再结合逆瞬态法对发生泄漏的管道的泄漏点进行定位。
S1：判断管道是否发生泄漏。
步骤A1，对实验的U形管道进行建模，尺寸设置与实验装置相同，两端管段长3.2m，中间管段1.6m，管径长0.0456m。泄漏孔位于距离管道首端2m处，泄漏孔径为1mm；网格划分采用结构化的四面体网格划分区域，并对泄漏口附近进行适当的网格加密处理，泄漏口处网格加密见图4。
步骤A2，控制方程和算法选择：
数值模型选用标准k-ε双方程模型，在模型中，用湍能k反映特征速度，用湍能耗散率ε反映特征长度尺度，利用Boussinesq假定简化而形成了k方程和ε方程，联合上述两个方程一道组成封闭的方程组对流动进行模拟，具体为：
k方程：
Figure BDA0002164450830000151
ε方程：
Figure BDA0002164450830000152
其中：xi，xj，xk为坐标分量；u为速度矢量；ρ为流体密度；μ为粘性系数；μt为湍动粘度；Gk、Gb为速度梯度、浮力引发的湍动能k的产生项；各常数C1ε＝1.44，C2ε＝1.92，C3ε＝0.09，σk＝1.0，σt＝1.3；Sk、Sε为用户定义源项；YM为可压湍流中脉动扩张项；
模型的数值计算算法选用SIMPLE算法，对于给定的压力，求解离散形式的动量方程，得出速度场，然后对给定的压力加以修正，反复迭代，直到获得收敛的解；
步骤A3，边界条件设置
入口边界条件设置：可压缩流体的入口边界条件应设为压力入口，本试验进行了三种入口压力下的数值模拟，分别设为0.3MPa、0.2MPa和0.1MPa；
出口边界条件设置：出口边界条件设为压力出口，并且管道末端出口和泄漏孔的末端出口都设置为外界环境压力；
流体介质设置：空气内介质为空气，空气的物理性质见表1；
管道材质设置：管道的材质为PE管，比热容为0.45cal/g℃，密度为0.956g/cm3，导热系数为0.5W/m·K。
表1常温下空气的物理性质表
Figure BDA0002164450830000161
步骤A4，在给定的多组进口压力下，分别对管道泄漏和非泄漏状态的管道内压力及流量参数的变化进行模拟，找出压力和流量参数的变化规律，对管道是否发生泄漏给出判断标准，然后将实际测得的管道数据与该标准进行对比，确定泄漏管道和无泄漏管道；
步骤A4.1，数值模拟管道压力、流速流场分布：
见附图5中的(a)图，在0.3MPa的进口压力下未泄漏时管道内的压力值从上游到下游总体呈线性分布，并保持着一定范围内的压力梯度随管道长度增加而减小；见附图5中的(b)图在同样的进口压力下，管道产生泄漏时，管道内总体压力下降幅度较大，最为明显的是在泄漏孔附近的压力相比未泄漏时下降了大约4.5×104Pa，并且在泄漏孔前后的压力梯度明显减小；见附图5中的(b)、(c)和(d)图，对比不同进口压力下管道泄漏状态的压力分布图，管道进口压力越大，管道内压力的变化梯度越大，在泄漏孔附近压力普遍减小并且泄漏孔前后管段的压力变化幅度减小。
见附图6中的a图，0.3MPa未泄漏时管道内流速趋于平稳状态，整体的变化幅度很小；见附图6中的b图，0.3MPa泄漏时管道外壁的流速明显增加，这是因为管道泄漏时，管道内外的压差加速了管道内气体的流动，从而导致管道外侧的流速与内侧的流速的差距变大；见附图6中的(b)、(c)和(d)图，对比不同进口压力下管道泄漏状态的速度矢量分布图，可知泄漏状态，管道进口压力越大，管道内整体流速越大，内外壁的流速差异越大。
步骤A4.2，实验管道压力、流速流场分布：
见附图7，正方形标记的折线表示上游流量，圆形标记的折线为下游流量，正三角标记的折线为上游压力，倒三角标记的折线为下游压力，如附图7的(a)、(b)和(c)图所示，在发生泄漏时，泄漏孔前端的流量大幅增高，泄漏孔后端的流量小幅度减小，过一段时间，泄漏孔前段流量减小，泄漏孔两端流量保持不变，且进口压力越大，流量差越大；产生泄漏时，整个管道压力减小，上下游压力差别不明显。
步骤A4.3，数值模拟和实验结果的比较分析
实验对模拟城市燃气管道的不同孔径泄漏前后管道状态变化的仿真和实验比较，可以得到总体相近的结论，但也有不同的结果。相近的结论：产生泄漏时，管道内整体压力降低，且压力梯度减小，管道内气体流速变快；不同的结论：实验中，泄漏孔前段管段流速加快，后段管段流速略微降低，模拟中整个管道内外壁流速相差较大，而泄漏孔前后的流速变化不明显。综合数值模拟和实验的分析，并结合相关文献，得到管道产生泄漏时参数变化的两个标准：(1)产生泄漏时，管道内整体压力降低，且压力梯度减小，管道内气体流速变快；(2)泄漏孔前段管段流速加快，后段管段流速略微降低，泄漏孔前后的流量差会在短时间内逐渐增大然后保持不变。
S2：管道泄漏实验：获取管道压力、流量数据
通过逐渐关闭阀门产生瞬态管流，阀门关闭时间持续足够长，以最小化不稳定和摩擦系数不确定性的影响，得到测量点X1、X2、X3、X4不同时段的压力、流量数据以及距离测量点两端Δx即1.6m处的压力、流量数据，见表1、2；
表10.3MPa泄漏各节点流量数据
Figure BDA0002164450830000171
Figure BDA0002164450830000181
表20.3MPa泄漏各节点压力数据
Figure BDA0002164450830000182
S3：建立管道泄漏控制方程。
假设节点为泄漏点，利用气体的连续性和运动方程求解出的非线性方程与泄漏孔泄漏量方程建立非金属管道气体泄漏状态的控制方程，如：
Figure BDA0002164450830000183
Figure BDA0002164450830000184
Figure BDA0002164450830000185
M1-M2＝ML (20)
P1＝P2＝PL (21)
其中PA、MA为节点前Δx处的压力，Pa、流量，kg·s-1；PB、MB为节点前Δx处的压力，Pa、流量，kg·s-1；P1、M1为流入泄漏点之前的一段距离的流体压力，Pa、流量，kg·s-1；P2、M2为流出泄漏点之后的一段距离的流体压力，Pa、流量，kg·s-1；PL为泄漏点处的压力，Pa；ML为泄漏孔泄漏流量，kg·s-1；Ae＝CA0为有效泄漏面积，m2；C为孔口系数，跟泄漏孔的形状有关；A0为泄漏孔面积，m2；Pd为管道起点压力，Pa；k为气体的绝热系数，无量纲；R为气体常数，J·kg-1·K-1；T为气体的温度，K；D为管内径，m；λ为摩阻系数；A为管道横截面积，m2。
S4：建立目标函数
S4.1：根据泄漏控制方程(17-21)表示出节点计算压力与真实压力的差值，如：
Figure BDA0002164450830000191
S4.2：确定泄漏面积的有效范围，从而确定泄漏量的上限，本文取泄漏孔半径为5mm时的有效泄漏面积作为泄漏面积的最大限制。
0≤Ae≤7.85×10-5m2(23)
S4.3：确定管道内的摩阻系数。
针对PE管道中气体不同的流动状态，对不同管材摩阻系数λ采用不同的计算公式，本文中采用PE管道计算公式是伯拉修斯公式。
伯拉修斯公式：
Figure BDA0002164450830000192
其中，雷诺系数计算公式采用经验公式：
Figure BDA0002164450830000201
Figure BDA0002164450830000202
式中，入为摩阻系数，Re为雷诺系数，Qv为体积流量m3/s，D为管道内径，m，v为运动粘度，m2/s，μ为动力粘度，Pa·s，ρ为气体密度，kg/m3。24摄氏度时空气的动力粘度为1.83×10-5Pa·s，运动粘度为1.4364×10-5m2/s。将未泄漏时稳态的数据及参数带入式(24-26)，求出管道内摩阻系数大约为0.01。
S5：算法优化
S5.1：将实验数据中第14至15组的节点X2的压力、流量数据及X2前后两端1.6m处的压力流量数据和已知参数代入式(22)，得到目标函数：
Figure BDA0002164450830000203
S5.2：定义非线性约束
M1-M2≥0(28)
0≤M1MB+Δt+MLmax(29)
MB+Δt≤M2MA+Δt(30)
MLmax是由泄漏量公式(3)代进最大限制的泄漏面积求出的最大泄漏量，为0.0258kg/s，得到：
Figure BDA0002164450830000204
Figure BDA0002164450830000205
S5.3：确定初始值X0＝[0.00343，0.00089]
S5.4：VLB＝[0，0.00089]，VUB＝[0.0293，0.00089]
S5.5：得出结果，M1＝0.0029；M2＝0.0009，即节点X2处前后两端计算出的流量差值为0.002kg/s，计算出泄漏面积6.09×10-7m2。按照步骤S5.1至5.5相应地计算出X1，X3，X4节点前后的流量差分别为0kg/s，0kg/s，0kg/s，并计算出泄漏面积分别为0m2，0m2，0m2。
S6：泄漏定位
S6.1，经过计算，节点X1，X3，X4前后的流量差为0，泄漏面积为0，故这三个节点都未产生泄漏；节点X2前后两端产生较明显的流量差，证明在X2节点前段1.6m至后段1.6m长的管段某个位置产生了泄漏，而X3节点计算结果正常，所以判断X3节点前段1.6m至后段1.6m长的管段没有产生泄漏，所以在X2节点前1.6m处至X2节点的管段再次进行节点的设置从而进行计算。
S6.2：继续划分节点X5，X6，见附图2，并根据步骤S5.1至5.4取Δx为0.4m时相应地计算出X5，X6节点前后的流量差，分别为0.002kg/s，0kg/s，计算出泄漏面积分别为6.09×10-7m2，0m2。证明在X5前后0.4m的管段存在泄漏，可依据需要再次设立节点继续计算，本实验X5节点为泄漏点，可以证明，该方法通过对监测节点压力、流量的设置计算，可以找取产生泄漏的管段，缺点是该方法过于依赖对节点的选取，并且模型的建立，参数的计算等方面都会对结果的准确性产生影响。
实施例二：
本实施例与实施例一的不同之处在于步骤S1中判断管道泄漏的方法不同，本实施例中采用基于马尔科夫链的方法实现非金属管道泄漏的判断。
S1：判断管道是否发生泄漏；基于马尔科夫链的非金属管道泄漏判断包括以下步骤：
步骤B1，利用安装在管线上的流量传感器采集到第一组五个流量：2.78、2.79、2.82、3.07、3.7，单位：m3/h，根据公式qi＝ΔQi/Qi×100％计算得到第一组五个流量变化率q1＝0％，q2＝3.58％，q3＝1.06％，q4＝8.14％，q5＝17.02％。
步骤B2，设置四个变化状态，第一状态：流量变化率qi的值为0，表示管道为非泄漏状态；第二状态：流量变化率qi的值为0-1(％)，表示管道为小型泄漏状态；第三状态：流量变化率qi的值为1-3(％)，表示管道为泄漏扩大状态；第四状态：流量变化率qi的值为3-100(％)，表示管道为大泄漏状态；
设定条件：管道的泄漏为突变过程，泄漏状态的变化为渐变过程，管道的工作状态由于其他各类因素的影响可能会直接从运行状态变为泄漏状态中某一程度的泄漏，且管道的泄漏状态不可逆，但是泄漏程度可逆，根据以上条件，得到管道泄漏的邻接矩阵A为：
Figure BDA0002164450830000221
其中，Aij＝1(i＝1，2...4；j＝1，2...4)表示第i状态与第j状态之间存在转移关系，Aij＝0表示第i状态与第j状态之间不存在转移关系。
将步骤1中获得第一组流量变化率的进行状态分类，将各流量变化率分别归为这四个状态分类，以时间为序，按时间顺序统计在这五个时间段中，流量变化率的状态转移情况，形成状态转移概率矩阵P中的各个元素。
状态转移概率矩阵P的计算过程如下：
步骤B2.1，根据分类得到处于第一状态的流量变化率数量为1，由第一状态至第一状态转移的数量为0，由第一状态至第二状态转移的数量为0，由第一状态至第三状态转移的数量为0，由第一状态至第四状态转移的数量为0，统计这5个流量变化率在下一时刻的状态转移情况：
第一状态至第一状态的转移概率p11＝0；
第一状态至第二状态的转移概率p12＝0；
第一状态至第三状态的转移概率p13＝1
第一状态至第四状态的转移概率p14＝0；
步骤B2.2，得到处于第二状态的流量变化率数量为0，由第二状态至第二状态转移的数量为0，由第二状态至第三状态转移的数量为0，由第二状态至第四状态转移的数量为0，统计这0个流量变化率在下一时刻的状态转移情况：
第二状态至第一状态的转移概率p21＝0；
第二状态至第二状态的转移概率p22＝0；
第二状态至第三状态的转移概率p23＝0；
第二状态至第四状态的转移概率p24＝0；
步骤B2.3，得到处于第三状态的流量变化率数量为1，由第三状态至第二状态转移的数量为0，由第三状态至第三状态转移的数量为0，由第三状态至第四状态转移的数量为1，统计这0个流量变化率在下一时刻的状态转移情况：
第三状态至第一状态的转移概率p31＝0；
第三状态至第二状态的转移概率p32＝0；
第三状态至第三状态的转移概率p33＝0；
第三状态至第四状态的转移概率p34＝1；
步骤B2.4，得到处于第四状态的流量变化率数量为3，由第四状态至第二状态转移的数量为0，由第四状态至第三状态转移的数量为0，由第四状态至第四状态转移的数量为2，统计这0个流量变化率在下一时刻的状态转移情况：
第四状态至第一状态的转移概率p41＝0；
第四状态至第二状态的转移概率p42＝0；
第四状态至第三状态的转移概率p43＝1/3
第四状态至第四状态的转移概率p44＝2/3
则该马尔科夫链中的状态转移概率矩阵P为：
Figure BDA0002164450830000241
步骤B3，收集第二组流量：4.79、6.62、7.16、7.66、8.22，单位：m3/h。得到第二组五个流量变化率：q1＝22.75q2＝27.64q3＝7.54％，q4＝6.52％，q5＝6.81％，同步骤2统计这五个流量变化率在每个区间内的数量，将此时实际流量变化率在每个区间的数量编为第一参数向量，设为α1＝[0 0 0 5]，其中α1表示第一参数向量，0表示在第一参数向量下，五个流量变化率中，在第一状态下的流量变化率的数量；0表示在第一参数向量下，五个流量变化率中，在第二状态下的流量变化率的数量；0表示在第一参数向量下，五个流量变化率中，在第三状态下的流量变化率的数量；5表示在第一参数向量下，五个流量变化率中，在第四状态下的流量变化率的数量；
由马尔科夫链计算，得到第三个时间段的基于马尔可夫链的预测流量变化率的参数向量β1，其公式为：
β1＝α1×P，得到的β1＝[0 0 0 5]，
由β1＝[0 0 0 5]，得到此时预测下一时刻处于第四状态的参数数量最多，含有5个参数，即
Figure BDA0002164450830000242
其中S＝4，则此时预测下一时刻状态为第四状态，即大泄漏状态。
步骤B4，收集到第三组流量：8.66、8.66、9.01、9.22、9.42，单位：m3/h，得到第三组五个流量变化率的值q1＝5.08％，q2＝0％，q3＝3.88％，q4＝2.33％，q5＝2.12％，同步骤2统计其在各个区间的匹配数量情况，设为α2＝[1 021]，其中得到此时处于第三状态的参数数量最多，含有2个参数，即
Figure BDA0002164450830000253
其中W＝2，即此时实际状态为第三状态。
步骤B5，比较
Figure BDA0002164450830000251
和
Figure BDA0002164450830000252
中S和W的数量值：
由于S＝4，W＝2，S＞W，得到处于开始泄漏后期阶段，且判断其状态为第三状态。
步骤B6，使第二组数据为第一组数据，第三组数据成为第二组数据，继续预测判定管道实时状态或泄漏情况。此时为泄漏第三状态，判定发生泄漏，触发报警程序，系统开始进行报警。
本例对泄漏定位的计算，即步骤S2至S6同实施例一，不再重复赘述说明。
总结：基于软件模拟的方法和马尔可夫法对管道的泄漏判断都有其优势和缺点，软件模拟的方法只能总结出泄漏和非泄漏两种状态流畅分布的规律，优点是对两种状态流场分布规律的分析，可以为泄漏的判断提供理论思路，给出理论依据，但在具体实践中，很难量化到具体数值，来判断管道是否发生微小泄漏；而马尔可夫法可以根据流量参数的捕捉对管道的泄漏情况进行判断，但是可能会由于其他因素对流量的影响从而产生误判，需要很多数据的支撑来保证对泄漏判断的准确性。
以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关的工作人员完全可以在不偏离本发明的范围内，进行多样的变更以及修改。本项发明的技术范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

Claims (6)

1.一种城市非金属管道泄漏定位方法，其特征在于：包括以下步骤：
S1：判断管道是否发生泄漏，通过数值模拟管道泄漏的特征规律分析或基于马尔科夫链的流量分析法对管道是否产生泄漏作出判断；
S2：管道泄漏实验：获取管道压力和流量数据；
通过逐渐关闭阀门产生瞬态管流，阀门关闭时间持续足够长，以最小化不稳定和摩擦系数不确定性的影响，得到各个监测点不同时段的压力和流量数据，即监测点两端△x处的压力PA、PB，流量MA、MB；
S3：建立管道泄漏控制方程；
利用气体的连续性和运动方程求解出的非线性方程与泄漏孔泄漏量方程建立非金属管道气体泄漏状态的控制方程；
S4：定义目标函数E；
将监测点处的压力计算值与测量值之差的最小值作为控制目标，控制方程和边界条件作为限制条件，对确定泄漏参数进行反问题分析，定义具有最小二乘目标准则的目标函数E，某个时间点距离监测点两端△x处测量点A、B的压力PA、PB和流量MA、MB由试验给出，单位时间步长△t后监测点两端的计算压力P1、P2和流量M1、M2是未知的，其中，P1＝P2＝PL，流量之差可以判断监测点是否产生泄漏，根据S3泄漏控制方程将已知的数据带入，即将步骤S2测得监测点两端△x处的压力PA、PB，流量MA、MB带入控制方程，得到关于监测点计算压力PL与流量M1、M2的方程，并且表示出监测点两端流量M1、M2关于计算压力PL的方程，以此为依据，根据逆瞬态计算原理，将某一时间点监测点计算压力PL与试验真实数据进行比较，并采用SQP算法对计算压力PL与真实压力的差值进行收敛达到最小值，得到符合真实情况的监测点流量数据，并将流量数据代入燃气泄漏量方程求出有效泄漏面积，从而确定监测点泄漏情况；并确定决策变量有效泄漏面积的限制范围，以及确定泄漏控制方程中摩阻系数；
S5：算法优化；
采取SQP方法在步骤S4的基础上设定监测点两端的流量限制范围，根据步骤S3所获得压力、流量数据将目标函数E最小化确定有效泄漏面积的泄漏参数；
S6：泄漏定位；
S6.1：在S5最小化结束时，如果发现节点的泄漏面积为非零，则将其视为泄漏；
S6.2：由于逆瞬态法的定位方法主要是判断监测点处是否泄漏，所以定位精度与监测点的设置数量有密切关系，并不能给出传统方法的定位公式，所以如果泄漏点并不在监测点处，即两个相邻的监测点的泄漏面积都不为零，则认为泄漏点在两个相邻且泄漏面积都不为零的监测点之间，并在两个监测点之间布置等距的多个测量节点，将测得的数据和计算的数据代入目标函数E重复S5，直至各节点距离小于管道检测长度的4％时终止重复步骤。
2.如权利要求1所述的城市非金属管道泄漏定位方法，其特征在于：步骤S1中通过数值模拟管道泄漏的特征规律分析管道发生泄漏判断方法具体包括：用AnsysIcem对管道进行建模，使用Fluent软件模拟气体在非金属管道泄漏和无泄漏状态下的压力、流量参数分布，并通过实验进行验证，分别找出气体在非金属管道泄漏和无泄漏状态下的压力、流量参数随时间和空间变化的规律，对管道是否发生泄漏给出判断标准；
步骤A1，对实验管道进行建模，采用结构化的四面体网格划分区域，并对泄漏口附近进行适当的网格加密处理；
步骤A2，控制方程和算法选择：
数值模型选用标准k-ε双方程模型，在模型中，用湍能k反映特征速度，用湍能耗散率ε反映特征长度尺度，利用Boussinesq假定简化而形成了k方程和ε方程，联合上述两个方程一道组成封闭的方程组对流动进行模拟，具体为：
k方程：
Figure FDA0003048061090000031
ε方程：
Figure FDA0003048061090000032
其中：xi，xj，xk为坐标分量；u为速度矢量；ρ为流体密度；μ为粘性系数；μt为湍动粘度；Gk、Gb为速度梯度、浮力引发的湍动能k的产生项；各常数C1ε＝1.44，C2ε＝1.92，C3ε＝0.09，σk＝1.0，σt＝1.3；Sk、Sε为用户定义源项；YM为可压湍流中脉动扩张项；
模型的数值计算算法选用SIMPLE算法，对于给定的压力，求解离散形式的动量方程，得出速度场，然后对给定的压力加以修正，反复迭代，直到获得收敛的解；
步骤A3，边界条件设置，包括对进口和出口的条件设置、流体和管道材质的物性参数设置；
步骤A4，在给定的多组进口压力下，分别对管道泄漏和非泄漏状态的管道内压力及流量参数的变化进行模拟，找出压力和流量参数的变化规律，对管道是否发生泄漏给出判断标准，然后将实际测得的管道数据与该标准进行对比，确定泄漏管道和无泄漏管道：
步骤A4.1，对数值模拟泄漏和非泄漏状态的管道内压力、流量流场进行分析，找出泄漏前后整个管道内压力、流量的变化规律；
步骤A4.2，对实验管道泄漏和非泄漏状态下的压力、流量流场进行分析，找出泄漏前后整个管道内压力、流量的变化规律；
步骤A4.3，将数值模拟和管道实验所得到的泄漏前后管道内压力、流量的变化规律进行对比分析，总结得到判断管道产生泄漏的具体标准。
3.如权利要求1所述的城市非金属管道泄漏定位方法，其特征在于：步骤S1中基于马尔科夫链的流量分析法判断管道发生泄漏方法具体包括：将非金属管道的运行状态分为正常运行状态和泄漏状态，管道由正常运行状态变为泄漏状态的过程为突变过程，且该过程不可逆，以此为依据，实时获取管道变化的转移概率矩阵，并预测管道流量变化的发展趋势，通过比较预测管道流量变化与实际管道流量变化，可判别管道的泄漏情况。
4.如权利要求1所述的城市非金属管道泄漏定位方法，其特征在于：步骤S3中所建立泄漏控制方程包括气体的连续性和运动方程求解出的非线性方程与泄漏孔泄漏量方程，具体如：
气体连续性和运动方程：
Figure FDA0003048061090000041
Figure FDA0003048061090000042
其中：M为气体流速，kg·S-1；A为管道横截面积，m2；P为压力，Pa；ρ为气体密度，kg·m-3；D为管内径，m；λ为摩阻系数；
采用特征线法对式(8)和(9)进行求解，得到两组特征线方程，用C+，C-进行表示，具体如：
Figure FDA0003048061090000043
Figure FDA0003048061090000044
Figure FDA0003048061090000045
Figure FDA0003048061090000051
其中，a为压力波速，m/s，气体管道等温流动时，a为定值；
公式(11)和(13)分别为C+，C-特征线AP和BP；公式(10)和(12)为满足各自特征线的相容性方程，将管道均分为N等份，步长为△x，时间步长为△t＝a△x，分别沿着C+特征线AP和C-特征线BP对公式(10)和(12)进行积分，其中对相容性方程的第三项即摩阻项采用二阶近似，得到C+特征线和C-特征线上的两个非线性方程，具体如：
Figure FDA0003048061090000052
Figure FDA0003048061090000053
管道泄漏视为小孔泄漏，则燃气泄漏公式取决于燃气在泄漏口处的流动速度，且燃气在泄漏口处的流动速度一般为亚音速流动，具体如：
Figure FDA0003048061090000054
其中：ML为泄漏孔泄漏流量，kg·s-1；Ae＝CA0为有效泄漏面积，m2；C为孔口系数，跟泄漏孔的形状有关；A0为泄漏孔面积，m2；Pd为管道起点压力，Pa；k为气体的绝热系数，无量纲；R为气体常数，J·kg-1·K-1；T为气体的温度，K；Mp为监测点处的流量，kg·s-1；Pp为监测点处的压力，Pa；
假设监测点为泄漏点，建立泄漏控制方程对监测点的压力、流量数据进行分析，以泄漏孔处为边界处，流入泄漏孔处之前△x距离内的流动特性满足公式(14)，流出泄漏孔处之后△x距离内的流动特性满足公式(15)；设流经泄漏孔处之前的流体参数下标为1，流经泄漏孔处之后的流体参数下标为2，则建立的泄漏控制方程具体如：
Figure FDA0003048061090000061
Figure FDA0003048061090000062
Figure FDA0003048061090000063
M1-M2＝ML (20)
P1＝P2＝PL (21)
其中P1、M1为流入泄漏点之前的一段距离的流体压力，Pa、流量，kg·s-1；P2、M2为流出泄漏点之后的一段距离的流体压力，Pa、流量，kg·s-1；PL为泄漏点处的压力，Pa。
5.如权利要求1所述的城市非金属管道泄漏定位方法，其特征在于：步骤S4中将监测点处压力的计算值与测量值之差的最小值作为控制目标，定义目标函数E并确定决策变量有效泄漏面积的限制范围，以及确定泄漏控制方程中摩阻系数，具体如：
S4.1：定义最小二乘准则目标函数：
Figure FDA0003048061090000064
其中：E是目标函数；M为时间步数；Pi是计算出的压力，Pa；P'i是测得的压力，Pa；通过将目标函数E最小化为零，从而产生最符合实际情况的有效泄漏面积，并根据有效泄漏面积的值确定该节点是否泄漏；
S4.2：确定有效泄漏面积的限制范围：
0≤Aei≤Aemax (23)
其中：Aei为节点i处的有效泄漏面积；Aemax为泄漏面积的最大限制，确定为管道横截面积的合理比例；
S4.3：摩阻系数根据实验数据由摩阻计算公式计算给出，而输气管道的流动一般处于紊流区，具体为：
伯拉修斯公式：
Figure FDA0003048061090000071
其中，雷诺系数计算公式采用经验公式：
Figure FDA0003048061090000072
Figure FDA0003048061090000073
式中，λ为摩阻系数，Re为雷诺系数，Qv为体积流量m3·s-1，D为管道内径，m，υ为运动粘度，m2·s-1，μ为动力粘度，Pa·s，ρ为气体密度，kg·m-3。
6.如权利要求1所述的城市非金属管道泄漏定位方法，其特征在于：步骤S5中使用SQP方法在S4的基础上设定的监测点两端的流量限制范围内对目标函数E进行收敛，算法的优化步骤具体如：
S5.1：建立M文件，根据实验数据定义目标函数E；
S5.2：建立M文件，定义约束条件中的非线性约束Ax≤b或Aeq·X＝Beq；
S5.3：确定迭代的初始值X0；
S5.4：确定变量的上下限VLB，VUB；
S5.5：建立主程序，非线性规划求解的函数为fmincon，运行求解。
CN201910742549.0A 2019-08-13 2019-08-13 一种城市非金属管道泄漏定位方法 Active CN110500511B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201910742549.0A CN110500511B (zh) 2019-08-13 2019-08-13 一种城市非金属管道泄漏定位方法
US17/254,285 US11592351B2 (en) 2019-08-13 2020-08-11 Urban non-metallic pipeline leakage location method
PCT/CN2020/108372 WO2021027803A1 (zh) 2019-08-13 2020-08-11 一种城市非金属管道泄漏定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910742549.0A CN110500511B (zh) 2019-08-13 2019-08-13 一种城市非金属管道泄漏定位方法

Publications (2)

Publication Number Publication Date
CN110500511A CN110500511A (zh) 2019-11-26
CN110500511B true CN110500511B (zh) 2021-06-22

Family

ID=68588085

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910742549.0A Active CN110500511B (zh) 2019-08-13 2019-08-13 一种城市非金属管道泄漏定位方法

Country Status (3)

Country Link
US (1) US11592351B2 (zh)
CN (1) CN110500511B (zh)
WO (1) WO2021027803A1 (zh)

Families Citing this family (25)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110500511B (zh) * 2019-08-13 2021-06-22 常州大学 一种城市非金属管道泄漏定位方法
CN112206455A (zh) * 2020-09-10 2021-01-12 云南省设计院集团有限公司 一种自喷管网埋地暗管漏损定位系统和方法
CN112115654B (zh) * 2020-09-27 2024-03-12 武汉第二船舶设计研究所(中国船舶重工集团公司第七一九研究所) 一种融合实测数据的管道流场压力分布评估方法及装置
CN113011636A (zh) * 2021-02-22 2021-06-22 北京市煤气热力工程设计院有限公司 城市热水供热管网水锤参数预测方法及装置
CN113052220A (zh) * 2021-03-16 2021-06-29 洛阳城市建设勘察设计院有限公司郑州工程分公司 直埋供热管道研究用密封性强度检测系统、终端、介质
CN113139353B (zh) * 2021-05-11 2023-09-29 东北大学 蒸汽管网动态计算及在线监测预警分析方法
CN113776952B (zh) * 2021-09-01 2024-06-11 武汉市政工程设计研究院有限责任公司 用于长距离大口径深层排水隧道的功能性试验方法及应用
CN113987696B (zh) * 2021-09-23 2024-04-30 西安交通大学 一种带破口高压气体容器临界流释放过程数值计算方法
CN114486108B (zh) * 2021-12-07 2023-03-31 南京航空航天大学 一种航空发动机短舱泄漏面积评估方法
KR102619454B1 (ko) * 2022-04-19 2023-12-29 한국건설기술연구원 이상 계측장치의 데이터 추정 시스템, 방법, 및 상기 방법을 실행시키기 위한 컴퓨터 판독 가능한 프로그램을 기록한 기록 매체
CN114777030B (zh) * 2022-05-23 2023-09-05 中用科技有限公司 一种基于nb-iot技术的危化气体监测方法
CN115164115A (zh) * 2022-07-05 2022-10-11 福州大学 基于物理模型驱动机器学习的燃气管道泄漏速率预测方法
CN115388343B (zh) * 2022-10-12 2024-04-16 广东海洋大学 一种高效的海洋油气管线泄漏检测与定位方法及系统
CN115618525B (zh) * 2022-11-11 2023-06-09 燕山大学 一种小型化非金属液压油箱设计方法
CN115978460B (zh) * 2022-11-28 2024-06-21 浙江英集动力科技有限公司 融合压力波和高精度授时的供热管道泄漏检测方法及系统
CN115824582B (zh) * 2023-02-09 2023-05-02 江苏新恒基特种装备股份有限公司 一种弯管流动性能测试方法、系统及存储介质
CN116300470B (zh) * 2023-04-06 2024-04-05 中国矿业大学 一种基于给水阀门通信故障的异步输出反馈控制方法
CN116562186B (zh) * 2023-05-10 2023-12-26 南京工程学院 基于模拟-优化的水下输气管道泄漏参数反演方法与系统
CN116557793B (zh) * 2023-07-10 2023-12-05 中建安装集团有限公司 一种融合压力传感和温度传感的供热管道运行状态监测系统及方法
CN116579272B (zh) * 2023-07-13 2023-09-08 天津市津安热电有限公司 一种基于稳态水力平差计算的管网渗漏定位方法和系统
CN116975767B (zh) * 2023-09-14 2023-12-08 长沙弘汇电子科技有限公司 基于大数据分析的智慧水务监测系统和灾害推演方法
CN117056864B (zh) * 2023-10-11 2024-01-23 广东力创信息技术有限公司 一种管道漏损预警方法及系统
CN117128777B (zh) * 2023-10-20 2024-01-19 湘潭新大粉末冶金技术有限公司 一种真空脱蜡烧结多气氛炉内安全预警系统
CN117451281B (zh) * 2023-11-15 2024-04-02 佛山市帆泰电器配件有限公司 一种空调连接管的密封性测试方法及系统
CN117781190A (zh) * 2024-01-03 2024-03-29 昆明理工大学 一种基于ARX-Laguerre模型管道泄漏检测与定位方法

Family Cites Families (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB0716027D0 (en) * 2007-08-16 2007-09-26 Ntnu Technology Transfer As Method and apparatus for opening a pipeline
CN101625071B (zh) * 2009-08-07 2012-11-28 天津大学 燃气管道泄漏检测和定位方法
CN101761780B (zh) * 2010-01-11 2012-12-26 中国石油大学(华东) 输气管道泄漏检测定位装置及其检测定位方法
CN102235575B (zh) * 2010-04-29 2013-12-25 国际商业机器公司 用于检查管道泄露的数据处理方法及系统
DE102010043482B4 (de) * 2010-11-05 2012-05-24 Siemens Aktiengesellschaft Leckageerkennung und Leckageortung in Versorgungsnetzen
US10161749B1 (en) * 2014-12-08 2018-12-25 Bentley Systems, Incorporated Optimizing water quality sensor placement for water distribution systems
CN105042339B (zh) * 2015-06-03 2017-06-23 中国石化销售有限公司华东分公司 一种基于无量纲的成品油管道泄漏量估计系统及方法
US10795382B2 (en) * 2016-08-02 2020-10-06 Sensus USA, Inc. Method and apparatus for model-based control of a water distribution system
CN106845437A (zh) * 2017-02-13 2017-06-13 常州大学 基于支持向量机回归的城市燃气管道泄漏定位方法
CN107013813B (zh) * 2017-05-27 2019-04-23 承德石油高等专科学校 一种供水管道泄漏量估算系统及方法
CN107590336B (zh) * 2017-09-13 2021-01-26 哈尔滨理工大学 燃气管道泄漏对内部流场影响的数值模拟方法
CN108019622B (zh) * 2018-02-05 2019-05-10 吉林大学 一种基于压力差的管道泄露定位的计算方法
CN108591836B (zh) * 2018-04-13 2020-06-26 中国石油大学(北京) 管道泄漏的检测方法和装置
CN108679458B (zh) * 2018-07-03 2021-03-19 安徽建筑大学 一种供水管网压力相关漏损定位方法
CN109033591B (zh) * 2018-07-14 2023-03-21 常州大学 基于逆瞬态模型的城市非金属管道泄漏定位方法
CN109210387B (zh) * 2018-09-27 2020-06-09 中国科学院合肥物质科学研究院 一种基于数学模型的气体管道泄漏检测定位的方法
CN109611696A (zh) * 2019-01-31 2019-04-12 西安建筑科技大学 一种管道泄漏检测与泄漏位置定位装置及方法
CN110500511B (zh) * 2019-08-13 2021-06-22 常州大学 一种城市非金属管道泄漏定位方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
分布式光纤油气长输管道泄漏检测及预警技术研究;曲志刚;《中国博士学位论文全文数据库 信息科技辑》;中国学术期刊（光盘版）电子杂志社;20090815;第I140-61页 *
基于SVR的城市燃气管道泄漏定位研究;王新颖等;《消防科学与技术》;20160215;第35卷(第2期);第283-286页 *

Also Published As

Publication number Publication date
US11592351B2 (en) 2023-02-28
US20220163421A1 (en) 2022-05-26
WO2021027803A1 (zh) 2021-02-18
CN110500511A (zh) 2019-11-26

Similar Documents

Publication Publication Date Title
CN110500511B (zh) 一种城市非金属管道泄漏定位方法
CN110197049B (zh) 一种基于瞬变反问题的非金属管道泄漏定位方法
CN105181040B (zh) 一种差压式流量计的数字化标定及优化方法
CN110513603B (zh) 一种基于逆瞬态分析法的非金属管道泄漏定位方法
Nakiboğlu et al. On the whistling of corrugated pipes: effect of pipe length and flow profile
Kostowski et al. Real gas flow simulation in damaged distribution pipelines
CN112628613A (zh) 一种管道泄漏监测、泄漏定位及泄漏量计算的方法及系统
CN109033591B (zh) 基于逆瞬态模型的城市非金属管道泄漏定位方法
CN105221933A (zh) 一种结合阻力辨识的管网泄漏检测方法
CN107014451A (zh) 基于广义回归神经网络推测超声波流量传感器系数的方法
Brahma Measurement and prediction of discharge coefficients in highly compressible pulsating flows to improve EGR flow estimation and modeling of engine flows
Tomor et al. Validation of a discrete model for flow distribution in dividing-flow manifolds: Numerical and experimental studies
Gan et al. k-factors for HVAC ducts: Numerical and experimental determination
CN109783972B (zh) 基于流固耦合分析计算的止回阀内泄漏流量的监测方法
Hao et al. The method for leakage detection of urban natural gas pipeline based on the improved ITA and ALO
Jiang et al. Study on pressure transients in low pressure water-hydraulic pipelines
US20210072060A1 (en) Real-time measurement of two-phase mass flow rate and enthalpy using pressure differential devices
Hao et al. An inverse transient nonmetallic pipeline leakage diagnosis method based on markov quantitative judgment
Mahood et al. Analytical and numerical investigation of transient gas blow down
Kiš et al. A CFD Analysis of Flow through a High-Pressure Natural Gas Pipeline with an Undeformed and Deformed Orifice Plate
CN106202743A (zh) 一种复杂汽、水管道压力取样位置的确定方法
Rezapour et al. Case study of leak detection based on Gaussian function in experimental viscoelastic water pipeline
Rup et al. Measurement of flow rate in square-sectioned duct bend
Shinde et al. Modelling and simulation of venturi parameters in relation to geometries and discharge coefficient with computational fluid dynamics techniques
Bartecki Transfer function models for distributed parameter systems: Application in pipeline diagnosis

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant