Research on the hydraulic parameters and leak rate of natural gas flow in natural gas pipelines based on the methods of characteristic line and successive approximation

Abstract

Pressure and flow rate are the most important hydraulic parameters in natural gas pipeline flow, and leak rate is the most crucial parameter after a leak accident occurs. The study of these parameters is vital to the safey of natural gas pipelines and risk assessment after accidents. In this research, based on the conservation laws universally applicable to the motion of objects, we establish the fundamental control equation group for natural gas flow. Then, for both normal and leak conditions of pipelines, we use the characteristic line method to derive the corresponding difference equations for the fundamental control equations, thereby providing calculation methods for pressure and flow rate. Finally, we investigate the calculation method of natural gas pipeline leak rate under different leak aperture sizes, and validate the accuracy of this method through simulation examples.

Keywords

Hydraulic parameters method of characteristic line normal and leak conditions leak rate natural gas pipelines

1. Introduction

In today’s rapidly developing world, with a high demand for energy and alarming levels of greenhouse gas emissions leading to severe environmental pollution, economic sustainability and human survival face enormous threats. The severity of this situation has led countries worldwide to place great emphasis on energy conservation and environmental protection. Natural gas, as a new energy source of the 21st century, has many advantages such as high calorific value, low environmental pollution, abundant reserves, and high economic benefits. It plays an vital role in improving the energy structure, alleviating energy supply-demand conflicts, and enhancing environmental quality. Natural gas has already become the third-largest energy source following coal and oil, and is hailed as the environmentally friendly energy source of the 21st century. Natural gas pipeline transportation is a complex and sophisticated system. The safety of this system lies on several factors, which including the design and construction of the pipelines, the quality of the materials used, the monitoring and maintenance of the pipelines, and the training and competence of the personnel involved. In recent years, there have been several incidents of natural gas pipeline leaks and explosions that have caused significant damage and loss of life. These incidents have highlighted the importance of safety in natural gas pipeline transportation and the need for rigorous safety measures to prevent and mitigate pipeline failures.

After a pipeline leakage accident, the importance of pipeline leakage for safety should not be ignored. The leaked natural gas can cause a range of problems, including fire, explosion, and environmental pollution, among others. Therefore, it is critical to have effective detection and mitigation systems in place to address pipeline leaks promptly. The importance of natural gas pipeline leakage for safety also extends to the response and recovery efforts after a leakage incident. The response should be swift, coordinated, and efficien. Therefore, studying the natural gas flow in pipelines under normal and leak conditions is of important practical significance for accident prevention and emergency response.

Through the literature search, we have identified several scholars who have made significant contributions to the study of hydraulic parameters and natural gas pipeline flow as well as the detection and estimation of pipeline leak rates. For instance, in their work, Smith et al. (2014) invested a hydraulic model for predicting the flow behavior of natural gas in pipelines [1], while Jones et al. (2017) focused on investigating advanced sensing technologies for pipeline leak detection and prevention [2]. Other notable researchers in this field include Brown et al. (2015), who developed several research on the impact of flow rate and viscosity on hydraulic parameters in natural gas pipeline flow [3], and Wang et al. (2018), who proposed a novel method for optimizing pipeline performance using the method of characteristic line [4]. Chen et al. (2015) put forward a method for leak detection of gas pipeline employing the method of characteristic line [5]. Similarly, Liu et al. (2016) utilized the method of characteristic line to calculate the leakage rate [6]. In their study, Li et al. (2018) proposed a method for detecting gas pipeline leaks by combining the wavelet transform and the method of characteristic line [7]. Zhou et al. (2019) suggested a technique for calculating the leakage rate using the method of characteristic line and machine learning approaches [8]. Furthermore, Yin et al. (2020) introduced a transient temperature and pressure variations based leak detection method, employing the method of characteristic line in conjunction with artificial neural networks [9]. Lastly, Wang et al. (2021) presented a leak detection method which based on spectral analysis of temperature and pressure variations using the method of characteristic line [10] Montiel et al. [11] initially presented the hole-pipe model and the application of various models for detecting leakage. Zhou et al. [12] performed an analysis of the gradual discharge of gas storage tanks, focusing on the thermal process, and developed a corresponding mathematical model. Levenspiel [13] formulated a leakage model specifically tailored for complete pipeline fractures. Their investigation involved analyzing the pressure descend along the pipeline, assuming constant pressure at the beginning point and adiabatic natural gas flow and leak processes within the pipe. Additionally, Dong et al. [14] calculated the steady gas leak rate in pipelines. Furthermore, several researchers have explored pipeline health monitoring and leak detection [15, 16, 17, 18, 19, 20, 21].

The above research reveals that the majority of studies on natural gas pipeline leak rate calculations are based on the research findings proposed by Montiel et al. [11], lacking new perspectives and innovations.

2. Fundamental control equations for natural gas flow

The dominating equations of control for the flow of natural gas encompass the continuity equation, momentum equation, energy equation, and the gas state equation when dealing with gases. These fundamental control equations are applicable to both normal and leakage conditions in natural gas pipelines. To establish the basic control equations for natural gas flow in pipelines, certain reasonable simplifications can be employed, based on the following three assumptions [11]:

(a)
One-dimensional gas flow within the pipeline;
(b)
The pipeline is rigid;
(c)
The slope of a specific section of the natural gas pipeline remains constant.

Under these above assumptions, the following set of basic control equations can be derived:

(1) Continuity equation

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u)}{\partial x% }=0$ (1)

where $t$ is the time ( $s$ ), $\rho$ is the natural gas density (kg/m³), $x$ is the distance from the beginning of the pipe ( $m$ ), $u$ is the flow velocity ( $m/s$ ).

(2) Momentum equation

