# An Inverse Transient Nonmetallic Pipeline Leakage Diagnosis Method Based on Markov Quantitative Judgment

Authors: Yongmei Hao; Yifei Ma; Juncheng Jiang; Zhixiang Xing; Lei Ni; Jian Yang
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

Aiming at the problems of early leakage monitoring of urban nonmetallic pipelines and the large positioning error, an inverse transient urban nonmetallic gas pipelines leakage location method based on Markov quantitative judgment was proposed. A Markov flow state transition probability matrix was established based on the flow data under different pressures obtained by experiments to quantitatively determine the pipeline leakage status. On this basis, an inverse transient leakage control equation suitable for urban nonmetallic gas pipeline leakage location was constructed according to the actual. The difference between the pressure and the calculated pressure was sought for the objective function. Finally, the objective function was optimized in conjunction with the sequential quadratic programming (SQP) method to obtain the actual leakage parameters and calculate the size and location of the leakage source. The results show that the inverse transient leakage localization method based on Markov’s quantitative judgment can more accurately determine the leakage status of the pipeline and calculate the early leakage parameters and leakage location of the gas pipeline, which improves the positioning accuracy.


---

## Body

## 1. Introduction

Nonmetallic pipelines have become the main materials for urban transportation pipelines in recent years due to their characteristics of corrosion resistance, high temperature resistance, low friction, and wear resistance. They have gradually replaced metal pipelines. According to relevant statistics, nonmetallic pipelines currently account for more than 70% of urban gas pipelines. Although nonmetallic pipelines have strong corrosion resistance, there are also problems with pipeline leakage caused by leakage and breakage at the pipe joints, fracture of the pipe body, shrinkage, and deformation of the inner liner of the pipeline. At present, the leakage monitoring of buried pipelines in most cities is still at the level of manual inspection. At the same time, due to the differences in pipeline materials, many methods for leak detection of metal pipelines are no longer applicable to nonmetallic pipelines. Therefore, it is necessary to study a leak diagnosis method and technology applicable to urban nonmetallic pipelines to strengthen the risk monitoring of urban nonmetallic pipelines.Leakage measurement of buried pipelines has always been a difficult subject to solve, especially in urban nonmetallic pipelines. There are still many problems to be studied in leak detection and positioning technology, which have caused widespread concern and achieved some results. In 2007, Li Heng et al. applied the Markov chain method to predict the aging of the anticorrosive layer of gas pipelines. The experiment proves that through long-term verification and regular summary, the accuracy of verification can be greatly improved, but the extracted features cannot accurately represent all useful information in the original data [1]; In 2008, Jung and Karney compared and pointed out that both the genetic algorithm and particle swarm optimization (PSO) algorithms tend to provide relatively accurate leakage calibration, but the calibration of the friction coefficient value is poor [2]. In 2008, Karney et al. combined the inverse transient method with the frequency domain analysis of system response, the analysis of reflected transient signals, and other methods that do not rely on transient models. The feasibility of the inverse transient method in the pipe network was verified by numerical simulation. It is mentioned that the positioning accuracy of the inverse transient method is related to the establishment of the model [3]. In 2011, Soares et al. combined hydraulic transient analysis with an optimized model through inverse transient analysis for nonmetallic pipe leakage monitoring and proposed that the key factor for positioning is the accurate calibration of the transient solver. The results show that, with the dispersion of the system model, the accuracy of the location of the leak is 4–15% of the total length of the pipeline [4]. In 2012, Haghighi et al. applied the central heuristic optimization (CFO) metaheuristic method to the detection of pipeline leakage specifications and pipeline friction coefficients, confirming that this method is better than genetic algorithm (GA), particle swarm algorithm (PSO), and ant colony optimization algorithm (ACO) [5]. In 2015, Huang et al. combined transient flow simulation and simulated annealing (SA) methods to develop the leak detection simulated annealing method (LDSA) and applied it to real-life pipelines and numerically simulated synthetic pipeline networks to verify its ability to identify leaks. The results show that this method can accurately detect the leakage of pores [6]. In 2016, Yu et al. combined gray theory with Markov theory to predict the corrosion state of pipelines. The results show that the method is simple in calculation and high in accuracy, but still has great limitations in engineering applications [7]. In 2016, Zhihong and Baomin aimed at the characteristics of nonlinear and random volatility of urban pipe network leakage rate data. Based on the gray Markov chain model, the pipeline leakage rate from 2005 to 2012 was the research object. A comparative analysis is performed with the gray Markov chain. The prediction results show that the prediction accuracy obtained by the modified Markov chain model is higher [8]. In 2017, Liu combined the correlation function method with ITA to determine the leakage point position by calculating the node sensitivity coefficient and then applied the GA algorithm to obtain the leakage amount through the inverse transient model. The feasibility of the combination of the two methods was verified through a large number of experiments [9]. In 2017, Cai Yongqiao et al. established a mathematical model of a gas-liquid two-phase flow leakage system based on a two-fluid model and confirmed that the numerical calculation results of the system are consistent with the experimental data trends, which proves that the model has a good applicability [10]. In 2017, Abdulshaheed proved that the biggest problem in the use of the inverse transient method is that the modeling of transients in the system does not accurately reflect the actual situation, and the calculated data and measured experimental data are very different. The difference is mainly because the model does not accurately consider the multidimensional effect of unstable friction [11]. In 2017, Chen et al. tested the transient model of pipeline leakage through the characteristic line method and confirmed that the positioning results calculated using the starting and ending pressure and flow results at the steady state after pipeline leakage are close to the real values [12]. In 2019, Liu et al. proposed a new method for leak detection based on Markov feature extraction and two-stage. The experimental results show that the method has higher accuracy and lower false alarm rate, but the calculation model needs to be adjusted for different pipelines, and the adaptability is not strong [13].The above studies have promoted the research of nonmetallic pipeline leakage monitoring, but most of them are still in the laboratory research stage, and there are also many problems such as incomplete acquisition of pipeline leakage information, the inability to quantitatively analyze and judge the pipeline leakage status, the inaccuracy of the inverse transient model leading to large calculation error, and poor applicability. In this paper, for gas pipelines in cities, a method for locating inverse transient nonmetallic pipelines based on quantitative judgment is explored. By establishing a Markov flow state transition probability matrix to quantitatively determine pipeline leakage, based on this, an inverse transient leakage control equation suitable for urban gas pipeline positioning analysis is constructed, and an objective function is established based on the difference between the actual pressure and the calculated pressure and finally combined with the sequential quadratic programming (SQP) method to accurately calculate the location of the pipeline leak source.


## 2. Markov Chain Quantitative Judgment Principle

The Markov Chain method can predict whether the future state of the pipeline will leak through the data of the current state of the pipeline. Since the leak state is irreversible, this method can be used for pipeline leak monitoring. The idea of applying it to the judgment of pipeline leakage is to convert the original flow data into a parameter vector of the Markov chain and extract the transfer characteristics of the Markov chain to build a state transition probability matrix to predict the trend of the pipeline flow. The rate of change of the flow rate is compared with the rate of change of the actual pipeline flow rate to determine whether the pipeline is leaking or not.


