Urban non-metallic pipeline leakage location method

Abstract

Disclosed is an urban non-metallic pipeline leakage location method, the method comprise the following steps: determining whether leakage occurs in a pipeline through numerical simulation law analysis or a Markov chain-based flow analysis method; for a leaking pipeline, establishing an inverse-transient control equation for non-metallic pipeline gas leakage, and obtaining pressure and flow rate data of each measuring point in different periods of time through experiments and substituting the data into the control equation to analyze experimental data; and defining a nonlinear programming problem of an objective function with a least squares criterion, and applying a sequential quadratic programming method to minimize the objective function, so as to determine the size and position of the leakage.

Claims

1. An urban non-metallic pipeline leakage location method, comprising the following steps: step S1: determining whether leakage occurs in a pipeline, which is determining whether the leakage occurs in the pipeline through characteristic law analysis based on numerical simulation of pipeline leakage or a Markov chain-based flow analysis method; step S2: conducting pipeline leakage experiment and obtaining pipeline pressure and flow rate data, which is producing transient pipe flow by gradually closing a valve, and maintaining a valve closing time long enough to minimize influence of instability and uncertainty of a friction coefficient, so as to obtain pressure and flow rate data of each of a plurality of measuring point in different periods of time; step S3: establishing a pipeline leakage control equation, which is establishing the pipeline leakage control equation of a non-metallic pipeline gas leakage state by using a nonlinear equation solved based on continuity and motion equations of a gas and a leakage amount equation of a leakage hole, substituting pressures P.sub.A and P.sub.B and flow rates M.sub.A and M.sub.B at both Δx-distance ends of a monitoring point measured in step S2 into the control equation to obtain equations for calculating pressure P.sub.L and flow rates M.sub.1 and M.sub.2 of the monitoring point, on this basis, comparing the calculated pressure P.sub.L of the monitoring point at a certain time point with real data of test according to a principle of inverse-transient calculation, and converging a difference between the two to a minimum value by using a Sequence Quadratic Programing (SQP) algorithm so as to obtain leakage parameters that conform to an actual situation; step S4: defining an objective function, which is performing inverse problem analysis on the obtained leakage parameters by taking a minimum value of a difference between a calculated value and a measured value of a pressure at the monitoring point as a control target and taking the control equation and boundary conditions as limiting conditions, to define the objective function with a least squares objective criterion, wherein pressures P.sub.A and P.sub.B and flow rates M.sub.A and M.sub.B of measuring points A and B at the both ends of the monitoring point α at a certain time point are given by the experiment, the calculated pressures P.sub.1 and P.sub.2 and flow rates M.sub.1 and M.sub.2 at the both ends of the monitoring point α are unknown after a unit time step Δt, and based on a difference in flow rate, whether a leakage occurs in the monitoring point α is determined; substituting data according to the leakage control equation in S3 and expressing an equation of the flow rates M.sub.1 and M.sub.2 on the pressure P.sub.L at the both ends of the monitoring point α, obtaining flow rate data of the monitoring point conforming to a real situation by converging a difference between the calculated pressure P.sub.L and a real pressure to the minimum, and substituting the flow rate data of the monitoring point into a gas leakage amount equation to obtain an effective leakage area, so as to determine a leakage situation of the monitoring point; step S5: performing algorithm optimization, which is setting flow constraint ranges at the both ends of the monitoring point in step S4 by using an SQP method, and minimizing the objective function according to the pressure and flow rate data obtained in step S3 to determine the leakage parameters of the effective leakage area; and step S6: leakage location positioning, which comprises the following step: S6.1: if a leakage area of a node is found to be non-zero at an end of a minimization in S5, regarding the node as a leaking node; and S6.2: if the leakage point is not at the monitoring point, that is, neither of the leakage areas of two adjacent monitoring points is zero, then considering that the leakage point is between two adjacent monitoring points of which neither of the leakage areas is zero, and equidistantly arranging a plurality of measuring nodes between the two adjacent monitoring points; substituting measured data and calculated data of a node of the plurality of measuring nodes as the monitoring point into the objective function; and repeating S5 and the equidistantly arranging and the substituting steps of S6.2 for a next plurality of measuring nodes between two previously arranged adjacent nodes until a distance between two adjacent nodes is less than 4% of a pipeline detection length, or until a leakage area of a node is found to be non-zero at an end of a minimization in S5.