$\displaystyle\frac{\partial(\rho u)}{\partial t}+\frac{\partial(\rho u^{2})}{% \partial x}+\frac{\partial P}{\partial x}+g\rho\sin\theta+\frac{\lambda}{D}% \frac{u^{2}}{2}\rho=0$ (2)

where $\theta$ is the angle of pipeline inclination (rad), $g$ is the gravitational accelevaluen ( $m/s^{2}$ ), $D$ is the diameter of pipe ( $m$ ), $P$ is the pressure of natural gas (Pa) $\lambda$ is the hydraulic friction coefficient.

(3) Energy equation

$\displaystyle\frac{\partial}{\partial t}\left[{\rho\left(e+\frac{u^{2}}{2}+gz% \right)}\right]+\frac{\partial}{\partial x}\left[{\rho u\left(h+\frac{u^{2}}{2% }+gz\right)}\right]+\rho u\frac{\partial q}{\partial x}=0$ (3)

where $q$ is the amount of heat emitted per unit flow (J/kg), $e$ is the internal energy of natural gas (J/kg), $z$ is the height of the natural gas pipeline ( $m$ ), $h$ is the enthalpy of natural gas (J/kg).

(4) Gas state equation

$\displaystyle\rho=\frac{P}{ZR_{m}T}$ (4)

where $T$ is the temperature ( $K$ ), $R_{m}$ is gas consatant (J/(kg.K)), $Z$ is the gas compression factor.

The flow state in a pipe can be divided into two categories: transient flow and steady-state flow. Transient flow refers to the flow parameters at any point along the pipeline not only being dependent on the location but also on time. Therefore, the above basic control equations for natural gas flow (Eqs (1) to (4)) can also be referred to as the transient flow equations for natural gas pipelines. Steady-state flow, on the other hand, refers to the flow parameters at any point along the pipeline being dependent only on the location and not on time. When it is steady-state flow, the partial derivative terms with respect to time in Eqs (1) to (4) are zero.
3. Calculation method for pressure and flow rate under normal condition

Pressure and flow rate are the most crucial hydraulic parameters in natural gas pipelines. Researching their calculation methods holds significant importance in ensuring pipeline safety and enabling gas leak detection in natural gas pipelines. Equations (1) to (4) accurately represent the transient flow of natural gas pipelines, but for practical engineering applications, further simplifications can be made. In the case of long, straight pipelines, the difference of temperature between the surrounding soil and pipe’s outer wall, heat between the wall and the gas in the pipe should be considered small. As a result, we can assume that the flow within the pipeline is isothermal, then calculations made under this assumption yield acceptable results. Consequently, when temperature variations inside the pipeline are not taken into account, the energy Eq. (3) can be disregarded. Moreover, the gas flow velocity and sound speed within the pipeline are much smaller compared to the fluid velocity. Hence, the convective part in Eq. (2) can be ingored. Additionally, if the discrepancy in elevation is small, the gravity part can be ignored as well. Furthermore, the Eqs (1) and (2) are derived based on the study of fluid with unit volume. For practical convenience, the mass flow rate is introduced by utilizing the control equation for fluid with unit mass. Therefore, the following simplified equations can be obtained:

$\displaystyle\frac{\partial P}{\partial t}+\frac{c^{2}}{A}\frac{\partial Q}{% \partial x}=0$ (5) $\displaystyle\frac{\partial Q}{\partial t}+A\frac{\partial P}{\partial x}+% \frac{\lambda c^{2}}{2\textit{DA}}\frac{Q\left|Q\right|}{P}=0$ (6)

where $Q$ is the average mass flow rate (kg/s), $A$ is the cross sectional area of the pipe ( $m^{2}$ ), $c$ is sound velocity ( $m/s$ ).

Equations (5) and (6) represent hyperbolic partial differential equations that are challenging to get an analytical solution. The method of characteristics is a numerical calculation technique based on the theory of characteristics, used to solve systems of hyperbolic partial differential equations. This method is highly mature and capable of providing highly accurate numerical solutions. With the advent of electronic computers, the method has become more refined and widely applied.Therefore, this method is employed to solve these equations in this paper, the basic concept behind applying the characteristics line method to derive the partial differential Eqs (5) and (6) involves several steps. Firstly, the equations are deformed into normal differential equations by using the characteristic lines. This deformation helps simplify the problem by reducing it to ordinary differential equations. Next, the ordinary differential equations is further transformed into a system of difference equations, allowing for discrete calculations at specific points along the characteristic lines. This discretization is a crucial step in the numerical solution process. Finally, mathematical methods are utilized to convert the system of difference equations into a system of nonlinear equations. These equations, along with the initial and boundary conditions, are then solved numerically using appropriate numerical algorithms. Overall, the method of characteristics provides a systematic approach to numerically solve hyperbolic partial differential equations like Eqs (5) and (6), offering a practical and effective solution method in engineering applications.

Equations (5) and (6) can be represented as $E_{1}$ and $E_{2}$ , respectively:

$\displaystyle\left.{\begin{array}[]{l}E_{1}=\displaystyle\frac{\partial P}{% \partial t}+\frac{c^{2}}{A}\frac{\partial Q}{\partial x}=0\\ E_{2}=\displaystyle\frac{\partial Q}{\partial t}+A\frac{\partial P}{\partial x% }+\frac{\lambda c^{2}}{2DA}\frac{Q\left|Q\right|}{P}=0\\ \end{array}}\right\}$ (7)

The two equations are linearly combined using an arbitrary unknown multiplier $\delta$ as follows:

$\displaystyle E=\delta E_{1}+E_{2}=\delta\left(\frac{\partial P}{\partial t}+% \frac{A}{\delta}\frac{\partial P}{\partial x}\right)+\left(\frac{\partial Q}{% \partial t}+\delta\frac{c^{2}}{A}\frac{\partial Q}{\partial x}\right)+\frac{% \lambda c^{2}}{2\textit{DA}}\frac{Q\left|Q\right|}{P}=0$ (8)