### 2.1. Establishment of a State Transition Probability Matrix

The first set of flow data is collected through the flow sensor on the pipeline and theK flow rate (q1,q2,q3,q4,…,qK) changes are calculated. The flow change rate is defined as follows:(1)qi=ΔQnQn×100%,qii=1,2,3,4,5,…,K.Among them,i is an ordering of the collected K flow change rates, △Qn is the change amount of the flow at the time point Tn and the flow at the time point Tn−1, and Qn is the real-time flow in the time period T. The qi includes four changing states. In the first state, the value of qi is 0, indicating a nonleak state; in the second state, the value of qi is 0-1 (%), indicating a small leakage state; in the third state, the value of qi is 1–3 (%), which indicates the leakage expansion state; in the fourth state, the value of qi is 3–100 (%), which indicates the large leakage state. The states of the first group of flow change rates are classified according to the above four states, and the probability of transition from one state to another among the K flow change rates is counted in chronological order as each element of the P matrix.The number of flow change rates for the four states isn1, n2, n3, and n4; the number of transitions from the first state to the first, second, third, and fourth states is, respectively, n11, n12, n13, and n14, and the transition probabilities are p11=n11/n1, p12=n12/n1, p13=n13/n1, and p14=n14/n1; the number of transitions from the second state to the second, third, and fourth states is, respectively, n22, n23, and n24, and the transition probabilities are p22=n22/n2, p23=n23/n2, and p24=n24/n2; the number of transitions from the third state to the second, third, and fourth states is, respectively, n32, n33, and n34, and the transition probabilities are p32=n32/n3, p33=n33/n3, and p34=n34/n3; the number of transitions from the fourth state to the second, third, and fourth states is, respectively, n42, n43, and n44, and the transition probabilities are p42=n42/4, p43=n43/n4, and p44=n44/4.Then, the state transition probability matrixP in the Markov chain is given by(2)P=p11p12p13p140p22p23p240p32p33p340p42p43p44,and satisfies: n1+n2+n3+n4=K.


### 2.2. Quantitative Judgment of Pipeline Leakage

(1)
As in Section2.1 and under the same environmental conditions, the second set of flow data is collected, and the number of actual flow rate changes in each interval at this time is set to the α1 parameter vector as follows:
(3)
α
1
=
N
1
1
N
2
1
N
3
1
N
4
1
.
Among them,α1 represents the first parameter vector, and N11, N21, N31, and N41 represent the quantity of the flow change rates in the K flow change rates in the first state, the second state, the third state, and the fourth state, respectively.
(2)
The Markov chain is calculated to obtain aβ1 parameter vector based on the predicted flow change rate based on the Markov chain in the future time period, and its formula is given by
(4)
β
1
=
α
1
×
P
=
m
1
1
m
2
1
m
3
1
m
4
1
,
wherem11, m21, m31, and m41, respectively, represent the quantity of flow change rates in the first state, the second state, the third state, and the fourth state, which are obtained by prediction.
Among theseK flow parameters, S is used as the predicted state serial number, (S = 1, 2, 3, 4), and the state with the largest number of flow parameters is selected, that is, NSmax is selected as the predicted pipeline state at this time, presented as the Smax-th state.
(3)
Similarly to [2], the third set of flow data is collected. Repeat the procedure described in Section 2.1, and the number of actual flow change rates in each interval at this time is set to the α2 parameter vector as follows:
(5)
α
2
=
N
1
2
N
2
2
N
3
2
N
4
2
.
Among theseK flow parameters, W is represented as the actual pipeline state serial number (W = 1, 2, 3, 4), and the state with the largest number of flow parameters is selected, that is, NWmax is selected as the pipeline is in the Wmax-th state.
(4)
The value of the state sequence numbersS in the predicted state NSmax is compared with the value of the state sequence numbers W in the actual state NWmax:
① IfS = W = 1, it is judged that the pipeline is in normal operation
② IfS < W, it is judged that the pipeline is in the early stage of starting leakage
③ IfS > W, it is judged that the pipeline is in the late stage of starting leakage
④ IfS = W ≠ 1, it is judged that the pipeline is in the leakage stabilization stage


## 2.1. Establishment of a State Transition Probability Matrix

The first set of flow data is collected through the flow sensor on the pipeline and theK flow rate (q1,q2,q3,q4,…,qK) changes are calculated. The flow change rate is defined as follows:(1)qi=ΔQnQn×100%,qii=1,2,3,4,5,…,K.Among them,i is an ordering of the collected K flow change rates, △Qn is the change amount of the flow at the time point Tn and the flow at the time point Tn−1, and Qn is the real-time flow in the time period T. The qi includes four changing states. In the first state, the value of qi is 0, indicating a nonleak state; in the second state, the value of qi is 0-1 (%), indicating a small leakage state; in the third state, the value of qi is 1–3 (%), which indicates the leakage expansion state; in the fourth state, the value of qi is 3–100 (%), which indicates the large leakage state. The states of the first group of flow change rates are classified according to the above four states, and the probability of transition from one state to another among the K flow change rates is counted in chronological order as each element of the P matrix.The number of flow change rates for the four states isn1, n2, n3, and n4; the number of transitions from the first state to the first, second, third, and fourth states is, respectively, n11, n12, n13, and n14, and the transition probabilities are p11=n11/n1, p12=n12/n1, p13=n13/n1, and p14=n14/n1; the number of transitions from the second state to the second, third, and fourth states is, respectively, n22, n23, and n24, and the transition probabilities are p22=n22/n2, p23=n23/n2, and p24=n24/n2; the number of transitions from the third state to the second, third, and fourth states is, respectively, n32, n33, and n34, and the transition probabilities are p32=n32/n3, p33=n33/n3, and p34=n34/n3; the number of transitions from the fourth state to the second, third, and fourth states is, respectively, n42, n43, and n44, and the transition probabilities are p42=n42/4, p43=n43/n4, and p44=n44/4.Then, the state transition probability matrixP in the Markov chain is given by(2)P=p11p12p13p140p22p23p240p32p33p340p42p43p44,and satisfies: n1+n2+n3+n4=K.


## 2.2. Quantitative Judgment of Pipeline Leakage