2. The urban non-metallic pipeline leakage location method according to claim 1, wherein, method for determining an occurrence of leakage of the pipeline through characteristic law analysis based on numerical simulation of pipeline leakage in step S1 comprises: modeling the pipeline by using a modeling software, simulating pressure and flow rate parameter distributions, of the gas in non-metallic pipeline leakage and non-leakage states by using a fluid simulation software, and verifying same through experiments to respectively find, out laws of pressure and flow rate parameters of the gas changing with time and space in the non-metallic pipeline leakage and non-leakage states, so as to give standards for determining whether leakage occurs in the pipeline, wherein the method for determining an occurrence of leakage of the pipeline through characteristic law analysis based on numerical simulation of pipeline leakage in step S1 specifically comprises the following steps: step A1: modeling an experimental pipeline, using a structured tetrahedral mesh for region division, and performing proper mesh refinement near a leakage port; step A2: control equation and algorithm selection: choosing a standard k−ε double equation model as a numerical model; in the model, using turbulent energy k to reflect a characteristic velocity, using a turbulent energy dissipation rate ε to reflect a characteristic length scale, using Boussinesq hypothesis for simplification to form k equation and ε equation, and combining the two equations together to form a closed equation set to simulate a flow, specifically: $\begin{array}{cc}k\mathrm{equation}\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho {\mathrm{ku}}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\mathrm{\rho \epsilon }-Y_M+S_K& \left(1\right)\end{array}$ $\begin{array}{cc}\epsilon \mathrm{equation}\frac{\partial \left(\mathrm{\rho \epsilon }\right)}{\partial t}+\frac{\partial \left(\rho {\mathrm{ku}}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_t}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\frac{C_{1\epsilon }\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}+S_{\epsilon }& \left(2\right)\end{array}$ where x.sub.i and x.sub.j are coordinate components; u.sub.i a velocity vector; ρ is a fluid density; μ is a viscosity coefficient; μ.sub.i is a turbulent viscosity; t represents time; G.sub.k and G.sub.b are production terms of velocity gradient and turbulent kinetic energy k caused by buoyancy; constants C.sub.1ε=1.44, C.sub.2ε=1.92, C.sub.3ε=0.09, σ.sub.k=1.0, and σ.sub.t=1.3, S.sub.k and S.sub.εare user-defined source items; and Y.sub.M is a pulsating expansion term in compressible turbulent flow; choosing SIMPLE algorithm as a numerical calculation algorithm of the model; for a given pressure, solving a momentum, equation in a discrete form to obtain a velocity field, and then correcting the given pressure and iterating repeatedly until a convergent solution is obtained; step A3: boundary condition settings, comprising condition settings for an inlet and an outlet, and physical parameter settings fora fluid and a pipe material; and step A4: simulating changes in internal pressure and flow rate parameters of the pipeline in the pipeline leakage and non-leakage states respectively under given multiple sets of inlet pressures to find out change laws of pressure and flow rate parameters, so as to give the standards for determining whether leakage occurs in the pipeline, and then comparing an actual measured pipeline data with the standards to determine leaking pipelines and non-leaking pipelines: step A4.1: analyzing numerically simulated flow fields of internal pressure and flow rate of the pipeline in the leakage and non-leakage states to find out the change laws of the internal pressure and flow rate of the entire pipeline before and after the leakage; step A4.2: analyzing the flow fields of pressure and flow rate of the experimental pipeline in the leakage and non-leakage states to find out change laws of the internal pressure and flow rate of the entire pipeline before and after the leakage; and step A4.3: comparatively analyzing the change laws of the internal pressure and flow rate of the pipeline before and after the leakage obtained through the numerical simulation and the experiments, and summarizing the change laws to obtain specific standards for determining an occurrence of leakage of the pipeline.

3. The urban non-metallic pipeline, leakage location method according to claim 1, wherein, determining the occurrence of leakage of the pipeline through a Markov chain-based flow analysis method in step S1 specifically comprises: dividing an operation state of a non-metallic pipeline into a normal operation state and a leakage state, wherein a process of the pipeline changing from the normal operation state to the leakage state is a sudden change process, and the process is irreversible; on this basis, obtaining a matrix of transition probability of pipeline change in real time, and predicting a development trend of pipeline flow rate change; and determining the occurrence of leakage of the pipeline by comparing predicted pipeline flow rate change with an actual pipeline flow rate change.

4. The urban non-metallic pipeline leakage location method according to claim 1, wherein, the pipeline leakage control equation established in step S3 comprises the nonlinear equation solved based on continuity and motion equations of the gas and the leakage amount equation of the leakage hole, specifically as follows: the continuity and motion equations of the gas: $\begin{array}{cc}\frac{\partial M}{\partial t}+\frac{1}{A}\frac{\partial }{\partial x}\left(\frac{M^2}{\rho }\right)+A\frac{\partial p}{\partial x}+\frac{\lambda }{2DA\rho }M.Math.M.Math.=0& \left(8\right)\end{array}$ $\begin{array}{cc}A\frac{\partial \rho }{\partial t}+\frac{\partial M}{\partial x}=0& \left(9\right)\end{array}$ where M is a flow velocity of the gas, whose dimensional unit is kg.Math.S.sup.−1; A is a cross sectional area of the pipeline, whose dimensional unit is m.sup.2; P is a pressure, whose dimensional unit is Pa; ρ is a gas density, whose dimensional unit is kg.Math.m.sup.−3; D is an inner diameter of the pipeline, whose dimensional unit is m; and λ is a friction resistance coefficient; equations (8) and (9) are solved by using a characteristic line method to obtain two sets of characteristic line equations, which are expressed by C+ and C−, specifically as follows: $\begin{array}{cc}C+:\frac{1}{a}\frac{\mathrm{dp}}{dt}+\frac{1}{A}\frac{\mathrm{dM}}{\mathrm{dt}}+\frac{\lambda a^{2}M.Math.M.Math.}{2D{A}^{2}P}=0& \left(10\right)\end{array}$ $\begin{array}{cc}\frac{\mathrm{dx}}{dt}=a& \left(11\right)\end{array}$ $\begin{array}{cc}C-:-\frac{1}{a}\frac{\mathrm{dp}}{dt}+\frac{1}{A}\frac{\mathrm{dM}}{dt}+\frac{\lambda a^{2}M.Math.M.Math.}{2D{A}^{2}P}=0& \left(12\right)\end{array}$ $\begin{array}{cc}\frac{\mathrm{dx}}{dt}=-a& \left(13\right)\end{array}$ where α is a velocity of pressure wave, whose dimensional unit is m/s; and when the gas in the pipeline is flowing isothermally, α is a constant value; formulas (11) and (13) are C+ and C− characteristic lines AP and BP, respectively; formulas (10) and (12) are compatibility equations satisfying respective characteristic lines; the pipeline is equally divided into N equal parts, with a step length being Δx and a time step length being Δt=αΔx; t represents time and x represents length; formulas (10) and (12) are integrated along C+ characteristic line AP and C− characteristic line BP, respectively; second-order approximation is used for a third term of a compatibility equation, namely a friction term, to obtain two nonlinear equations on the C+ characteristic line and the C− characteristic line, specifically as follows: $\begin{array}{cc}{P}_p-P_A+\frac{a}{A}\left(M_P-M_A\right)+\frac{\lambda a^2\Delta x}{4D{A}^2}\left(\frac{M_A.Math.M_A.Math.}{P_A}+\frac{M_P.Math.M_P.Math.}{P_P}\right)=0& \left(14\right)\end{array}$ $\begin{array}{cc}{P}_B-P_P+\frac{a}{A}\left(M_P-M_B\right)+\frac{\lambda a^2\Delta x}{4D{A}^2}\left(\frac{M_B.Math.M_B.Math.}{P_B}+\frac{M_P.Math.M_P.Math.}{P_P}\right)=0& \left(15\right)\end{array}$ where M.sub.P is a flow rate at the monitoring point requiring a calculation of the leakage parameters, and P.sub.P is a pressure at the monitoring point requiring a calculation of the leakage parameters; if the pipeline leakage is regarded as small hole leakage, a gas leakage formula depends on the flow velocity of the gas at the leakage port, and the flow velocity of the gas at the leakage port is generally subsonic, specifically as follows: $\begin{array}{cc}{M}_L=A_e{P}_d\sqrt{\frac{k}{RT}{\left(\frac{2}{k+1}\right)}^{\frac{k+1}{k-1}}}& \left(16\right)\end{array}$ where M.sub.L is a leakage flow rate of the leakage hole, whose dimensional unit is kg.Math.s.sup.−1; A.sub.e=CA.sub.0 is an effective leakage area, whose dimensional unit is m.sup.2; C is an orifice coefficient, which is related to the shape of the leakage hole; A.sub.0 is an area of the leakage hole, whose dimensional unit is m.sup.2; P.sub.d is an initial pressure of the pipeline, whose dimensional unit is Pa; k is an adiabatic coefficient of the gas, dimensionless; R is a gas constant, whose dimensional unit is J.Math.kg.sup.−1; and T is a temperature of the gas, whose dimensional unit is K; assuming that the monitoring node is a leakage point, the pipeline leakage control equation is established to analyze the pressure and flow rate data of the monitoring node; taking the leakage hole as a boundary, flow characteristics within the Δx distance before flowing into the leakage hole satisfy formula (14) , and the flow characteristics within the Δx distance after flowing out of the leakage hole satisfy formula (15); supposing that a subscript of the fluid parameter before flowing through the leakage hole is 1, and a subscript of the fluid parameter after flowing through the leakage hole is 2, then the established pipeline leakage control equation is specifically as follows: $\begin{array}{cc}{P}_1-P_A+\frac{a}{A}\left(M_1-M_A\right)+\frac{\lambda a^2\Delta x}{4D{A}^2}\left(\frac{M_A.Math.M_A.Math.}{P_A}+\frac{M_1.Math.M_1.Math.}{P_1}\right)=0& \left(17\right)\end{array}$ $\begin{array}{cc}{P}_B-P_2+\frac{a}{A}\left(M_2-M_B\right)+\frac{\lambda a^2\Delta x}{4D{A}^2}\left(\frac{M_B.Math.M_B.Math.}{P_B}+\frac{M_2.Math.M_2.Math.}{P_2}\right)=0& \left(18\right)\end{array}$ $\begin{array}{cc}{M}_L=A_e{P}_d\sqrt{\frac{k}{RT}{\left(\frac{2}{k+1}\right)}^{\frac{k+1}{k-1}}}& \left(19\right)\end{array}$ $\begin{array}{cc}{M}_1-M_2=M_L& \left(20\right)\end{array}$
P.sub.1=P.sub.2=P.sub.L  (21) where P.sub.1 is a fluid pressure, whose dimensional unit is Pa, at the Δx distance before flowing into the leakage point, M.sub.1 is a flow rate, whose dimensional unit is kg.Math.s.sup.−1, at the Δx distance before flowing into the leakage point; P.sub.2 is a fluid pressure, whose dimensional unit is Pa, at the Δx distance after flowing out of the leakage point, and M.sub.2 is a flow rate, whose dimensional unit is kg.Math.s.sup.−1, at the Δx distance after flowing out of the leakage point; and P.sub.L is a pressure at the leakage point, whose dimensional unit is Pa.

5. The urban non-metallic pipeline leakage location method according to claim 1, wherein, taking a minimum value of the difference between the calculated value and the measured value of the pressure at the monitoring point as the control target to define the objective function in step S4 is specifically as follows: S4.1: defining an objective function with a least squares criterion: $\begin{array}{cc}E=\min \sqrt{\underset{i=1}{\overset{M}{.Math.}}{\left(P_i-P_i^{\prime }\right)}^2}& \left(22\right)\end{array}$ where E is an objective function; M is a number of time steps; P.sub.i is a calculated pressure, Pa; P′.sub.i is a measured pressure, Pa; and the effective leakage area that best conforms to the actual situation is generated by minimizing the objective function E, and whether a node is leaking is determined according to a value of the effective area; S4.2: determining a limit range of the effective leakage area:
0≤A.sub.ei≤A.sub.emax  (23) where A.sub.ei , is an effective leakage area at node i; and A.sub.emax is a maximum limit of the leakage area, which is determined as a reasonable proportion of a cross sectional area of the pipeline; S4.3: calculating a friction resistance coefficient by a friction resistance calculation formula according to an experimental data, wherein a flow of the gas in the pipeline is generally in a;turbulent region, specifically: Blasius formula $\begin{array}{cc}\lambda =\frac{0.3164}{{\mathrm{Re}}^{0.25}}& \left(24\right)\end{array}$ wherein Reynolds coefficient calculation formula adopts an empirical formula: $\begin{array}{cc}\mathrm{Re}=0.354\times \frac{Q_v}{D\upsilon }& \left(25\right)\\ \upsilon =\frac{\mu }{\rho }& \left(26\right)\end{array}$ where λ is a friction resistance coefficient; Re is Reynolds coefficient; Q.sub.νis a volume flow rate, whose dimensional unit is m.sup.3.Math.s.sup.−1; D is an inner diameter of the pipeline, whose dimensional unit is m; ν is a kinematic viscosity, whose dimensional unit is m.sup.2.Math.s.sup.−1; μ is a dynamic viscosity, whose dimensional unit is Pa.Math.s; and ρ is a gas density, whose dimensional unit is kg.Math.m.sup.−3.

6. The urban non-metallic pipeline leakage location method according to claim 1, wherein, the SQP algorithm is used in step S5 to converge the objective function E within the flow constraint ranges at the both ends of the monitoring point set in S4, and performing the algorithm optimization comprises specifically as follows: S5.1: establishing a file and defining an objective function E(x) according to experimental data; S5.2: establishing a file and defining a nonlinear constraint in a constraint condition; S5.3: determining an initial value X.sub.0 of iteration; S5.4: determining upper and lower limits VLB and VUB of variable; and S5.5: establishing a main program comprising a function of a nonlinear programming solution, and running the solution.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The present invention is further descried blow in conjunction with the accompanying drawings and embodiments.

(2) FIG. 1 is a schematic flowchart according to the present invention.

(3) FIG. 2 is a schematic structural diagram of a pipeline system according to the present invention.

(4) FIG. 3 is a schematic diagram of a characteristic mesh with a specific time interval.

(5) FIG. 4 is a schematic diagram of leakage mesh division.

(6) FIG. 5 is a simulated pipeline pressure profile in two states.

(7) FIG. 6 is a simulated pipeline velocity vector profile in two states.

(8) FIG. 7 is a diagram of pipeline experiment data collection.

(9) In the drawings: 1—upstream flow sensor, 2—upstream pressure sensor, 3—upstream infrasonic sensor, 4—leakage valve, 5—downstream infrasonic sensor, 6—downstream pressure sensor, 7—downstream flow sensor.

DETAILED DESCRIPTION

(10) The present invention is now descried in detail in conjunction with the accompanying drawings. The drawings are simplified schematic diagrams, which only illustrate the basic structure of the present invention in a schematic way, so the drawings only show the structures related to the present invention.

(11) As shown in FIG. 2, the experimental pipe section of the present invention adopts three straight pipe sections connected in sequence to form a U-shaped structure, wherein X.sub.1-X.sub.2 pipe section is 3.2 m long, X.sub.2-X.sub.3 pipe section is 1.6 m long, and X.sub.3-X.sub.4 pipe section is 3.2 m long. The leakage hole is at X.sub.5. At 2 m from the pipe inlet, the diameter of the leakage hole is 1 mm. X.sub.1 is positioned at the inlet of the pipeline, X.sub.4 is positioned at the outlet of the pipeline, and X.sub.6 is 2.8 m away from the inlet. Moreover, an upstream flow sensor 1, an upstream pressure sensor 2 and an upstream infrasonic sensor 3 are arranged at the inlet end of the pipeline. A downstream flow sensor 7, a downstream pressure sensor 6 and a downstream infrasonic sensor 5 are sequentially arranged at the outlet end of the pipeline. The leakage valve 4 is provided at X.sub.3.

Embodiment 1

(12) in this embodiment, Fluent software is used for simulating law analysis to perform leakage determination.

(13) As shown in FIG. 1, the urban non-metallic pipeline leakage location method according to the present invention includes: using Fluent software to simulate pressure and flow rate parameter distributions of a gas in non-metallic pipeline leakage and non-leakage states, studying the characteristics of gas leakage in a non-metallic pipeline to give standards for determining whether leakage occurs in the pipeline, and comparing the experimental data with the standards to determine leaking pipeline and non-leaking pipelines, then locating leaking points of the leaking pipelines in combination with an inverse-transient method.

(14) S1: Determine whether leakage occurs in a pipeline.

(15) Step A1: model an experimental U-shaped pipeline. The size settings are the same as those of the experimental device. Two end pipe sections are 3.2 m long, the middle pipe section is 1.6 m long, and the pipe diameter is 0.0456 m. The leakage hole is positioned 2 m distant from the head end of the pipeline, and the diameter of the leakage hole is 1 mm. A structured tetrahedral mesh is used for mesh division to obtain divided regions, and proper mesh refinement is performed near a leakage port. See FIG. 4 for the mesh refinement at the leakage port.

(16) Step A2: control equation and algorithm selection:

(17) choosing a standard k−ε double equation model as a numerical model; in the model, using turbulent energy k to reflect the characteristic velocity, using a turbulent energy dissipation rate ε to reflect the characteristic length scale, using Boussinesq hypothesis for simplification to form k equation and ε equation, and combining the two equations together to form a closed equation set to simulate the flow, specifically:

(18) $\begin{array}{cc}k\mathrm{equation}\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho {\mathrm{ku}}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\mathrm{\rho \epsilon }-Y_M+S_K& \left(1\right)\end{array}$$\begin{array}{cc}\epsilon \mathrm{equation}\frac{\partial \left(\mathrm{\rho \epsilon }\right)}{\partial t}+\frac{\partial \left(\rho {\mathrm{ku}}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_t}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\frac{C_{1\epsilon }\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}+S_{\epsilon }& \left(2\right)\end{array}$

(19) where x.sub.i, x.sub.j, and x.sub.k are coordinate components; u is the velocity vector; ρ is the fluid density; μ is the viscosity coefficient; μ.sub.i is the turbulent viscosity; G.sub.k and G.sub.b are production terms of velocity gradient and turbulent kinetic energy k caused by buoyancy; constants C.sub.1ε=1.44, C.sub.2ε=1.92, C.sub.3ε=0.09, σ.sub.k=1.0, and σ.sub.t=1.3; S.sub.k and S.sub.ε are user-defined source items; and Y.sub.M is a pulsating expansion term in compressible turbulent flow;

(20) choosing a SIMPLE algorithm as the numerical calculation algorithm of the model; for a given pressure, solving a momentum equation in a discrete form to obtain a velocity field, and then correcting the given pressure and iterating repeatedly until a convergent solution is obtained.

(21) Step A3: boundary condition settings:

(22) Inlet boundary condition settings: the inlet boundary condition for a compressible fluid should be set as pressure inlet; in this experiment, numerical simulations are performed under three inlet pressures, which are set to 0.3 MPa, 0.2 MPa and 0.1 MPa, respectively.

(23) outlet boundary condition settings: the outlet boundary condition is set as a pressure outlet, and both the end outlet of the pipeline and the end outlet of the leakage hole are set to the external environmental pressure;

(24) fluid medium settings: the medium in the air is air, and the physical properties of the air are shown in Table a;

(25) pipeline material settings: the material of the pipeline is PE, the specific heat capacity is 0.45 cal/g° C., the density is 0.956 g/cm.sup.3, and the thermal conductivity is 0.5 W/m.Math.K.

(26) TABLE-US-00001 TABLE a Table of physical properties of air at room temperature specific heat capacity Prandtl density kg/m.sup.3 viscosity 10.sup.−5 Pa .Math. s kJ/(kg .Math. K) number Pr 1.225 1.81 1.005 0.703

(27) step A4: simulate the changes in internal pressure and flow rate parameters of the pipeline in the pipeline leakage and non-leakage states respectively under given multiple sets of inlet pressures to find out the change laws of pressure and flow rate parameters, so as to give standards for determining whether leakage occurs in the pipeline, and then compare the actually measured pipeline data with the standards to determine leaking pipelines and non-leaking pipelines.

(28) step A4.1: numerical simulation of flow field distribution of pipeline pressure and flow velocity.

(29) Referring to figure a in FIG. 5, when there is no leakage under the inlet pressure of 0.3 MPa, the pressure values in the pipeline is generally linearly distributed from upstream to downstream, and the pressure gradient within a certain range is maintained to decrease with the increase of pipeline length. Referring to figure b in FIG. 5, under the same inlet pressure, when leakage occurs in the pipeline, the overall pressure in the pipeline decreases significantly, the most obvious is that the pressure in the vicinity of the leakage hole increases about 4.5×10.sup.4 Pa as compared with that under non-leakage, and the pressure gradients in front of and behind the leakage hole decrease significantly. Referring to figures b, c, and d in FIG. 5, by comparing the pressure profiles of the pipeline leakage state under different inlet pressures, the greater the pipeline inlet pressure is, the greater the pressure gradient in the pipeline is, and the pressure in the vicinity of the leakage hole generally decreases and the pressure change range of the pipe sections in front of and behind the leakage hole decreases.

(30) Referring to figure a in FIG. 6, when there is no leakage at 0.3 MPa, the flow velocity in the pipeline tends to be stable, and the overall change range is very small. Referring to figure b in FIG. 6, when there is leakage at 0.3 MPa, the flow velocity on the outer wall of the pipeline increases significantly, because when the pipeline is leaking, the pressure difference between the inside and outside of the pipeline accelerates the flow of the gas in the pipeline, resulting in a larger difference between the flow velocity outside the pipeline and the flow velocity inside. Referring to figures b, c, and d in FIG. 6, by comparing the velocity vector profiles of the pipeline leakage state under different inlet pressures, it can be known that in the leakage state, the greater the inlet pressure of the pipeline is, the greater the overall flow velocity in the pipeline is, and the greater the difference between flow velocities on the inner and outer walls are.

(31) Step A4.2: flow field distribution of pressure and flow velocity of an experimental pipeline:

(32) Referring to FIG. 7, the broken line based on square mark indicates the upstream flow rate, and the broken line based on circular mark is the downstream flow rate, the broken line based on regular angular mark is the upstream pressure, and the broken line based on inverted triangular mark is the downstream pressure. As shown in figures a, b, and c in FIG. 7, when leakage occurs, the flow rate at the front end of the leakage hole greatly increases, and the flow rate at the rear end of the leakage hole slightly decreases. After a period of time, the flow rate in the front section of the leakage hole decreases, and the flow rates at both ends of the leakage hole remain unchanged. The greater the inlet pressure is, the greater the difference in flow rate is. When leakage occurs, the pressure of the entire pipeline decreases, and the difference between upstream and downstream pressures is not obvious.

(33) Step A4.3: comparative analysis of numerical simulation and experimental results

(34) In the experiment, through simulation and experimental comparison on pipeline state changes before and after leakage due to holes with difference diameters of a simulated urban fuel pipeline, substantially similar conclusions can be obtained, but there are also different results. Similar conclusions: when leakage occurs, the overall pressure in the pipeline decreases, the pressure gradient decreases, and the gas flow rate in the pipeline becomes faster. Different conclusions: in the experiment, the flow velocity in the front pipe section of the leakage hole is accelerated, and the flow velocity in the back pipe section is slightly reduced, and in the simulation, the flow velocities on the inner and outer walls of the entire pipeline are quite different, but the flow velocities in front of and behind the leakage hole are not changed significantly. Through comprehensive numerical simulation and experimental analysis, and in combination with related articles, two standards for parameter changes during the occurrence of leakage in the pipeline are obtained: (1) when leakage occurs, the overall pressure in the pipeline decreases, the pressure gradient decreases, and the gas flow rate in the pipeline becomes faster; (2) the flow velocity in the front section of the leakage hole increases, the flow velocity in the back section decreases slightly, and the flow rate differences in front of and behind the leakage hole gradually increase in a short time and then remain unchanged.

(35) S2: pipeline leakage experiment: obtaining pipeline pressure and flow rate data.

(36) Transient pipe flow is produced by gradually closing a valve, the valve closing time lasts long enough to minimize the influence of instability and the uncertainty of a friction coefficient, so as to obtain pressure and flow rate data of measuring points X1, X2, X3, and X4 in different periods of time, and pressure and flow rate data at Δx, i.e., 1.6 m, from both ends of the measuring points, see Tables 1 and 2.

(37) TABLE-US-00002 TABLE 1 Flow rate data of each node under leakage at 0.3 Mpa Sequence X1 flow rate/ X2 flow rate/ X3 flow rate/ X4 flow rate/ number m3/h m3/h m3/h m3/h 1 2.82 2.79 2.79 2.79 2 3.07 2.74 2.74 2.73 3 3.7 2.67 2.67 2.66 4 4.79 2.58 2.58 2.56 5 6.62 2.49 2.49 2.48 6 7.16 2.49 2.49 2.48 7 7.66 2.49 2.49 2.48 8 8.22 2.49 2.48 2.48 9 8.66 2.47 2.47 2.47 10 8.66 2.46 2.46 2.46 11 9.01 2.49 2.49 2.48 12 9.22 2.49 2.49 2.48 13 9.42 2.49 2.49 2.48 14 9.54 2.49 2.49 2.48 15 9.64 2.49 2.49 2.48

(38) TABLE-US-00003 TABLE 2 Pressure data of each node under leakage at 0.3 Mpa Sequence X1 X2 X3 X4 number pressure/MPa pressure/MPa pressure/MPa pressure/MPa 1 0.26 0.26 0.26 0.25 2 0.24 0.24 0.24 0.23 3 0.22 0.22 0.22 0.21 4 0.2 0.2 0.2 0.19 5 0.19 0.19 0.19 0.17 6 0.18 0.18 0.18 0.17 7 0.17 0.17 0.17 0.16 8 0.16 0.16 0.16 0.16 9 0.16 0.16 0.16 0.15 10 0.15 0.15 0.15 0.15 11 0.15 0.15 0.15 0.14 12 0.15 0.15 0.15 0.14 13 0.14 0.14 0.14 0.13 14 0.14 0.14 0.14 0.13 15 0.14 0.14 0.14 0.13

(39) S3: Establish a pipeline leakage control equation.

(40) Assuming that the node is a leakage point, a control equation of a non-metallic pipeline gas leakage state is established by using a nonlinear equation solved based on continuity and motion equations of a gas and a leakage amount equation of a leakage hole, for example:

(41) $\begin{array}{cc}{P}_1-P_A+\frac{a}{A}\left(M_1-M_A\right)+\frac{\lambda a^2\Delta x}{4D{A}^2}\left(\frac{M_A.Math.M_A.Math.}{P_A}+\frac{M_1.Math.M_1.Math.}{P_1}\right)=0& \left(17\right)\end{array}$$\begin{array}{cc}{P}_B-P_2+\frac{a}{A}\left(M_2-M_B\right)+\frac{\lambda a^2\Delta x}{4D{A}^2}\left(\frac{M_B.Math.M_B.Math.}{P_B}+\frac{M_2.Math.M_2.Math.}{P_2}\right)=0& \left(18\right)\end{array}$$\begin{array}{cc}{M}_L=A_e{P}_d\sqrt{\frac{k}{RT}{\left(\frac{2}{k+1}\right)}^{\frac{k+1}{k-1}}}& \left(19\right)\end{array}$$\begin{array}{cc}{M}_1-M_2=M_L& \left(20\right)\end{array}$$\begin{array}{cc}{P}_1=P_2=P_L& \left(21\right)\end{array}$

(42) where P.sub.A and M.sub.A are the fluid pressure (Pa) and flow rate (kg.Math.s.sup.−1), at Δx in front of the node; P.sub.B and M.sub.B are the fluid pressure (Pa) and flow rate (kg.Math.s.sup.−1), at Δx behind the node; and P.sub.L is the pressure at the leakage point (PaA); P.sub.1 and M.sub.1 are the fluid pressure (Pa) and flow rate (kg.Math.s.sup.−1) , at a distance before flowing into the leakage point; P.sub.2 and M.sub.2 are the fluid pressure (Pa) and flow rate (kg.Math.s.sup.−1), at a distance after flowing out of the leakage point; P.sub.L is the pressure (Pa) at the leakage point; A.sub.e=CA.sub.0 is the effective leakage area (m.sup.2); C is the orifice coefficient, which is related to the shape of the leakage hole; A.sub.0 is the area of the leakage hole (m.sup.2); P.sub.d is the initial pressure (Pa) of the pipeline; k is the adiabatic coefficient of the gas, dimensionless; R is the gas constant, J.Math.kg.sup.−1.Math.K.sup.−1; and T is the temperature (K) of the gas; D is the inner diameter (m) of the pipeline; λ is the friction resistance coefficient; and A is the cross sectional area (m.sup.2) of the pipeline.

(43) S4: Define an objective function:

(44) S4.1: show the difference between the calculated pressure and the real pressure of the node according to the leakage control equations (17-21), for example:

(45) $\begin{array}{cc}\mathrm{Min}E=.Math.\frac{\frac{\lambda a^2\Delta x\left(M_1^2+M_2^2\right)}{4{\mathrm{DA}}^2}}{\begin{array}{c}{P}_A-P_B+\frac{a}{A}\left(M_A+M_B\right)-\\ \frac{a}{A}\left(M_1+M_2\right)-\frac{\lambda a^2\Delta x}{4{\mathrm{DA}}^2}\left(\frac{M_A^2}{P_A}+\frac{M_B^2}{P_E}\right)\end{array}}-P_{\mathrm{real}}.Math.& \left(22\right)\end{array}$

(46) S4.2: Determine the effective range of the leakage area, so as to determine the upper limit of the leakage amount. Herein, the effective leakage area when the radius of the leakage hole is 5 mm is taken as the maximum limit of the leakage area.
0≤A.sub.e≤7.85×10.sup.−5m.sup.2  (23)

(47) S4.3: Determine the friction resistance coefficient in the pipeline.

(48) According to the different flow states of the gas in the PE pipe, different calculation formulas are used for the friction resistance coefficients λ of different pipe materials. A PE pipe calculation formula used in this article is the Blasius formula.

(49) Blasius formula:

(50) $\begin{array}{cc}\lambda =\frac{0.3164}{{\mathrm{Re}}^{0.25}}& \left(24\right)\end{array}$

(51) wherein the Reynolds coefficient calculation formula adopts an empirical formula:

(52) $\begin{array}{cc}\mathrm{Re}=0.354\times \frac{Q_v}{D\upsilon }& \left(25\right)\end{array}$$\begin{array}{cc}\upsilon =\frac{\mu }{\rho }& \left(26\right)\end{array}$

(53) where λ is the friction resistance coefficient; Re is the Reynolds coefficient; Q.sub.ν is the volume flow rate (m.sup.3/s); D is the inner diameter (m) of the pipeline; ν is the kinematic viscosity (m.sup.2/s); μ is the dynamic viscosity (Pa.Math.s); and ρ is the gas density (kg/m.sup.3). At 24 degrees Celsius, the dynamic viscosity of air is 1.83×10.sup.−5 Pa.Math.s, and the kinematic viscosity is 1.4364×10.sup.−5m.sup.2/s. By substituting the steady-state data and parameters under non-leakage into equations (24-26), it is calculated that the friction resistance coefficient in the pipeline is approximately 0.01.

(54) S5: Algorithm optimization.

(55) S5.1: Substitute the pressure and flow rate data of the node X2 in the 14th to 15th sets in the experimental data, the pressure and flow rate data at 1.6 m in front of and behind X2, and known parameters into equation (22) to obtain an objective function:

(56) 0$\begin{array}{cc}\mathrm{Min}E=.Math.\frac{3.80587305773\times {10}^9\left(M_1^2+M_1^2\right)}{781.91-208295.35\left(M_1+M_2\right)}-P_{\mathrm{real}}.Math.& \left(27\right)\end{array}$

(57) S5.2: Define a nonlinear constraint
M.sub.1−M.sub.2≥0  (28)
0≤M.sub.1≤M.sub.B+Δt+M.sub.Lmax  (29)
M.sub.B+Δt≤M.sub.2≤M.sub.A+Δt  (30)

(58) M.sub.Lmax is the maximum leakage amount, which is 0.0258 kg/s, calculated by substituting the leakage area of the maximum limit into the leakage amount formula (3), and the followings are obtained:

(59) $\begin{array}{cc}s.t.\left(\begin{array}{cc}-1& 1\\ 0& 0\end{array}\right)\left(\begin{array}{c}{M}_1\\ M_2\end{array}\right)\le \left(\begin{array}{c}0\\ 0\end{array}\right)& \left(31\right)\end{array}$$\begin{array}{cc}\left(\begin{array}{c}0\\ 0.00089\end{array}\right)\le \left(\begin{array}{c}{M}_1\\ M_2\end{array}\right)\le \left(\begin{array}{c}0.0293\\ 0.00089\end{array}\right)& \left(32\right)\end{array}$

(60) S5.3: Determine an initial value X.sub.0=[0.00343, 0.00089]

(61) S5.4: VLB=[0, 0.00089], and VUB=[0.0293, 0.00089]

(62) S5.5: Obtain a result, M.sub.1=0.0029, and M.sub.2=0.0009, that is, the calculated difference in flow rate between the front and back ends of node X.sub.2 is 0.002 kg/s, and it is calculated that the leakage area is 6.09×10.sup.−7 m.sup.2. According to steps S5.1 to S5.5, the differences in flow rate in front of and behind the X.sub.1, X.sub.3, and X.sub.4 nodes are correspondingly calculated as 0 kg/s, 0 kg/s, and 0 kg/s respectively, and the leakage areas are calculated as 0 m.sup.2, 0 m.sup.2, and 0 m.sup.2 respectively.

(63) S6: Leakage location:

(64) S6.1: After calculation, the difference in flow rate in front of and behind nodes X.sub.1, X.sub.3, and X.sub.4 is 0, and the leakage area is 0, and therefore, none of the three nodes is leaking; obvious difference in flow rate between the front and rear ends of node X.sub.2 proves that leakage occurs at node X.sub.2 in a certain position of the pipe section from 1.6 m long front section to 1.6 m long rear section, while the calculation result of node X.sub.3 is normal, so it is determined that there is no leakage at node X.sub.3 on the pipe section from 1.6 m long front section to 1.6 m long rear section. Therefore, on the pipe section 1.6 m in front of node X.sub.2 to node X.sub.2, the node is set again for calculation.

(65) S6.2: Division of nodes X.sub.5 and X.sub.6 is continued, see FIG. 2, and the leakage areas are respectively calculated to 6.09×10.sup.−7 m.sup.2 and 0 m.sup.2 according to the differences in flow rate in front of and behind nodes X.sub.5 and X.sub.6, which are 0.002 kg/s and 0 kg/s respectively, are calculated correspondingly when Δx is 0.4 m according to steps S5.1 to S5.4. It is proved that there is leakage in the pipe section at 0.4 m in front of and behind X.sub.5, and the node can be set up again according the needs to continue the calculation. In this experiment, node X.sub.5 is a leakage point, and it can be proved that the leaking pipe section can be found by calculating the pressure and flow rate of the monitoring node in this method. The shortcomings are that this method relies too much on the selection of nodes, and the establishment of the model, the calculation of parameters, etc. all have influence on the accuracy of the results.

Embodiment 2

(66) This embodiment is different from Embodiment 1 in that the method for determining pipeline leakage in step S1 is different. In this embodiment, a Markov chain-based method is adopted to realize non-metallic pipeline leakage determination.

(67) S1: Determine whether leakage occurs in the pipeline. The non-metallic pipeline leakage determination based on Markov chain includes the following steps:

(68) Step B1: collect a first set of five flow rates by using a flow sensor installed on the pipeline: 2.78, 2.79, 2.82, 3.07, and 3.7, in units of m.sup.3/h, and calculating a first set of five flow change rates, q.sub.i=ΔQ.sub.i/Q.sub.i×100%, according to the formula q.sub.1=0%, q.sub.2=3.58%, q.sub.3=1.06%, q.sub.4=8.14%, q.sub.5=17.02%.

(69) step B2: set four change states: a first state: the value of the flow change rate q.sub.i is 0, indicating that the pipeline is in a non-leakage state; a second state: the value of the flow change rate q.sub.i is 0-1(%), indicating that the pipeline is in a small leakage state; a third state: the value of the flow change rate q.sub.i is 1-3(%), indicating that the pipeline is in a leakage expanded state; and a fourth state: the value of the flow change rate q.sub.i is 3-100(%), indicating that the pipeline is in a large leakage state;

(70) setting conditions: the leakage of the pipeline is a sudden change process, the change of the leakage state is a gradual process; due to the influence of various other factors, the working state of the pipeline may be directly changed from an operation state to a certain degree of leakage in the leakage state due to the influence of various other factors, and the leakage state of the pipeline is irreversible, but the degree of leakage is reversible; according to the above-mentioned conditions, an adjacency matrix A of pipeline leakage is obtained as follows:

(71) $A=\left[\begin{array}{ccc}{A}_{11}& .Math.& A_{1j}\\ .Math.& \ddots & .Math.\\ A_{i1}& .Math.& A_{\mathrm{ij}}\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& 1& 1& 1\\ 0& 1& 1& 1\\ 0& 1& 1& 1\end{array}\right]$

(72) where A.sub.ij=1=1, 2 . . . 4; j=1, 2 . . . 4) indicates that there is a transition relationship between the i-th state and the j-th state, and A.sub.ij=0 indicates that there is no transition relationship between the i-th state and the j-th state.

(73) State classification is performed on the first set of flow change rates obtained in step 1 to classify the flow change rates into these four states. In order of time, these are counted in the five periods of time in chronological order. Each element in the state transition probability matrix P is formed according to the state transition situation of the flow change rates.

(74) The calculation process of the state transition probability matrix P is as follows:

(75) step B2.1: according to the classification, deriving that the number of flow rate changes in the first state is 1, the number of transitions from the first state to the first state is 0, the number of transitions from the first state to the second state is 0, the number of transitions from the first state to the third state is 0, and the number of transitions from the fourth state to the second state is 0, and collecting statistics about the state transition situation of these five flow change rates at the next moment:

(76) the transition probability from the first state to the first state is p.sub.11=0;

(77) the transition probability from the first state to the second state is p.sub.12=0;

(78) the transition probability from the first state to the third state is p.sub.13=1; and

(79) the transition probability from the first state to the fourth state is p.sub.14=0;

(80) step B2.2: deriving that the number of flow rate changes in the second state is 0, the number of transitions from the second state to the second state is 0, the number of transitions from the second state to the third state is 0, and the number of transitions from the second state to the fourth state is 0, and collecting statistics about the state transition situation of the zero flow change rate at the next moment:

(81) the transition probability from the second state to the first state is p.sub.21=0;

(82) the transition probability from the second state to the second state is p.sub.22=0;

(83) the transition probability from the second state to the third state is p.sub.23=0; and

(84) the transition probability from the second state to the fourth state is p.sub.24=0;

(85) step B2.3: deriving that the number of flow rate changes in the third state is 1, the number of transitions from the third state to the second state is 0, the number of transitions from the third state to the third state is 0, and the number of transitions from the third state to the fourth state is 1, and collecting statistics about the state transition situation of the zero flow change rate at the next moment:

(86) the transition probability from the third state to the first state is p.sub.31=0;

(87) the transition probability from the third state to the second state is p.sub.32=0;

(88) the transition probability from the third state to the third state is p.sub.33=0; and

(89) the transition probability from the third state to the fourth state is p.sub.34=1;

(90) step B2.4: deriving that the number of flow rate changes in the fourth state is 3, the number of transitions from the fourth state to the second state is 0, the number of transitions from the fourth state to the third state is 0, and the number of transitions from the fourth state to the fourth state is 2, and collecting statistics about the state transition situation of the zero flow change rate at the next moment:

(91) the transition probability from the fourth state to the first state is p.sub.41=0;

(92) the transition probability from the fourth state to the second state is p.sub.42=0;

(93) the transition probability from the fourth state to the third state is p.sub.43=⅓; and

(94) the transition probability from the fourth state to the fourth state is p.sub.44=⅔;

(95) then the state transition probability matrix P in the Markov chain is:

(96) $P=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1/3& 2/3\end{array}\right]$