We can select two values for $\delta$ to transform Eq. (8) into a special ordinary differential equation expressed in terms of variables $P$ and $Q$ . $P=P(x,t)$ and $Q=Q(x,t)$ be solutions of Eq. (8), then their total derivatives are:

$\displaystyle\frac{dP}{dt}=\frac{\partial P}{\partial t}+\frac{\partial P}{% \partial x}\frac{dx}{dt}$ (9) $\displaystyle\frac{dQ}{dt}=\frac{\partial Q}{\partial t}+\frac{\partial Q}{% \partial x}\frac{dx}{dt}$ (10)

If we compare the first term of Eq. (8) with Eq. (9), and the second term within the parentheses of Eq. (8) with Eq. (10), and satisfy the following conditions:

$\displaystyle\frac{dx}{dt}=\frac{A}{\delta}$ (11) $\displaystyle\frac{dx}{dt}=\delta\frac{c^{2}}{A}$ (12)

Therefore, Eq. (8) can be transformed into an ordinary differential equation:

$\displaystyle\delta\frac{dP}{dt}+\frac{dQ}{dt}+\frac{\lambda c^{2}}{2DA}\frac{% Q\left|Q\right|}{P}=0$ (13)

Equation (11) is equal to Eq. (12), that is:

$\displaystyle\frac{A}{\delta}=\delta\frac{c^{2}}{A}$ (14)

Solving the above equations, we obtain that $\delta=\pm A\mathord{\left/{\vphantom{Ac}}\right.\kern-1.2pt}c$ ï¼Œindicating that $\delta$ is two distinct real numbers. Substituting $\delta$ into Eqs (11) and (12), we have:

$\displaystyle\frac{dx}{dt}=\pm c$ (15)

Equation (15) is the two conditions that Eq. (13) must satisfy.

By substituting the values of $\delta$ into Eq. (13) and corresponding combination with Eq. (15), two equations of ordinary differential equations can be obtained:

$\displaystyle C^{+}:\left\{{\begin{array}[]{l}\displaystyle\frac{dP}{dt}+\frac% {c}{A}\frac{dQ}{dt}+\frac{\lambda c^{3}}{2DA^{2}}\frac{Q\left|Q\right|}{P}=0\\ \displaystyle\frac{dx}{dt}=+c\\ \end{array}}\right.$ (16) $\displaystyle C^{-}:\left\{{\begin{array}[]{l}\displaystyle\frac{dP}{dt}-\frac% {c}{A}\frac{dQ}{dt}-\frac{\lambda c^{3}}{2DA^{2}}\frac{Q\left|Q\right|}{P}=0\\ \displaystyle\frac{dx}{dt}=-c\\ \end{array}}\right.$ (17)

The characteristic equations $\frac{dx}{dt}=\pm c$ in the ordinary differential equation system Eqs (16) and (17) are called characteristic line equations, which have the geometric meaning of two straight lines on the $x-t$ plane. The spatiotemporal plane of the pipeline is discretized with space interval $\Delta x$ and time interval $\Delta t$ , and the characteristic line grid is established, as shown in Fig. 1. MO and NO are two characteristic lines, which integrate $C^{+}$ and $C$ along the two characteristic lines, respectively. By introducing first-order approximation, the difference equations can be geted:

$\displaystyle C^{+}:(P_{i,j}-P_{i-1,j-1})+\frac{c}{A}(Q_{i,j}-Q_{i-1,j-1})+% \frac{\lambda c^{3}\Delta t}{4DA^{2}}$ $\displaystyle\quad\left(\frac{Q_{i,j}\left|{Q_{i,j}}\right|}{P_{i,j}}+\frac{Q_% {i-1,j-1}\left|{Q_{i-1,j-1}}\right|}{P_{i-1,j-1}}\right)=0$ (18) $\displaystyle C^{-}:(P_{i,j}-P_{i+1,j-1})-\frac{c}{A}(Q_{i,j}-Q_{i+1,j-1})-% \frac{\lambda c^{3}\Delta t}{4DA^{2}}$ $\displaystyle\quad\left(\frac{Q_{i,j}\left|{Q_{i,j}}\right|}{P_{i,j}}+\frac{Q_% {i+1,j-1}\left|{Q_{i+1,j-1}}\right|}{P_{i+1,j-1}}\right)=0$ (19)

Figure 1.

Pipeline characteristic lines under normal condition.

Equations (3) and (3) represent two equations that satisfy the $O$ point at a specific time $j\Delta t$ . By combining these equations, it becomes possible to solve for the flow rate and pressure at each segmented point on the $O$ point at that time $j\Delta t$ . A similar difference equation can be established for each segmented point along the pipeline, except for the two endpoints. Therefore, calculating the flow rate and pressure from time $t=0$ to the desired time is feasible, given that the initial conditions are known.

In Fig. 1, there exists only one characteristic equation at each ending of the pipeline, which implies that the flow rate and pressure at the endpoints cannot be directly determined. Hence, it is essential to provide boundary conditions related to these two parameters at two ends of the pipe to obtain these two parameters at each time on the boundary. These boundary conditions play a crucial role in solving the differential equations accurately and comprehensively.

4. Calculation method for pressure and flow rate under leakage condition

Similar to the calculation method under normal conditions, when considering a sudden leak occurring at a specific point on the pipeline, we can still utilize the characteristics line method to solve the simplified Eqs (5) and (6). As depicted in Fig. 2, let’s assume that one sudden leak occur at node $i$ , and at time $j$ .

Figure 2.

Pipeline characteristic lines under leakage condition.

By introducing first-order approximation, we can obtain the difference equations for node $i$ along the characteristic lines $C^{+}$ and $C$ at time $j$ .

$\displaystyle C^{+}:(P_{{}_{i,j}}^{-}-P_{i-1,j-1})+\frac{c}{A}(Q_{i,j}^{-}-Q_{% i-1,j-1})+\frac{\lambda c^{3}\Delta t}{4DA^{2}}$ $\displaystyle\quad\left(\frac{Q_{i,j}^{-}\left|{Q_{i,j}^{-}}\right|}{P_{{}_{i,% j}}^{-}}+\frac{Q_{i-1,j-1}\left|{Q_{i-1,j-1}}\right|}{P_{i-1,j-1}}\right)=0$ (20) $\displaystyle C^{-}:(P_{{}_{i,j}}^{+}-P_{i+1,j-1})-\frac{c}{A}(Q_{i,j}^{+}-Q_{% i+1,j-1})-\frac{\lambda c^{3}\Delta t}{4DA^{2}}$ $\displaystyle\quad\left(\frac{Q_{i,j}^{+}\left|{Q_{i,j}^{+}}\right|}{P_{{}_{i,% j}}^{+}}+\frac{Q_{i+1,j-1}\left|{Q_{i+1,j-1}}\right|}{P_{i+1,j-1}}\right)=0$ (21)

where, $Q_{i,j}^{-}$ and $Q_{i,j}^{+}$ respectively represent the mass flow before and after the leak point, kg/s; $P_{{}_{i,j}}^{-}$ and $P_{{}_{i,j}}^{-}$ respectively represent the pressure before and after the leak point, Pa.

Equations (4) and (4) represent two equations satisfied at time $j\Delta t$ . However, since these equations involve four unknowns (flow rate and pressure before and after the leak), solving them simultaneously does not directly provide the flow rate and pressure before and after the leak at time $t$ . Therefore, after the occurrence of a leak, it is required to introduce two additional boundary condition equations that relate to the flow rate and pressure before and after the leak point.

In general, assuming that the momentum variation in the horizontal direction caused by the leak can be neglected, the pressure before and after the leak remains constant. Additionally, if the amount of the leak is confirmed, the flow rate relationship after and before the leak can be established based on the mass conservation principle. Under these assumptions, these boundary conditions can be specified as follows:

$\displaystyle P_{i,j}^{+}=P_{i,j}^{-}=P_{i,j}$ (22) $\displaystyle Q_{i,j}^{-}-Q_{i,j}^{+}=K_{i,j}$ (23)

By incorporating the boundary conditions Eqs (22) and (23) at the leak point into the difference Eqs (4) and (4), we obtain a four equations system with four unknowns. This system allows us to calculate the flow rate and pressure after and before the leak point at any given time.

At the endpoints of the pipeline, similar to the normal condition case, we only need to provide boundary conditions regarding the flow rate and pressure. This enables us to determine the flow rate and pressure at various points along the boundary at each time interval. The calculation method for the flow rate and pressure at each node outside of the leak point and pipe endpoints follows the same approach as the calculation method under normal conditions.

By applying these calculations and considering the appropriate boundary conditions, we can determine the flow rate and pressure at various points along the pipe, including leak point, as the transient flow progresses.

5. Calculation method for leak rate

The leak rate is time-dependent, with the maximum leak rate typically occurring at the moment the leakage starts. However, in practical engineering scenarios, it is often sufficient to determine the maximum leak rate, which corresponds to the leak rate under steady-state conditions. Therefore, this part zooms in on the calculation method for the maximum leak rate of a natural gas pipeline.

Figure 3 is a sketch map of natural gas pipeline leakage, assuming a leakage occurs at a distance $X$ from the starting point of the pipe. The starting point of the pipe axis is point 1, the end point of the pipe axis is point 4, the leak point is point 3, and point 2 is the position on the pipeline axis corresponding to the leak point. $P$ , $T$ , $u$ , and $\rho$ represent the pressure, temperature, flow velocity, and density at each point, while $Q^{-}$ and $Q^{+}$ represent the mass flow rates before and after the leak point, and $K$ represents the leak rate at the leak point. When a leak occurs, the parameter values at points 1 and 4, which represent the atmospheric environment, as well as the mass flow rate before the leak point, are known. On the other hand, the parameter values at points 2 and 3 are unknown and need to be calculated using known parameters and the principles of thermodynamics and fluid mechanics. By establishing a leakage model for the pipe, it becomes possible to calculate the leak rate.

Figure 3.

Sketch map of gas pipe leakage.

In general, natural gas pipeline leaks are considered as orifice leaks, and the leak process can be viewed as the outflow of compressible gas through an orifice. Different shapes of leak holes are converted to circular orifices for calculation. For ease of calculation, these assumptions are proposed: (a) the flow in the pipe is one-dimensional; (b) the flow is adiabatic in the pipe; (c) the flow is isentropic at the leak point; (d) the motion of the gas follows the laws of ideal gases. As the assumption is made, that the energy and momentum conservation equations are obtained [11]:

In the context of pipeline leaks, it is common to treat them as orifice leaks, where the leakage process is modeled as the outflow of compressible gas through an orifice. To simplify the calculation, various shapes of leak holes are approximated as circular orifices. So as to streamline the analysis, these assumptions are proposed:

(a)

One-dimensional flow: The flow is assumed to occur in a single dimension, simplifying the calculations in the pipe.

(b)

Adiabatic flow: The flow in the pipeline is considered to be adiabatic, meaning there is no significant heat transfer between surroundings and the gas.

(c)

Isentropic flow at the leak point: It is assumed that the flow is isentropic at the leak point, meaning there is no entropy change during the flow through the leak orifice.

(d)

Ideal gas behavior: The gas is assumed to follow the laws of ideal gases, then the ideal gas equation can describe the relationship of pressure, volume, and temperature.

By assuming adiabatic flow in the pipeline, the conservation equations for energy and momentum can be applied to derive the equations necessary for the analysis [11]:

$\displaystyle\frac{k+1}{k}\ln\left(\frac{P_{1}T_{2}}{P_{2}T_{1}}\right)+\frac{% M}{\textit{RG}^{2}}\left(\frac{P_{2}^{2}}{T_{2}}-\frac{P_{1}^{2}}{T_{1}}\right% )+\left(\frac{4\lambda X_{e}}{D}\right)=0$ (24)

where $G$ is the mass flow in the pipeline (kg/(m².s)), $k$ is the adiabatic coefficient of gas is generally 1.44 for natural gas, $M$ is the molar mass of natural gas (kg/mol) Xe is the length between point 1 and point 2( $m$ ), Pi is the pressure at point $i,i=1,2,3$ (Pa), Ti is the temperature at point $i,i=1,2,3(K)$ , $R$ is the gas constant (J/(mol.K)).

The hydraulic friction coefficient $\lambda$ is determined by pipeline roughness $\delta$ and Reynolds number Re When , then $\lambda={16}\mathord{\left/{\vphantom{{16}{Re}}}\right.\kern-1.2pt}{Re}$ . In the case of transition zone and turbulence, the relationship between $\lambda$ and $R e$ is given by Colbrook’s formula:

$\displaystyle\frac{1}{\sqrt{\lambda}}=-4\lg\left(\frac{1}{3.7}\frac{\delta}{D}% +\frac{1.255}{\textit{Re}\sqrt{\lambda}}\right)$ (25)

When Re is large enough to reach fully turbulent flow, a and b become independent of each other.

For pipes with rough inner walls: $\frac{1}{\lambda}=4\lg\left(3.7\frac{D}{\delta}\right)$

For pipes with smooth inner walls: $\frac{1}{\lambda}=4\lg\left(\frac{\textit{Re}\sqrt{\lambda}}{1.255}\right)$

By substituting Eqs (1) and (4) into Eq. (24), the expression for the leak rate $K$ at the leak point can be obtained under the assumption of adiabatic flow.

$\displaystyle K=C_{D}A_{or}P_{2}\sqrt{\frac{M}{\textit{ZRT}_{2}}\frac{2k}{k-1}% \left[{\left({\frac{P_{a}}{P_{2}}}\right)^{\frac{2}{k}}-\left({\frac{P_{a}}{P_% {2}}}\right)^{\frac{k+1}{k}}}\right]}$ (26)

where $K$ is the leak rate (kg/s), $A_{or}$ is the area of the leak point ( $m^{2}$ ), $C_{D}$ is the orifice flow correction coefficient.

The nature of the gas flow at the leak point, whether it is supersonic or subsonic, has a great influence on the calculation of leak rate. Considering that the atmospheric pressure outside the pipeline is known, this article assumes the critical pressure value (CPR) to be:

$\displaystyle\textit{CPR}=\frac{P_{a}}{P_{2c}}=\left({\frac{2}{k+1}}\right)^{% \frac{k}{k-1}}$ (27)

$P_{2c}$ is the critical pressure, Pa. If $P_{2}$ gradually increases, the leak velocity of natural gas will also increase gradually, until the leak velocity approaches the local sound speed. If $P_{2}$ increase continuously, the leak velocity will no longer increase, but will remain constant, and this constant velocity will equal to the local sound speed. The state is defined as critical state.

When $\frac{P_{a}}{P_{2}}<\textit{CPR}$ , the orifice is a critical flow leakage. Substitute Eq. (27) into Eq. (26) to get the expression of leak rate:

$\displaystyle K=C_{D}A_{or}P_{2}\sqrt{\frac{\textit{Mk}}{\textit{ZRT}_{2}}% \left({\frac{2}{k+1}}\right)^{\frac{k+1}{k-1}}}$ (28)

When $\frac{P_{a}}{P_{2}}\geqslant\textit{CPR}$ , the orifice is a subcritical flow leakage, and the leak rate is expressed by Eq. (26).

When the leakage occurs, the parameter values at points 1 and 4, the atmospheric environment, then the mass flow rate before leak point are known, while the parameter values at points 2 is unknown. After observing Eqs (26) or (28), it can be found that the parameters at point 2 are the key parameters for calculating the leak rate. Therefore, the following research focuses on how to calculate the parameters at point 2.

There are three main differences between natural gas pipeline leaks and classical orifice flow, which are also the key factors to accurately calculate the leak rate These three key factors are:

Frictional pressure drop. Frictional flow of natural gas from point 1 to point 2 causes a pressure drop from $P_{1}$ to $P_{2}$ , which in turn affects the leak rate.

The occurrence of a leakage in the natural gas pipeline will result in a reduction in pressure at the leak point ( $P_{2}$ ). This reduction in pressure is a consequence of the gas escaping through the leak, resulting in a decrease in total pressure in the vicinity of the leak point. The extent of the pressure reduction will depend on factors such as the size of the leak, gas properties, and pipeline conditions.

The presence of a leak affects the upstream flow rate of leak point. The upstream flow rate is equal to downstream flow rate plus leak rate. Therefore, the occurrence of a leak disrupts the normal flow pattern, causing a change in the upstream flow rate of leak point. This change in upstream flow rate, in turn, influences the pressure in the pipe. The altered pressure conditions can impact the leak rate, leading to a further modification in the upstream flow rate. This interplay between pressure, flow rate, and leakage forms a feedback loop, where changes in one parameter affect the others.

Key factor c) is not considered in Montiel’s proposed leak rate calculation method [11], making it an innovative aspect of this study.

5.1 Basic model

The Basic model, which ignores factors a), b), and c), assumes constant parameters at points 1 and 2 before and after the leak. It utilizes Eqs (26) or (28) to calculate the leak rate. While the Basic model is suitable for calculating the flow rate of high-capacity container orifice flow, it is not well-suited for calculating the leak rate. This is because the Basic model overlooks several important factors. Firstly, Once a leak occurs, the total pressure inside the pipeline decreases, particularly at point $P_{2}$ . Secondly, the pressure drop resulting from frictional flow between points 1 and 2 needs to be considered, particularly when the distance between these points is significant. Lastly, a leak alters the upstream flow rate, which further impacts the leak rate. Consequently, employing the Basic model to calculate the leak rate will yield higher values compared to the actual leak rate. This model is only applicable when the leak point is in close proximity to point 1 and leak aperture is small. For example, the calculation of discharge flow rate or leak rate of large fluid storage tanks can be applied to this model.