(1)
As in Section2.1 and under the same environmental conditions, the second set of flow data is collected, and the number of actual flow rate changes in each interval at this time is set to the α1 parameter vector as follows:
(3)
α
1
=
N
1
1
N
2
1
N
3
1
N
4
1
.
Among them,α1 represents the first parameter vector, and N11, N21, N31, and N41 represent the quantity of the flow change rates in the K flow change rates in the first state, the second state, the third state, and the fourth state, respectively.
(2)
The Markov chain is calculated to obtain aβ1 parameter vector based on the predicted flow change rate based on the Markov chain in the future time period, and its formula is given by
(4)
β
1
=
α
1
×
P
=
m
1
1
m
2
1
m
3
1
m
4
1
,
wherem11, m21, m31, and m41, respectively, represent the quantity of flow change rates in the first state, the second state, the third state, and the fourth state, which are obtained by prediction.
Among theseK flow parameters, S is used as the predicted state serial number, (S = 1, 2, 3, 4), and the state with the largest number of flow parameters is selected, that is, NSmax is selected as the predicted pipeline state at this time, presented as the Smax-th state.
(3)
Similarly to [2], the third set of flow data is collected. Repeat the procedure described in Section 2.1, and the number of actual flow change rates in each interval at this time is set to the α2 parameter vector as follows:
(5)
α
2
=
N
1
2
N
2
2
N
3
2
N
4
2
.
Among theseK flow parameters, W is represented as the actual pipeline state serial number (W = 1, 2, 3, 4), and the state with the largest number of flow parameters is selected, that is, NWmax is selected as the pipeline is in the Wmax-th state.
(4)
The value of the state sequence numbersS in the predicted state NSmax is compared with the value of the state sequence numbers W in the actual state NWmax:
① IfS = W = 1, it is judged that the pipeline is in normal operation
② IfS < W, it is judged that the pipeline is in the early stage of starting leakage
③ IfS > W, it is judged that the pipeline is in the late stage of starting leakage
④ IfS = W ≠ 1, it is judged that the pipeline is in the leakage stabilization stage


## 3. Inverse Transient Method for Nonmetallic Pipeline Leak Location Method

Based on the Markov chain method to quantitatively determine the leakage of the pipeline, the inverse transient method is used to further determine whether each monitoring node has a leakage, so as to delineate the pipeline where the leakage is at the node or between two nodes, and continuously shorten the node distance for precise positioning. The principle of inverse transient method positioning is as follows.


### 3.1. Governing Equation

Gas continuity and equation of motion [14]:(6)∂M∂t+1A∂∂xM2ρ+A∂p∂x+λ2DAρMM=0,(7)A∂ρ∂t+∂M∂x=0,where M is the gas flow rate, kg·s−1; A is the cross-sectional area of the pipe, m2; P is the pressure, Pa; ρ is the gas density, kg·m−3; D is the inner diameter of the pipe, m; λ is the friction coefficient.Solve the equations (6) and (7) by the characteristic line method, and obtain two sets of characteristic line equations, which are represented by C+ and C−:(8)C+:1adpdt+1AdMdt+λa2MM2DA2P=0,(9)dxdt=a,(10)C−:−1adpdt+1AdMdt+λa2MM2DA2P=0,(11)dxdt=−a,where a is the pressure wave velocity, m/s, and when the gas pipe is isotherm, a is a fixed value.The characteristic line grid is shown in Figure1. Equations (9) and (11) are C+ and C− characteristic lines AP and BP, respectively; equations (8) and (10) are the compatibility equations for satisfying the respective characteristic lines. The pipes are equally divided into N equal parts, the step size is Δx, and the time step is Δt = aΔx. Equations (8) and (10) are integrated along the C+ characteristic line AP and the C− characteristic line BP, respectively, wherein the third term of the compatibility equation is the second-order approximation, and the C+ characteristic line and the C− characteristic line are obtained. Two nonlinear equations on the feature line are as follows:(12)Pp−PA+aAMP−MA+λa2Δx4DA2MAMAPA+MPMPPP=0,PB−PP+aAMP−MB+λa2Δx4DA2MBMBPB+MPMPPP=0.Figure 1
Schematic diagram of the feature grid.The early leakage of urban pipelines can be regarded as small hole leakage, and the gas leakage formula depends on the flow velocity of the gas at the leakage port. Since the pressure of the city gas pipeline is relatively low, the flow velocity of the gas at the leak port is generally subsonic flow:(13)ML=AePdkRT2k+1k+1/k−1,where ML is the leaking leakage flow, kg·s−1; Ae is the effective leakage area, m2, and Ae=C×A0; C is the orifice coefficient, which is related to the shape of the leakage hole; A0 is the leaking hole area, m2; Pd is the starting point pressure of the pipeline, Pa; k is the adiabatic coefficient of the gas, dimensionless; R is the gas constant, J·kg−1·K−1; T is the temperature of the gas, K.Assuming that the monitoring node is a leak point, a leak control equation is established to analyze the pressure and flow data of the monitoring node. Taking the leak hole as a node, the flow characteristic within the Δx distance before flowing into the leak hole satisfies formula (8), and the flow characteristic within the Δx distance after flowing out of the leak hole satisfies formula (10), before flowing through the leak hole. The fluid parameter is subscripted as 1, and the fluid parameter after flowing through the leak hole is subscripted as 2, and the established leakage control equation is as follows:(14)P1−PA+aAM1−MA+λa2Δx4DA2MAMAPA+M1M1P1=0,(15)PB−P2+aAM2−MB+λa2Δx4DA2MBMBPB+M2M2P2=0,(16)ML=AePdkRT2k+1k+1/k−1,(17)M1−M2=ML,(18)P1=P2=PL,where P1 and M1 are the fluid pressures at a distance before flowing into the leak point, Pa, flow, kg·s−1; P2 and M2 are the fluid pressures at a distance after the leak point, Pa, flow, kg·s−1; PL is the pressure at the point of leakage, Pa.


### 3.2. Establishment of the Objective Function

(1)
The least squares criterion objective function is defined as follows:
(19)
E
=
min
∑
i
=
1
M
P
i
−
P
i
′
2
,
whereE is the objective function; M is the time step; Pi is the calculated pressure, Pa; Pi′ is the measured pressure, Pa. By minimizing the objective function E to zero, the effective leakage area that best fits the actual situation is generated, and whether the node leaks is determined according to the value of the effective area.
(2)
The friction coefficient in the pipe is determined as follows:
The friction coefficient is calculated based on the experimental data from the friction calculation formula:
Blasius formula:
(20)
λ
=
0.3164
Re
0.25
.
Among them, empirical formula is used in the Reynolds coefficient calculation formula.
Empirical formula:
(21)
Re
=
0.354
×
Q
v
D
υ
,
υ
=
μ
ρ
,
where λ is the friction coefficient, Re is the Reynolds coefficient,Qv is the volume flow rate, m3/s; D is the inner diameter of the pipe, m; ν is the kinematic viscosity, m2/s; μ is the dynamic viscosity, Pa·s; ρ is the gas density, kg·m−3.
(3)
The reasonable range of leakage area is determined as follows:
(22)
0
≤
A
e
i
≤
A
e
max
,
whereAei is the effective leakage area at node i; Aemax is the maximum limit of the leakage area, determined as a reasonable proportion of the cross-sectional area of the pipeline.


### 3.3. Algorithm Optimization

Sequential quadratic programming (SQP) algorithm is one of the most effective methods for dealing with nonlinear programming problems. It has the advantages of high computing efficiency and good convergence. Because the objective function contains the programming problem of the nonlinear function, it is optimized and converged by the SQP algorithm. The SQP algorithm steps are as follows:(1)
Define the objective functionE (x).(2)
Define the nonlinear constraints in the constraints.(3)
Determine the initial valueX0 of the iteration.(4)
Determine the upper and lower limits of the variables and express them as VUB and VLB, respectively.(5)
Establish the main program. The function of nonlinear programming is fmincon, and run the solution.(6)
At the end of the algorithm optimization, if the leak area of the node is found to be nonzero, it is regarded as a leaked state; otherwise it is regarded as a nonleakage state.(7)
Accurate positioning of pipeline leaks.