(97) Step B3: collect a second set of flow rates: 4.79, 6.62, 7.16, 7.66, and 8.22, in units of m3/h. A second set of five flow change rates are obtained: q.sub.1=22.75, q.sub.2=27.64, q.sub.3=7.54%, q.sub.4=6.52%, q.sub.5=6.81%; like step 2, the number of these five flow change rates in each interval is counted, and the actual flow change rate in each interval at this time is edited as the first parameter vector and is set as α.sub.1=[0 0 0 5], wherein α.sub.1 represents the first parameter vector; 0 represents the number of flow change rates in the first state among the five flow rate change rates under the first parameter vector; 0 represents the number of flow change rates in the second state among the five flow rate change rates under the first parameter vector; 0 represents the number of flow change rates in the third state among the five flow rate change rates under the first parameter vector; and 5 represents the number of flow change rates in the fourth state among the five flow rate change rates under the first parameter vector;

(98) the parameter vector β.sub.1 of the predicted flow change rate based on the Markov chain in the future time period is obtained through Markov chain calculation, and the formula thereof is:
β.sub.1=α.sub.1×P, the obtained β.sub.1=[0 0 0 5].

(99) It can be derived from β.sub.1=[0 0 0 5] that, at this time, it is predicted that the number of parameters in the fourth state at the next moment is the largest, containing five parameters, wherein S=4, then at this time, it is predicted that the state at the next moment is the fourth state, i.e., the large leakage state.