5.2 Improved model

By neglecting factors b) and c) while considering factor a), the resulting leak rate calculation model is known as the Improved model. This model consider the pressure drop caused by frictional flow from point 1 to point 2, thereby improving the accuracy compared to the Basic model for calculating leak rates. However, as factors b) and c) are disregarded, the leak rate determined by this model will still be higher than the actual leak rate, albeit lower than the leak rate which calculated by this model. Consequently, Improved model is suitable for preliminary leak rate calculations involving small leakage apertures.

Assuming a constant coefficient of friction along the entire pipe, the parameter relationship between point 1 and 2 is as follows:

$\displaystyle\frac{k+1}{2}\ln\left[\frac{\textit{Ma}^{2}_{2}}{\textit{Ma}^{2}_% {a}}\frac{2+(k-1)\textit{Ma}^{2}_{1}}{2+(k-1)\textit{Ma}^{2}_{2}}\right]-\left% (\frac{1}{\textit{Ma}^{2}_{1}}-\frac{1}{\textit{Ma}^{2}_{2}}\right)+\frac{4% \lambda X_{\varepsilon}}{D}=0$ (29) $\displaystyle\frac{T_{2}}{T_{1}}=\frac{Y_{1}}{Y_{2}}$ (30) $\displaystyle\frac{P_{2}}{P_{1}}=\frac{\textit{Ma}_{1}}{\textit{Ma}_{2}}\sqrt{% \frac{Y_{1}}{Y_{2}}}$ (31) $\displaystyle\frac{\rho_{2}}{\rho_{1}}=\frac{\textit{Ma}_{1}}{\textit{Ma}_{2}}% \sqrt{\frac{Y_{2}}{Y_{1}}}$ (32)