## 3.1. Governing Equation

Gas continuity and equation of motion [14]:(6)∂M∂t+1A∂∂xM2ρ+A∂p∂x+λ2DAρMM=0,(7)A∂ρ∂t+∂M∂x=0,where M is the gas flow rate, kg·s−1; A is the cross-sectional area of the pipe, m2; P is the pressure, Pa; ρ is the gas density, kg·m−3; D is the inner diameter of the pipe, m; λ is the friction coefficient.Solve the equations (6) and (7) by the characteristic line method, and obtain two sets of characteristic line equations, which are represented by C+ and C−:(8)C+:1adpdt+1AdMdt+λa2MM2DA2P=0,(9)dxdt=a,(10)C−:−1adpdt+1AdMdt+λa2MM2DA2P=0,(11)dxdt=−a,where a is the pressure wave velocity, m/s, and when the gas pipe is isotherm, a is a fixed value.The characteristic line grid is shown in Figure1. Equations (9) and (11) are C+ and C− characteristic lines AP and BP, respectively; equations (8) and (10) are the compatibility equations for satisfying the respective characteristic lines. The pipes are equally divided into N equal parts, the step size is Δx, and the time step is Δt = aΔx. Equations (8) and (10) are integrated along the C+ characteristic line AP and the C− characteristic line BP, respectively, wherein the third term of the compatibility equation is the second-order approximation, and the C+ characteristic line and the C− characteristic line are obtained. Two nonlinear equations on the feature line are as follows:(12)Pp−PA+aAMP−MA+λa2Δx4DA2MAMAPA+MPMPPP=0,PB−PP+aAMP−MB+λa2Δx4DA2MBMBPB+MPMPPP=0.Figure 1
Schematic diagram of the feature grid.The early leakage of urban pipelines can be regarded as small hole leakage, and the gas leakage formula depends on the flow velocity of the gas at the leakage port. Since the pressure of the city gas pipeline is relatively low, the flow velocity of the gas at the leak port is generally subsonic flow:(13)ML=AePdkRT2k+1k+1/k−1,where ML is the leaking leakage flow, kg·s−1; Ae is the effective leakage area, m2, and Ae=C×A0; C is the orifice coefficient, which is related to the shape of the leakage hole; A0 is the leaking hole area, m2; Pd is the starting point pressure of the pipeline, Pa; k is the adiabatic coefficient of the gas, dimensionless; R is the gas constant, J·kg−1·K−1; T is the temperature of the gas, K.Assuming that the monitoring node is a leak point, a leak control equation is established to analyze the pressure and flow data of the monitoring node. Taking the leak hole as a node, the flow characteristic within the Δx distance before flowing into the leak hole satisfies formula (8), and the flow characteristic within the Δx distance after flowing out of the leak hole satisfies formula (10), before flowing through the leak hole. The fluid parameter is subscripted as 1, and the fluid parameter after flowing through the leak hole is subscripted as 2, and the established leakage control equation is as follows:(14)P1−PA+aAM1−MA+λa2Δx4DA2MAMAPA+M1M1P1=0,(15)PB−P2+aAM2−MB+λa2Δx4DA2MBMBPB+M2M2P2=0,(16)ML=AePdkRT2k+1k+1/k−1,(17)M1−M2=ML,(18)P1=P2=PL,where P1 and M1 are the fluid pressures at a distance before flowing into the leak point, Pa, flow, kg·s−1; P2 and M2 are the fluid pressures at a distance after the leak point, Pa, flow, kg·s−1; PL is the pressure at the point of leakage, Pa.


## 3.2. Establishment of the Objective Function

(1)
The least squares criterion objective function is defined as follows:
(19)
E
=
min
∑
i
=
1
M
P
i
−
P
i
′
2
,
whereE is the objective function; M is the time step; Pi is the calculated pressure, Pa; Pi′ is the measured pressure, Pa. By minimizing the objective function E to zero, the effective leakage area that best fits the actual situation is generated, and whether the node leaks is determined according to the value of the effective area.
(2)
The friction coefficient in the pipe is determined as follows:
The friction coefficient is calculated based on the experimental data from the friction calculation formula:
Blasius formula:
(20)
λ
=
0.3164
Re
0.25
.
Among them, empirical formula is used in the Reynolds coefficient calculation formula.
Empirical formula:
(21)
Re
=
0.354
×
Q
v
D
υ
,
υ
=
μ
ρ
,
where λ is the friction coefficient, Re is the Reynolds coefficient,Qv is the volume flow rate, m3/s; D is the inner diameter of the pipe, m; ν is the kinematic viscosity, m2/s; μ is the dynamic viscosity, Pa·s; ρ is the gas density, kg·m−3.
(3)
The reasonable range of leakage area is determined as follows:
(22)
0
≤
A
e
i
≤
A
e
max
,
whereAei is the effective leakage area at node i; Aemax is the maximum limit of the leakage area, determined as a reasonable proportion of the cross-sectional area of the pipeline.


## 3.3. Algorithm Optimization

Sequential quadratic programming (SQP) algorithm is one of the most effective methods for dealing with nonlinear programming problems. It has the advantages of high computing efficiency and good convergence. Because the objective function contains the programming problem of the nonlinear function, it is optimized and converged by the SQP algorithm. The SQP algorithm steps are as follows:(1)
Define the objective functionE (x).(2)
Define the nonlinear constraints in the constraints.(3)
Determine the initial valueX0 of the iteration.(4)
Determine the upper and lower limits of the variables and express them as VUB and VLB, respectively.(5)
Establish the main program. The function of nonlinear programming is fmincon, and run the solution.(6)
At the end of the algorithm optimization, if the leak area of the node is found to be nonzero, it is regarded as a leaked state; otherwise it is regarded as a nonleakage state.(7)
Accurate positioning of pipeline leaks.


## 4. Nonmetallic Pipeline Leak Diagnosis Models and Procedures

### 4.1. Model Establishment

According to the characteristics of urban gas pipelines, a Markov chain flow state transition probability matrix is established for pipelines suspected of leaking status, and the predicted flow status is compared with the actual flow status to determine the pipeline leakage quantitatively. An inverse transient analysis model suitable for leakage location of nonmetallic gas pipelines in cities is established, and the objective function is established based on the difference between actual pressure and calculated pressure and is optimized with the sequential quadratic programming method (SQP) to ensure the calculation efficiency and convergence. The specific model is shown in Figure2 [15, 16].Figure 2
Leak diagnosis analysis model of urban buried nonmetallic pipelines.


### 4.2. Implementation Steps