(100) Step B4: collect a third set of flow rates: 8.66, 8.66, 9.01, 9.22, and 9.42, in units of m.sup.3/h, and obtain a third set of flow change rate values q.sub.1=5.08%, q.sub.2=0%, q.sub.3=3.88%, q.sub.4=2.33%, q.sub.5=2.12%. Like steps, the number of matches in each interval is counted and is set as α.sub.2=[1 021], wherein it is derived that the number of parameters in the third state at this time is the largest, containing two parameters, i.e., N.sub.W.sub.max=3, where W=2, that is, the actual state is the third state at this time.

(101) Step B5: Compare the number values of S in N.sub.S.sub.max and Win N.sub.W.sub.max.

(102) Because S=4, W=2, S>W, it is derived that the pipeline is in the late stage of the beginning of leakage, and it is determined that the state is the third state.

(103) Step B6: set the second set of data as the first set of data and the third set of data as the second set, of data, and continue to predict and determine the real-time state or leakage situation of the pipeline. In this case, the pipeline is in the third state of leakage, it is determined that leakage occurs, an alarm program is triggered, and the system starts to given an alarm.

(104) The calculation of leakage location in this example, i.e., steps S2 to S6, is the same as that in Embodiment 1. It is not necessary to repeat the description.

(105) Summary: both methods for pipeline leakage determination based on software simulation and Markov method have their advantages and disadvantages, the software simulation method can only summarize the law of smooth distribution of leaking and non-leaking states, the advantage is that the analysis of the flow field distribution laws in two states can provide theoretical ideas for the determination of leakage and give a theoretical basis; however, in specific practice, it is difficult to quantify the specific value to determine whether tiny leakage occurs in the pipeline; moreover, the Markov method can determine the leakage situation of the pipeline according to the capture of flow rate parameters, but it may cause misjudgments due to the influence of other factors on the flow rate, and requires a lot of data to ensure the accuracy of the leakage determination.

(106) Taking the above-mentioned ideal embodiment according to the present invention as enlightenment, through the above-mentioned description content, relevant workers can make various changes and modifications without departing from the scope of the present invention. The technical scope of the present invention is not limited to the content on the description, and its technical scope must be determined according to the scope of the claims.