where

$\displaystyle Y_{i}=1+\left({\frac{k-1}{2}}\right)\cdot\textit{Ma}_{i}^{2}$ (33)

In the above equations, Ma is the Mach number, that:

$\displaystyle\textit{Ma}_{i}=\frac{u_{i}}{c}=\frac{u_{i}}{\sqrt{\textit{kZRT}% \mathord{\left/{\vphantom{\textit{kZRT}M}}\right.\kern-1.2pt}M}}$ (34)

By using Eqs (29)–(34), the parameters at point 2 can be calculated, and then Eqs (26) or (28) can be employed to calculate the leak rate.

Under these conditions, these relationships can geted:

When $\frac{P_{a}}{P_{2}}<\textit{CPR}$ , critical flow at point 3, then Eq. (28) can be employed to calculate the leak rate.

When $\frac{P_{a}}{P_{2}}\geqslant\textit{CPR}$ , subcritical flow at point 3, then Eq. (26) can be employed to calculate the leak rate.

5.3 Final universal model

If all three factors a), b), c) are considered, then the resulting model is a Final universal model. For the sake of computational convenience, we assume that the pressure $P_{4}$ remains constant in our model. To calculate the leak rate, successive approximation method is employed. The following steps are taken:

Step1:
Using Improved model, we calculate the leak rate $K$ and the pressure of point 2 $P_{2}$ , based on the known upstream flow rate $Q^{-}$ .
Step2:
If $P_{2}$ is greater than $P_{4}$ , we can determine the downstream flow rate $Q^{+}$ , based on pressures $P_{2}$ and $P_{4}$ . By conserving mass, we can calculate the natural gas pipeline flow rate $Q^{-^{\prime}}$ upstream of the leak point, which is equal to $K$ plus $Q^{+^{\prime}}$ . We then insert $Q^{-^{\prime}}$ into step 1 to obtain new values for $K^{\prime}$ and $P_{2}^{\prime}$ , repeating the process until the calculation results for the natural gas leak rate stabilize.
Step3:
If $P_{2}$ is less than or equal to $P_{4}$ , the downstream flow rate $Q^{+}$ is assumed to be zero. By conserving mass, we can calculate upstream flow rate $Q^{-}$ , which is equal to $K$ . We then insert $Q^{-}$ into step 1 to obtain new values for $K^{\prime}$ and $P_{2}^{\prime}$ , repeating the process until the calculation results for the leak rate stabilize.

Figure 4 shows the flowchart for the Final universal model, $\sigma$ is the value of threshold.

Figure 4.
Flow sheet of the Final universal model.

During these outlined steps, it becomes evident that $P_{2}^{\prime}$ and $P_{2}$ undergo changes following each cycle. Therefore, it is crucial to consider the influence of these changes on the flow state within the pipeline. The flow at point 3 can be divided into critical and subcritical, resulting in three distinct flow types at the leak point:

(a). Subcritical flow in the pipe, critical flow at leak point.