Step 1: Markov quantitative leak determination
(1)
By gradually closing the pipeline valve to generate a transient flow, the first set of pipeline flow data is collected to establish a Markov state transition probability matrix
(2)
Under the same environment, the second set of flow data is collected, parameter vectors are calculated, and the flow state of the pipeline in the next period is predicted by combining the state transition probability matrix
(3)
The predicted pipeline flow state is compared with the actual state to judge whether the pipeline is leaking
Step 2: dividing nodes
If it is judged that the pipeline is leaking, the leaking pipe section will be divided into multiple monitoring nodes to obtain the flow and pressure parameters at the nodes
Step 3: establishment of the inverse transient positioning model
(1)
The inverse transient leakage control equation is established
(2)
An objective function based on near-steady-state flow data in pipeline leak monitoring is established
Step 4: algorithm optimization
The SQP algorithm is used to converge the objective function, to find the flow difference between the two ends of each node, and to calculate the effective leakage area of the monitoring node
Step 5: positioning calculation
If the effective leakage area of only a single monitoring node is nonzero, the monitoring node is judged to be a leakage point. If the effective leakage areas of multiple monitoring nodes are nonzero, the monitoring nodes of equal steps continue to be divided between nonzero nodes. Repeat steps 2 to 4 until the distance between the two nodes meets the requirements for locating the leakage of the pipeline.
The positioning accuracy is closely related to the number of monitoring nodes. If it is judged that a single node leaks, the leak point is considered to be located at that node. If it is calculated that there are multiple monitoring nodes leaking, continue to arrange multiple monitoring nodes with equal distances in the pipe section between these monitoring nodes, and repeat the calculation until the termination step is reached when the distance between each node meets the accuracy requirements. Find nodes with nonzero effective leakage area for accurate positioning.


## 4.1. Model Establishment

According to the characteristics of urban gas pipelines, a Markov chain flow state transition probability matrix is established for pipelines suspected of leaking status, and the predicted flow status is compared with the actual flow status to determine the pipeline leakage quantitatively. An inverse transient analysis model suitable for leakage location of nonmetallic gas pipelines in cities is established, and the objective function is established based on the difference between actual pressure and calculated pressure and is optimized with the sequential quadratic programming method (SQP) to ensure the calculation efficiency and convergence. The specific model is shown in Figure2 [15, 16].Figure 2
Leak diagnosis analysis model of urban buried nonmetallic pipelines.


## 4.2. Implementation Steps

Step 1: Markov quantitative leak determination
(1)
By gradually closing the pipeline valve to generate a transient flow, the first set of pipeline flow data is collected to establish a Markov state transition probability matrix
(2)
Under the same environment, the second set of flow data is collected, parameter vectors are calculated, and the flow state of the pipeline in the next period is predicted by combining the state transition probability matrix
(3)
The predicted pipeline flow state is compared with the actual state to judge whether the pipeline is leaking
Step 2: dividing nodes
If it is judged that the pipeline is leaking, the leaking pipe section will be divided into multiple monitoring nodes to obtain the flow and pressure parameters at the nodes
Step 3: establishment of the inverse transient positioning model
(1)
The inverse transient leakage control equation is established
(2)
An objective function based on near-steady-state flow data in pipeline leak monitoring is established
Step 4: algorithm optimization
The SQP algorithm is used to converge the objective function, to find the flow difference between the two ends of each node, and to calculate the effective leakage area of the monitoring node
Step 5: positioning calculation
If the effective leakage area of only a single monitoring node is nonzero, the monitoring node is judged to be a leakage point. If the effective leakage areas of multiple monitoring nodes are nonzero, the monitoring nodes of equal steps continue to be divided between nonzero nodes. Repeat steps 2 to 4 until the distance between the two nodes meets the requirements for locating the leakage of the pipeline.
The positioning accuracy is closely related to the number of monitoring nodes. If it is judged that a single node leaks, the leak point is considered to be located at that node. If it is calculated that there are multiple monitoring nodes leaking, continue to arrange multiple monitoring nodes with equal distances in the pipe section between these monitoring nodes, and repeat the calculation until the termination step is reached when the distance between each node meets the accuracy requirements. Find nodes with nonzero effective leakage area for accurate positioning.


## 5. Test Flow Field Analysis

The material of the test pipe section is a PPC pipe. The total length of the pipe is 8 m and it has a U-shaped structure. The diameter of the pipe is 0.0456 m, and the working pressures are 0.3 MPa, 0.2 MPa, and 0.1 MPa. The transmission medium is air at room temperature. The pipeline structure is shown in Figure3. An upstream flow sensor 1, an upstream pressure sensor 2, and an upstream infrasound sensor 3 are sequentially arranged at the inlet end of the pipeline. A downstream flow sensor 6, a downstream pressure sensor 5, and a downstream infrasound sensor 4 are sequentially arranged at the outlet end of the pipeline. The structure diagram is shown in Figure 4. X1 to X6 are monitoring nodes, the length of the X1-X2 tube section is 3.2 m, the length of the X2-X3 tube section is 1.6 m, and the length of the X3-X4 tube section is 3.2 m. The leakage hole is 0.6 m below the distance from X3, and the leakage aperture is 0.004 m. The X1 is at the pipe inlet and the X4 is at the pipe outlet. The X5 and X6 are 0.4 m and 1.2 m above the X2, respectively. The test ambient temperature is 24 degrees Celsius.Figure 3
Test pipeline.Figure 4
Schematic diagram of the test pipeline structure.Under the inlet pressure of 0.1 MPa, 0.2 MPa, and 0.3 MPa, the pressure and flow changes at the head and end of the pipeline under leaking and nonleaking conditions were collected. The collected data are shown in Figures5–7. The polygonal line marked with a square represents upstream flow, the polygonal line marked with a circle is downstream flow, the polyline marked with a positive triangle is upstream pressure, and the polyline marked with an inverted triangle is downstream pressure.Figure 5
Pressure and flow data collection at 0.1 MPa leakage.Figure 6
Pressure and flow data collection at 0.2 MPa leakage.Figure 7
Pressure and flow data collection at 0.3 MPa leakage.From Figures5–7, and repeat the above test several times: (1) When the pipeline is not leaking, the pressure at the beginning and end of the pipeline is similar, and the pressure of the entire pipeline is stable; when leaking, the overall pressure in the pipeline decreases, and the difference between upstream and downstream pressures is not obvious. (2)When there is no leakage, the flow velocity at the head and end of the pipeline is the same, and the flow velocity of the entire pipeline is stable; when leaking, the flow rate at the front end of the leak hole increases significantly, and the flow rate at the rear end of the leak hole decreases slightly. After a period of time, the flow at the front of the leak hole decreases, and the flow at both ends of the leak hole remains the same, and the greater the inlet pressure, the greater the flow difference.


## 6. Markov Chain Quantitative Leak Judgment

### 6.1. Establishment of a State Transition Probability Matrix

The first set of flow data was collected from the leakage test under an inlet pressure of 0.3 MPa in Section5 which includes 2.66, 2.79, 2.81, 3.16, and 3.76 (m3/h). It was obtained that the flow change rates at each time were q1 = 0, q2 = 4.66%, q3 = 0.71%, and q4 = 11.06%; the flow change rates in the four leakage states were n1 = 1, n2 = 1, n3 = 0, and n4 = 3.Then, the state transition probability matrixP in the Markov chain was as follows:(23)P=00010001000001/302/3.


### 6.2. Quantitative Determination of Pipeline Leakage

(1)
Under the same environmental conditions, the second set of flow data was collected which includes 4.36, 4.89, 5.33, 5.71, and 5.99 (m3/h). It was obtained that the flow change rates at each time were q1 = 13.67%, q2 = 10.84%, q3 = 8.26%, q4 = 6.65%, and q5 = 4.67%; the flow change rates in the four leakage states were n1 = 0, n2 = 0, n3 = 0, and n4 = 5. Coding the number of actual flow rate changes in each interval at this time as a parameter vector, we obtain
(24)
α
1
=
0
0
0
5
.
(2)
Calculated from the Markov chain, the parameter vectorβ1 of the predicted flow change rate based on the Markov chain in the future period was as follows:
(25)
β
1
=
0
5
/
3
0
10
/
3
.
N
S
max
=
N
4, it was predicted that the state of the pipeline will be in the fourth state.
(3)
As above, the third set of flow data was collected which includes 6.24, 6.65, 6.81, 7.02, and 7.12 (m3/h). It was obtained that the flow change rates at each moment were q1 = 4.01%, q2 = 6.17%, q3 = 2.35%, q4 = 2.99%, and q5 = 1.40%; the flow change rate in the four leakage states were n1 = 0, n2 = 0, n3 = 3, and n4 = 2. We obtain:
(26)
α
2
=
0
0
3
2
.
N
W
max
=
N
3, it was predicted that the state of the pipeline will be in the third state.
(4)
The value of the state sequence numbersS in the predicted state NSmax was compared with the value of the state sequence numbers W in the actual state NWmax:
Since S > W, it was judged that the state was the later stage of starting leakage. Similarly, the flow data under the inlet pressure of 0.1 MPa and 0.2 MPa were analyzed and calculated, and the results wereNSmax=N4 and NWmax=N3, respectively. That was, S > W, it was determined that the test pipeline was in the later stage of leakage.


## 6.1. Establishment of a State Transition Probability Matrix

The first set of flow data was collected from the leakage test under an inlet pressure of 0.3 MPa in Section5 which includes 2.66, 2.79, 2.81, 3.16, and 3.76 (m3/h). It was obtained that the flow change rates at each time were q1 = 0, q2 = 4.66%, q3 = 0.71%, and q4 = 11.06%; the flow change rates in the four leakage states were n1 = 1, n2 = 1, n3 = 0, and n4 = 3.Then, the state transition probability matrixP in the Markov chain was as follows:(23)P=00010001000001/302/3.


## 6.2. Quantitative Determination of Pipeline Leakage

(1)
Under the same environmental conditions, the second set of flow data was collected which includes 4.36, 4.89, 5.33, 5.71, and 5.99 (m3/h). It was obtained that the flow change rates at each time were q1 = 13.67%, q2 = 10.84%, q3 = 8.26%, q4 = 6.65%, and q5 = 4.67%; the flow change rates in the four leakage states were n1 = 0, n2 = 0, n3 = 0, and n4 = 5. Coding the number of actual flow rate changes in each interval at this time as a parameter vector, we obtain
(24)
α
1
=
0
0
0
5
.
(2)
Calculated from the Markov chain, the parameter vectorβ1 of the predicted flow change rate based on the Markov chain in the future period was as follows:
(25)
β
1
=
0
5
/
3
0
10
/
3
.
N
S
max
=
N
4, it was predicted that the state of the pipeline will be in the fourth state.
(3)
As above, the third set of flow data was collected which includes 6.24, 6.65, 6.81, 7.02, and 7.12 (m3/h). It was obtained that the flow change rates at each moment were q1 = 4.01%, q2 = 6.17%, q3 = 2.35%, q4 = 2.99%, and q5 = 1.40%; the flow change rate in the four leakage states were n1 = 0, n2 = 0, n3 = 3, and n4 = 2. We obtain:
(26)
α
2
=
0
0
3
2
.
N
W
max
=
N
3, it was predicted that the state of the pipeline will be in the third state.
(4)
The value of the state sequence numbersS in the predicted state NSmax was compared with the value of the state sequence numbers W in the actual state NWmax:
Since S > W, it was judged that the state was the later stage of starting leakage. Similarly, the flow data under the inlet pressure of 0.1 MPa and 0.2 MPa were analyzed and calculated, and the results wereNSmax=N4 and NWmax=N3, respectively. That was, S > W, it was determined that the test pipeline was in the later stage of leakage.


## 7. Calculation of Pipeline Leakage Location by Inverse Transient Method

### 7.1. Dividing Nodes

According to Figure4, four nodes X1, X2, X3, and X4 are first divided, and the step length was 1.6 m. Calculate the leakage parameters of 4 nodes, respectively.


### 7.2. Establishment of the Objective Function

(1)
According to the leakage control equation (15–19), the difference between the calculated pressure of the node and the real pressure was expressed as follows [17–19]:
(27)
Min
E
=
λ
a
2
Δ
x
M
1
2
+
M
2
2
/
4
D
A
2
P
A
−
P
B
+
a
/
A
M
A
+
M
B
−
a
/
A
M
1
+
M
2
−
λ
a
2
Δ
x
/
4
D
A
2
M
A
2
/
P
A
+
M
B
2
/
P
B
−
P
.
(2)
Determination of the range of the effective leakage area:
If the test leak hole radius was 2 mm, then take the effective leak area when the leak hole radius was 3 mm as the maximum limit of the leak area:
(28)
0
≤
A
e
≤
2.82
×
10
−
5
.
(3)
Determination of the friction coefficient in the pipeline:
At 24 degrees Celsius, the dynamic viscosity of air is 1.83 × 10−5 Pa·s and the kinematic viscosity is 1.4364 × 10−5 m2/s. Bring the steady state data and parameters without leakage into equation (20), and find the friction coefficient in the pipeline was about 0.07.


### 7.3. Algorithm Optimization

(1)
The data at the nodeX2 were calculated, and the pressure and flow data of the node X2, the pressure and flow data at △x (1.6 m) at both ends of X2, and the known parameters were substituted into equation (27). The test data at 0.3 MPa were taken as an example, and the transient data were brought into equation (27) to obtain the objective function:
(29)
Min
F
=
26641111404.11
×
M
1
2
+
M
2
2
733.26
−
208295.35
×
M
1
+
M
2
−
0.14
×
10
6
.
(2)
Definition of nonlinear constraints:
(30)
M
1
−
M
2
≥
0
,
M
2
≤
M
1
≤
M
B
+
Δ
t
+
M
L
max
,
0
t
≤
M
2
≤
M
A
+
Δ
t
.
M
Lmax is the maximum leakage amount obtained by substituting the maximum restricted leakage area into the leakage amount formula (16), which was 0.092 kg/s:
(31)
−
1
1
0
0
M
1
M
2
≤
0
0
,
0.000635
0
≤
M
1
M
2
≤
0.00893
0.000635
.
(3)
The flowMA and MB of the 14th group of data was taken as the initial calculated value, with the initial value X0 = [0.00198, 0.000635].
(4)
Determination of the upper and lower limits of variables:
The data were calculated according to Step 2: VLB = [0.000635, 0]; VUB = [0.00893, 0.000635].
(5)
The flow difference betweenX2 and X3 was calculated to be 0.0033 kg/s, and the data were brought into formula (16) to get the effective leakage area of 1.01 × 10−5m2. According to the above steps, the flow differences before and after nodes X1 and X4 were calculated, respectively, as 0 and 0, so the effective leakage areas of them were 0 and 0. The SQP algorithm iterative diagrams of 4 nodes are shown in Figures 8(a)–8(d).Figure 8
Iterative graph of algorithm for 4 nodes. (a)X1 node algorithm iteration, (b) X2 node algorithm iteration, (c) X3 node algorithm iteration, and (d) X4 node algorithm iteration.
(a)
(b)
(c)
(d)