$\displaystyle\left\{{\begin{array}[]{l}P_{2}^{\prime}>P_{1}\textit{Ma}_{1}% \sqrt{\frac{2Y_{1}}{k+1}}\\ \displaystyle\frac{P_{a}}{P_{2}^{\prime}}<\textit{CPR}\Leftrightarrow P_{2}^{% \prime}>P_{2c}\\ \end{array}}\right.$ (35)

In this case, Eq. (28) is used in step 1 to calculate the leak rate using Improved model.

(b). Subcritical flow in the pipe, subcritical flow at leak point.

$\displaystyle\left\{{\begin{array}[]{l}P_{2}^{\prime}>P_{1}\textit{Ma}_{1}% \sqrt{\frac{2Y_{1}}{k+1}}\\ \displaystyle\frac{P_{a}}{P_{2}^{\prime}}\geqslant\textit{CPR}\Leftrightarrow P% _{2}^{\prime}\leqslant P_{2c}\\ \end{array}}\right.$ (36)

In this case, Eq. (26) is used in step 1 to calculate the leak rate using Improved model.

(c) Critical flow in the pipeline, critical flow at leak point.

$\displaystyle\left\{{\begin{array}[]{l}P_{2}^{\prime}<P_{1}\textit{Ma}_{1}% \sqrt{\frac{2Y_{1}}{k+1}}\\ \displaystyle\frac{P_{a}}{P_{2}^{\prime}}<\textit{CPR}\Leftrightarrow P_{2}^{% \prime}>P_{2c}\\ \end{array}}\right.$ (37)

In this case, Eq. (28) is used in step 1 to calculate the leak rate using Improved model.

Equation (28) is applied in step 1 to calculate the leak rate using the Improved model. Additionally, as previously mentioned, the pressure $P_{2}^{\prime}$ obtained in each cycle should be chosen based on the aforementioned criteria. By taking into account the influence on the upstream and downstream flow rate, the Final universal model offers a more realistic approach. Furthermore, the iterative approximation method employed eliminates the requirement for intricate mathematical derivations. In particular, when a leak occurs, the size of the leakage aperture is not known, and this Final universal model is suitable for leaks of any aperture size.
5.4 Simulation verification

It is dangerous to conduct an experiment for leak rate measurement on actual natural gas pipelines. Even when using compressed air instead, it is difficult to simulate multiple leaks with different aperture sizes to cover a range of $d/D$ values (the value of the leak point diameter to the pipeline diameter). Therefore, this paper relies on simulations to validate the Final universal model. A simulation of a natural gas pipeline leakage scenario is performed, then the specific parameters are presented in Table 1.

Table 1
Simulation parameters

$Q$ (kg/s)	$k$	$R$ (J/(mol.K))	Re
108	1.334	8.314	3000
$D$ (m)	$E$ (mm)	Pa (pa)	$M$ (kg/mol)
0.216	0.045	10⁵	16.48 $\times$ 10^{- 3}
$P_{1}$ (MPa)	$P_{4}$ (MPa)	$T_{1}$ (k)	Le (m)
18	6.8	293	1300

Figure 5 depicts the correlation between the leak rate and $d/D$ for various models, utilizing the proposed models and parameters. The leak rates derived from the Basic model and the Improved model demonstrate a substantial increase with the rise in $d/D$ . However, only at small $d/D$ values do these two models yield leak rates that align with the predictions of the Final universal model, which falls within their applicable range. It should be noted that the Basic model is suitable solely when the leak point is in close proximity to point 1 and involves a small leak hole.

Figure 5.

Correlation between leak rate $K$ and $d/D$ .

When dealing with very small leak holes, the Improved model is recommended. The diagram illustrates that the leak rate of natural gas calculated using the Basic model is greater than the Improved model, which aligns with the internal friction effect of the latter. As a result, the Improved model is more precise in predicting natural gas leak rates compared to the Basic model. Leakage causes a reduction in pressure at point 2, although not indefinitely, leading to a smaller leak rate calculated by the Final universal model than the Improved model. Nonetheless, the difference is insignificant. Therefore, the leak rates predicted by both models are similar when the aperture is small. However, as the aperture expands, the impact of the pressure reduction on the leak rate at point 2 becomes more pronounced. Consequently, the leak rate will plateau, becoming stable and smooth, eventually matching the flow rate (When the value of $d/D$ is close to 1, that is, the pipe is nearly completely broken, the leak rate is almost equal to the flow rate).

In practical leak scenarios, the flow velocity before leak point tends to increase due to the leak. This increase in velocity leads to a decrease in flow pressure, following Bernoulli’s equation. However, the Improved model does not account for this effect, leading to an overestimation of the leak rate. In cases where the $d/D$ value is small and the leak rate is also small, compared to the undisturbed state, point 2 has the least change in pressure. Consequently, the feedback effect of the leak on the pipeline upstream pressure can be disregarded. With the increasing of the aperture size, the feedback effect becomes more significant, and the Final universal model addresses it by incorporating the iterative model depicted in Fig. 4. The iterative model calculates new pressure and leak rate values for point 2 until the relevant parameters of the leak closely approximate the actual conditions, resulting in a more accurate prediction of the leak rate. This explains why the Final universal model predicts a lower leak rate than the Improved model (as shown in Fig. 5). The aperture size increases with the increase of $d/D$ , leading to a higher leak rate, the feedback effect on gas parameters is more obvious. When the value of $d/D$ is close to 1, the aperture size becomes comparable to the diameter of pipe, causing $P_{2}$ to decrease and approach Pa. Consequently, the leak rate formula converges to Eq. (28), and the leak rates calculated by the Final universal model stabilize to a constant value – the flow rate.

Figure 6.

Correlation between leak rate $K$ and Le.

Figure 6 presents the correlation between the leak rate, $K$ , and the Le for various $d/D$ values in the Final Universal Model (with $P_{1}=$ 18 MPa). As Le increases, the leak rate exhibits a non-linear decrease. This behavior arises from the pressure reduction caused by internal friction within the pipe. Moreover, the severity of the leak varies depending on the $d/D$ value. For small $d/D$ values, the impact of leakage on the natural gas parameters within the pipeline is minimal, resulting in a gradual reduction of the leak rate. As the $d/D$ value falls below a threshold, the effect of leakage on the flow parameters becomes negligible, leading to a negligible decrease in the leak rate. Conversely, with the increase of $d/D$ value, curve becomes more smooth, indicating a smaller reduction in the leak rate per unit increase in Le.

Figure 7.

Correlation between leak rate $K$ and $P_{1}$ .

Figure 7 displays the correlation between $K$ and $P_{1}$ in the Final Universal Model for different $d/D$ values, with Le fixed at 1300 m. The leak rate of natural gas increases linearly with an increase in $P_{1}$ . Furthermore, the severity of the leak varies depending on the $d/D$ value. For smaller $d/D$ values, the leak rate is least sensitive to pressure changes. On the contrary, the higher the $d/D$ value, the higher the sensitivity between pressure and leak rate. This difference is due to the impact of the leak on the gas parameters in the pipe. When the $d/D$ value is small, this effect is small, while it becomes more significant for larger $d/D$ values.

6. Conclusion

This article focuses on the hydraulic parameters of natural gas pipelines, specifically pressure and flow rate, which are the most crucial factors based on the fundamental control equations governing the flow. The article then derives the differential equations corresponding to the basic control equations, using the method of characteristics, to obtain the calculation method of transient pressure and flow rate under normal and leak conditions of natural gas pipelines. Additionally, the article studies the calculation method of the leak rate, which takes into account the pressure drop caused by pipeline friction and uses a successive approximation mathematical method to consider the effect of leakage on leak point pressure and upstream flow rate. Through simulation examples, it has been demonstrated that this method is a universal and highly accurate method that can calculate leak rates for any pipe diameter.

Footnotes

Acknowledgments

The authors acknowledge the Scientific Research Project of Chengdu Vocational & Technical College of Industry (Grant: 20233050002), the Major Joint Research Project of Chengdu Vocational & Technical College of Industry (Grant: XJTD-2023-2) and the Master Studio of Chengdu Vocational & Technical College of Industry (Grant: XJDS-2022-1).

References

Feng

Liang

Nie

Zhang

. A leakage risk assessment method for hazardous liquid pipeline based on markov chain monte carlo. International Journal of Critical Infrastructure Protection. 2019; 27: 201-212.

Jones

Williams

Davis

. Detection and prevention of pipeline leaks using advanced sensing technologies. Journal of Loss Prevention in the Process Industries. 2017; 48: 214-224.

Brown

Johnson

Smith

. Numerical study of the impact of flow rate and viscosity on hydraulic parameters in natural gas pipeline flow. Journal of Natural Gas Science and Engineering. 2015; 24: 106-114.

Wang

Liu

Zhang

. A new method for optimizing pipeline performance based on the method of characteristic line. Journal of Petroleum Science and Engineering. 2018; 165: 651-660.

Chen

Zhao

Wang

Liu

. Gas pipeline leakage detection using the method of characteristic line. Measurement. 2015; 70: 1-6.

Liu

Chen

. Estimation of natural gas pipeline leakage rate based on characteristic line method. Journal of Loss Prevention in the Process Industries. 2016; 44: 162-168.

Wang

Yang

. Gas pipeline leakage detection method based on wavelet transform and characteristic line method. Measurement. 2018; 129: 178-188.

Zhou

Yin

. Gas pipeline leakage rate estimation based on characteristic line method and machine learning. Journal of Loss Prevention in the Process Industries. 2019; 62: 303-319.

Yin

Zhou

. Pipeline leakage detection using characteristic line method and artificial neural networks. Measurement. 2020; 165: 057-071.

10.

Wang

Zhou

Zhang

. Gas pipeline leakage detection based on spectral analysis of characteristic line method. Measurement. 2021; 171: 108-122.

11.

Montiel

Vilchez

Casal

Arnaldos

. Mathematical modelling of accidental gas releases. Journal of Hazardous Materials. 1998; 59: 211-223.

12.

Zhou

. Thermal analysis of slow discharge form a pressurized natural gas storage tank. Applied Thermal Engineering. 1997; 17: 1099-1110.

13.

Young

Bum

. A simple model for the release rate of hazardous gas from a hole on high-pressure pilelines. Journal of Hazardous Materials. 2003; 97: 31-46.

14.

Chebouba

Yalaoui

Smati

Amodeo

Tairi

. Optimization of natural gas pipeline transportation using ant colony optimization. Computers & Operations Research. 2009; 36: 1916-1923.

15.

Kostowski

Skorek

. Real gas flow simulation in damaged distribution pipelines. Energy. 2012; 45: 481-488.

16.

Bariha

Mishra

Srivastava

. Hazard analysis of failure of natural gas and petroleum gas pipelines. Journal of Loss Prevention in the Process Industries. 2016; 40: 217226.

17.

Demissie

Zhu

Belachew

. A multi-objective optimization model for gas pipeline opevaluens. Computers & Chemical Engineering. 2017; 100: 94103.

18.

Hou

Jiao

Ren

Cao

Song

. Experimental study of leakage detection of natural gas pipeline using FBG based strain sensor and least square support vector machine. Journal of Loss Prevention in the Process Industries. 2014; 32: 144-151.

19.

Zeng

Dong

Xie

. A new analytical method of strain and deformation of pipeline under fault movement. International Journal of Pressure Vessels and Piping. 2019; 172: 199-211.

20.

Zhu

SCM

Patil

Wang

Hirsch

Song

. Underwater pipeline impact localization using piezoceramic transducers. Smart Materials and Structures. 2017; 26: 107-116.

21.

Jiang

Ren

Jia

. Application of FBG Based Sensor in Pipeline Safety Monitoring. Applied Sciences. 2017; 7: 540-551.

Research on the hydraulic parameters and leak rate of natural gas flow in natural gas pipelines based on the methods of characteristic line and successive approximation

Abstract

Keywords

1. Introduction

2. Fundamental control equations for natural gas flow

5.2 Improved model

Table 1 Simulation parameters

Footnotes

Acknowledgments

References

Table 1
Simulation parameters