### 7.4. Leak Location Calculation

(1)
According to the calculated effective leakage areas of the four nodes, the leakage point was located between nodesX2 and X3.
(2)
Since the leak point was between nodesX2 and X3, two monitoring points X5 and X6 were again set in the pipe section between nodes X2 and X3, and △x was taken as 0.4 m. According to the above process, the flow differences before and after X5 and X6 were calculated to be 0 and 0.0033 kg/s, and effective leakage areas were 0 and 1.01 × 10−5 m2. According to the requirements of 4.2 node division and the actual situation of the test pipeline, it was determined that the leak point was on a pipe section 0.8 m long before and after X6, that was, the final location of the leak point was at X6. Similarly, the data under the inlet pressure of 0.2 MPa and 0.1 MPa were analyzed according to the above steps, and finally the leak point was located at X6. The leakage parameters of this node are shown in Table 1.Table 1
The first test ofX5 node leakage parameters under different pressures.
Inlet pressure (MPa)
Flow difference (kg/s)
Node pressure (MPa)
Effective node leakage area (m2)
Real leakage area (m2)
Absolute error (m2)
0.3
0.0033
0.14
1.01 × 10−5
1.256 × 10−5
0.246 × 10−5
0.2
0.0029
0.11
1.124 × 10−5
1.256 × 10−5
0.132 × 10−5
0.1
0.0012
0.04
1.279 × 10−5
1.256 × 10−5
0.023 × 10−5Repeat the test two more times under the same conditions as the first test, and the exact calculation results of the leak point are shown in Tables2 and 3. The average absolute calculation errors of the effective leakage area at the leakage point under different inlet pressures obtained by 3 tests are shown in Table 4.Table 2
The second test ofX5 node leakage parameters under different pressures.
Inlet pressure (MPa)
Flow difference (kg/s)
Node pressure (MPa)
Effective node leakage area (m2)
Real leakage area (m2)
Absolute error (m2)
0.3
0.0027
0.1
1.151 × 10−5
1.256 × 10−5
0.105 × 10−5
0.2
0.0031
0.1
1.32 × 10−5
1.256 × 10−5
0.064 × 10−5
0.1
0.0017
0.06
1.208 × 10−5
1.256 × 10−5
0.048 × 10−5Table 3
The third test ofX5 node leakage parameters under different pressures.
Inlet pressure (MPa)
Flow difference (kg/s)
Node pressure (MPa)
Effective node leakage area (m2)
Real leakage area (m2)
Absolute error (m2)
0.3
0.0032
0.12
1.137 × 10−5
1.256 × 10−5
0.119 × 10−5
0.2
0.0024
0.09
1.137 × 10−5
1.256 × 10−5
0.119 × 10−5
0.1
0.0012
0.04
1.279 × 10−5
1.256 × 10−5
0.023 × 10−5Table 4
Absolute error under different pressures in three experiments.
Inlet pressure (MPa)
Absolute error (m2)
Average error (m2)
0.3
0.246 × 10−5
0.105 × 10−5
0.119 × 10−5
0.156 × 10−5
0.2
0.132 × 10−5
0.064 × 10−5
0.119 × 10−5
0.105 × 10−5
0.1
0.023 × 10−5
0.048 × 10−5
0.023 × 10−5
0.031 × 10−5It can be known from Tables1 to 4 that the absolute error for the effective area of the pipeline leakage is the lowest 0.023 × 10−5 m2 and the highest 0.156 × 10−5 m2. Under inlet pressures of 0.3 MPa, 0.2 MPa, and 0.1 MPa, the average absolute errors of effective leakage area calculated through multiple tests are 0.49 × 10−5 m2, 0.105 × 10−5 m2, and 0.031 × 10−5 m2. Compared with the result of the average absolute error of the antitransient method modified by Gaussian function in [20] being 0.52 × 10−5m2, it can be seen that this method can be effectively applied to urban nonmetallic natural gas pipelines. However, the flowing medium in the pipeline is a compressible fluid. The greater the pressure in the pipeline, the greater the change in the density of the fluid in the pipeline when a leak occurs, resulting in errors in the calculation of the model. Therefore, in practical applications, data that are close to steady state in the transient process should be selected for calculation.


## 7.1. Dividing Nodes

According to Figure4, four nodes X1, X2, X3, and X4 are first divided, and the step length was 1.6 m. Calculate the leakage parameters of 4 nodes, respectively.


## 7.2. Establishment of the Objective Function

(1)
According to the leakage control equation (15–19), the difference between the calculated pressure of the node and the real pressure was expressed as follows [17–19]:
(27)
Min
E
=
λ
a
2
Δ
x
M
1
2
+
M
2
2
/
4
D
A
2
P
A
−
P
B
+
a
/
A
M
A
+
M
B
−
a
/
A
M
1
+
M
2
−
λ
a
2
Δ
x
/
4
D
A
2
M
A
2
/
P
A
+
M
B
2
/
P
B
−
P
.
(2)
Determination of the range of the effective leakage area:
If the test leak hole radius was 2 mm, then take the effective leak area when the leak hole radius was 3 mm as the maximum limit of the leak area:
(28)
0
≤
A
e
≤
2.82
×
10
−
5
.
(3)
Determination of the friction coefficient in the pipeline:
At 24 degrees Celsius, the dynamic viscosity of air is 1.83 × 10−5 Pa·s and the kinematic viscosity is 1.4364 × 10−5 m2/s. Bring the steady state data and parameters without leakage into equation (20), and find the friction coefficient in the pipeline was about 0.07.


## 7.3. Algorithm Optimization

(1)
The data at the nodeX2 were calculated, and the pressure and flow data of the node X2, the pressure and flow data at △x (1.6 m) at both ends of X2, and the known parameters were substituted into equation (27). The test data at 0.3 MPa were taken as an example, and the transient data were brought into equation (27) to obtain the objective function:
(29)
Min
F
=
26641111404.11
×
M
1
2
+
M
2
2
733.26
−
208295.35
×
M
1
+
M
2
−
0.14
×
10
6
.
(2)
Definition of nonlinear constraints:
(30)
M
1
−
M
2
≥
0
,
M
2
≤
M
1
≤
M
B
+
Δ
t
+
M
L
max
,
0
t
≤
M
2
≤
M
A
+
Δ
t
.
M
Lmax is the maximum leakage amount obtained by substituting the maximum restricted leakage area into the leakage amount formula (16), which was 0.092 kg/s:
(31)
−
1
1
0
0
M
1
M
2
≤
0
0
,
0.000635
0
≤
M
1
M
2
≤
0.00893
0.000635
.
(3)
The flowMA and MB of the 14th group of data was taken as the initial calculated value, with the initial value X0 = [0.00198, 0.000635].
(4)
Determination of the upper and lower limits of variables:
The data were calculated according to Step 2: VLB = [0.000635, 0]; VUB = [0.00893, 0.000635].
(5)
The flow difference betweenX2 and X3 was calculated to be 0.0033 kg/s, and the data were brought into formula (16) to get the effective leakage area of 1.01 × 10−5m2. According to the above steps, the flow differences before and after nodes X1 and X4 were calculated, respectively, as 0 and 0, so the effective leakage areas of them were 0 and 0. The SQP algorithm iterative diagrams of 4 nodes are shown in Figures 8(a)–8(d).Figure 8
Iterative graph of algorithm for 4 nodes. (a)X1 node algorithm iteration, (b) X2 node algorithm iteration, (c) X3 node algorithm iteration, and (d) X4 node algorithm iteration.
(a)
(b)
(c)
(d)


## 7.4. Leak Location Calculation

(1)
According to the calculated effective leakage areas of the four nodes, the leakage point was located between nodesX2 and X3.
(2)
Since the leak point was between nodesX2 and X3, two monitoring points X5 and X6 were again set in the pipe section between nodes X2 and X3, and △x was taken as 0.4 m. According to the above process, the flow differences before and after X5 and X6 were calculated to be 0 and 0.0033 kg/s, and effective leakage areas were 0 and 1.01 × 10−5 m2. According to the requirements of 4.2 node division and the actual situation of the test pipeline, it was determined that the leak point was on a pipe section 0.8 m long before and after X6, that was, the final location of the leak point was at X6. Similarly, the data under the inlet pressure of 0.2 MPa and 0.1 MPa were analyzed according to the above steps, and finally the leak point was located at X6. The leakage parameters of this node are shown in Table 1.Table 1
The first test ofX5 node leakage parameters under different pressures.
Inlet pressure (MPa)
Flow difference (kg/s)
Node pressure (MPa)
Effective node leakage area (m2)
Real leakage area (m2)
Absolute error (m2)
0.3
0.0033
0.14
1.01 × 10−5
1.256 × 10−5
0.246 × 10−5
0.2
0.0029
0.11
1.124 × 10−5
1.256 × 10−5
0.132 × 10−5
0.1
0.0012
0.04
1.279 × 10−5
1.256 × 10−5
0.023 × 10−5Repeat the test two more times under the same conditions as the first test, and the exact calculation results of the leak point are shown in Tables2 and 3. The average absolute calculation errors of the effective leakage area at the leakage point under different inlet pressures obtained by 3 tests are shown in Table 4.Table 2
The second test ofX5 node leakage parameters under different pressures.
Inlet pressure (MPa)
Flow difference (kg/s)
Node pressure (MPa)
Effective node leakage area (m2)
Real leakage area (m2)
Absolute error (m2)
0.3
0.0027
0.1
1.151 × 10−5
1.256 × 10−5
0.105 × 10−5
0.2
0.0031
0.1
1.32 × 10−5
1.256 × 10−5
0.064 × 10−5
0.1
0.0017
0.06
1.208 × 10−5
1.256 × 10−5
0.048 × 10−5Table 3
The third test ofX5 node leakage parameters under different pressures.
Inlet pressure (MPa)
Flow difference (kg/s)
Node pressure (MPa)
Effective node leakage area (m2)
Real leakage area (m2)
Absolute error (m2)
0.3
0.0032
0.12
1.137 × 10−5
1.256 × 10−5
0.119 × 10−5
0.2
0.0024
0.09
1.137 × 10−5
1.256 × 10−5
0.119 × 10−5
0.1
0.0012
0.04
1.279 × 10−5
1.256 × 10−5
0.023 × 10−5Table 4
Absolute error under different pressures in three experiments.
Inlet pressure (MPa)
Absolute error (m2)
Average error (m2)
0.3
0.246 × 10−5
0.105 × 10−5
0.119 × 10−5
0.156 × 10−5
0.2
0.132 × 10−5
0.064 × 10−5
0.119 × 10−5
0.105 × 10−5
0.1
0.023 × 10−5
0.048 × 10−5
0.023 × 10−5
0.031 × 10−5It can be known from Tables1 to 4 that the absolute error for the effective area of the pipeline leakage is the lowest 0.023 × 10−5 m2 and the highest 0.156 × 10−5 m2. Under inlet pressures of 0.3 MPa, 0.2 MPa, and 0.1 MPa, the average absolute errors of effective leakage area calculated through multiple tests are 0.49 × 10−5 m2, 0.105 × 10−5 m2, and 0.031 × 10−5 m2. Compared with the result of the average absolute error of the antitransient method modified by Gaussian function in [20] being 0.52 × 10−5m2, it can be seen that this method can be effectively applied to urban nonmetallic natural gas pipelines. However, the flowing medium in the pipeline is a compressible fluid. The greater the pressure in the pipeline, the greater the change in the density of the fluid in the pipeline when a leak occurs, resulting in errors in the calculation of the model. Therefore, in practical applications, data that are close to steady state in the transient process should be selected for calculation.


## 8. Conclusion

(1)
By extracting the transfer characteristics of Markov chains, the transition probability matrix of pipeline changes can be obtained in real time, and the trend development of pipeline flow changes can be predicted. By comparing the predicted changes in the pipeline flow rate with the actual changes in the pipeline flow rate, the leakage of the pipeline can be judged. The Markov chain can monitor the changes in pipeline flow in real time and reduce the false alarm rate through prediction and comparison.
(2)
According to the characteristics of the gas pipeline gas medium, according to the gas motion and continuity equations, specific boundary conditions are given to establish an inverse transient location analysis mathematical model suitable for gas pipeline leakage monitoring. Selecting the characteristic line method combined with the sequential quadratic programming method to solve the pipeline leakage model has certain convergence and improves the accuracy of the model calculation.
(3)
The judgment error for the effective area of pipeline leakage can reach 0.023 × 10−5 − 0.156 × 10−5m2. This is because the flowing medium in the pipeline is a compressible fluid. The greater the pressure in the pipeline, the greater the change in the density of the fluid in the pipeline when a leak occurs, causing errors in the calculation of the model. Therefore, in practical applications, data close to steady state during the transient process should be selected for calculation.
(4)
Inverse transient leakage detection method is proposed based on Markov's quantitative judgment that when conducting pipeline leak detection, it is necessary to set as many equidistant nodes as possible to ensure high positioning accuracy.



---

*Source: 9527836-2020-04-22.xml*