Pipeline Leak Detection by Transient-Based Method Using MATLAB R © Functions

This paper shows a method for pipeline leak detection using a transient-based method with MATLAB R © functions. The simulation of a pipeline systems in the time domain are very complex. In the case of the dissipative model, transfer functions are hyperbolic Bessel functions. Simulating a pipeline system in the frequency domain using a dissipative model we could find an approximate transfer function with equal frequency domain response to in order get the pipeline system’s time domain response. The method described in this paper can be used to detect, by comparison, to detect a leak in a pipeline system model.


Introduction
Pipeline networks are widely used worldwide for the transportation and distribution of natural gas, water and other "light" petroleum products.Natural gas and petroleum are transported long distances from Oil platforms and refineries to customers and consumers.These pipelines cross villages, industrial parks, deserts, forests, etc. ; wherein the detection of leakage is not easy but their presence is potentially dangerous.One of the biggest difficulties that affect the operational safety of these lines is leakage.Leaks may be caused by corrosion, pressure surges, earthquakes, etc. Therefore,leak detection and location are of utmost importance.Several authors have developed leak detection and / or leak location systems.Leak detection methods are based on "hardware" methods or, where appropriate, on "software" methods (Clasification based on Murvay, P.S. [1]).Between the systems using the "hardware" method we found the "pigging" (Furness, RA [2]), acoustic methods (Kim, M [3]), gas tracing methods (Lowry, W.E. [4]), cable sensor method (Sandberg, C. [5]), optical fiber method (Tapanes, E. [6]) and methods based on infrared thermography (Weil, G.J. [7]).These systems provide a precise location of the leak, however its implementation is costly and complex.Moreover, these systems are used to assess the continuity of the system and are not real time methods.Software based methods, as the name states, have software programs at their core.The implemented algorithms continuously monitor the state of pressure, temperature, flow rate or other pipeline parameters and can infer, based on the evolution of these quantities, if a leak has occurred.The software methods can use different approaches to detect leaks: mass/volume balance (Liu, J. [8]), real time transient modeling (Hauge, E. [9]), acoustic/negative pressure wave (Mpesha, W. [10]), pressure point analysis (Scott,S [11]), statistics (Zhang, J. [12]) or digital signal processing (US Department of Transportation [13]).From a technical point of view, pneumatic and hydraulic transmission lines often belongs to a dynamic system that should be analyzed in the time domain.The complexity of the equations governing these physical phenomena makes it difficult to find an analytical solution to the problem.This article proposes a methodology to find a model formed by a transfer function obtained from frequential analysis of the equations.The time domain simulation of a transient fluid system is quite complex.In the case of a dissipative model the transfer function consists of hyperbolic Bessel functions.The whole system simulation is combines all the system equations, including the transitional functions, in a single transfer function.Simulating equations transformed into frequency domain, we could obtain a transfer function that shows equal response both in the time and the frequency domain.A simulation of a simple break system is used to illustrate the method.For that, a theoretical model is developed in order to show the behavior of the brake system simulated with and without a leak.

Modeling a single Pipeline
A single pipeline (Fig. 1) could be modeled as nonlinear system with distributed parameters.The analytical solution for this unsteady flow problem is reached from continuity, momentum and energy equations in a cylindrical coordinate system.These equations are: Momentum: Continuity: Energy: A compressibility equation (for liquids) or an equation of state (for gases/steams) are also required: Details on the model derivation can be found in GOODSON, R.E.[14] , NURSILO, W. S. [15] and KING, J. D. [16].Laplace transform of the partial derivative equations lead to: where: Any model derived from these expressions uses a propagation operator Γ(s) and the characteristic impedance Z c .The Propagation operator governs the pressure between two points of the line through the following relation: and the characteristic impedance governs the flow due to the following expression: For the purpose of frequency response comparison, we use normalized parameters respect to time are used: Normalized Laplace Operator s where ω ν represents the Viscous Frequency Replacing the Laplace operator with the normalized operator s = sω ν the normalized Propagation Γ( s) remains as follows: where D n is the Dissipative Number or Dissipative constant and ω c is the characteristic frequency For Bessel parameter B r , J 0 is the Zero order Bessel Function and J 1 is the First order Bessel Function Parameter Bessel B rσ (where σ is the Prandtl Number): The Dimensionless Dissipative Number D n is used as a reference parameter to compare different frequency responses.To evaluate the effect of this dissipative constant an m-file code was prepared.This code shows the Magnitude and phase plots of outlet pressure respect to normalized frequency as shown in Laplace transform of the partial derivative expressions 7 and 8 could be rewritten in several matrix forms depending on the inputs and outputs selected.Using x = 0 as pipeline inlet and x = L as the outlet, the matrix form for a pressure input and a flow output is as follows: Equation 19 represent 2 equations and four variables.Knowing flows, solving Matko [17] introduced a simplification for the one-dimensional case yields.Introducing the mass flow rate ṁ, the Continuity equation could be written as: A (Pipe Section) c (Sound Speed) The Momentum equation (neglecting the effect of gravity) is as follows: where λ( ṁ) is the pipe friction coefficient as a function of the mass flow rate.Both equations 21 and 22 are linearized and the corresponding PDE system is: where L = 1/A is the system inductance, R = λ( ṁ) ṁ/A 2 ρd is the Fluid Resistance and C = A/c 2 is the System Capacitance.The analytical solution of a pipeline system model is obtained by derivation of this linear equations.Using this linearized solution, the equivalent matrix of eq.19 is: where n = (Ls + R) Cs, l is the length of pipe, and Z k = (Ls + R) / Cs .This simplification introduced by Matko has been tested and the results are similar to the results obtained using the equations 19.Using these linearized equations do not result in computational resources savings, thus Equation 19 was used to solve the model.
For a fluid with ρ = 870kg/m 3 , µ = 46 * 10 −6 m 2 /s, β e = 1.07 * 10 7 Pa, L = 2m and inner pipe diameter of d = 0.01m, the solution of this dissipative model leads to the frequency response magnitude and phase plots shown in Figure 3.
Using invfreqs MATLAB R function we obtain the folowing transfer function:

Modeling a single Pipeline with leak
A pipeline with a leak can be modeled as two leak free pipelines connected in series as shown in Figure 4.The leak is modeled in the two pipelines joint as an atmospheric flow sink, considering its flux as directly proportional to the pressure gradient and inversely proportional to leak resistance.R l is a linear approximation for leak flow resistance, ∂p/∂Q.Higher resistances correspond to smaller leaks.Also, mass conservation equation must be fulfilled at the leak joint:  From equation 19 for each line on this scheme we obtain: (Note:The variable X ij means the value of property X in the pipeline end i of the line j) : Line a: Line b: Equations at leak: Using the same fluid and pipeline parameters of the previous section, including a leak at the center of the line L a = L b = 1 2 • L, the obtained frequency response of the model is shown in Figure 5.It can be concluded that one effect of the leak is to modify the second order modes of the response, reducing their amplitude and their frequency.Step and impulse response of the system with and without leak Obtaining the corresponding rational polynomial transfer functions representing a linear ordinary differential equations of a simple line with and without leak, we could excite the models by a step or impulse signal.Obtained results are shown in Figure 6.

Oil Brake System
In this section a theoretical, approximated model for a n oil brake system is developed based on previously published results by Wongputorn et al. [18].In this article, a time domain simulation for hydraulic transients of a brake system is presented for an oil at two different temperatures.Results presented show that in certain circumstances the brake response is slow due to the temperature and suggests an appropriate pipe diameter to overcome this problem.We are going to use the same brake system to simulate the effect of a leak.
The scheme of the break system without leak is shown in the Figure 7 The oil brake system is composed by a pipe system, a brake mechanism (that actuate when receives a pressure increment) and a pedal pump.The pipe system is the group of pipes that connect the pedal to the four brake mechanisms (one for each wheel).The break mechanism is modeled as a tank able to accumulate oil.The pedal is represented as a small tank with a mechanism to get small oil flow Q 11 and increase the pressure (as a pressure wave) P s when braking.The equations that characterize the system, following the nomenclature of the Figure 7 are (Note : The variable X ij means the value of property X in the pipeline end i of the line j): Line 3: Line 4: Mass conservation in nodes s1, s2, s3: The equation system could be also written in the matrix form: A solution for the pressure of the brake mechanism P 24 could be found solving this equation system with the input P s .The solution was reached using Matlab's R solve function as is shown in the m-file appendix A. This solution could be represented as its frequency response magnitude and phase plots.On figure 8 the frequency response of the analytical model (blue) and the approximated model according the Transfer Function (58) (in red) are represented.Also, in time domain, we could represent the excitation of the approximated transfer function with a step and impulse signals (Fig. 9).
The transfer function P 24 (s)/P S (s) = H( s) with the approximated frequency response as the analytical model is:

Oil Brake System with leak
In this section a model for the previous brake system including a leak is presented.The leak is placed in pipeline 2 as is shown in Figure 10.Pipeline is divided at leak point in 2 sub-pipelines, pipeline 2a and pipeline 2b.Two additional equations must be included, mass conservation equation at leak and leak flow equation: Line 2 is split in Line 2a and line 2b: Mass conservation at leak and leak flow equations are as follows: The equation system for the break circuit with a leak in matrix form is not represented due to its size.MATLAB c code was coded to calculate the frequency response for the analytical model, in order to find the approximate transfer function and to plot the results for four different leak rates.In the figure 11    Looking at Figures 11 and 12 it can be observed that the leak affects at the frequency response of the system but do not affect the time response velocity of the brake itself.Keep in mind that we are analyzing transients, and then, we are including the effect of loosing oil after some time.Then we will realize that there is a leak in a brake system when we lose some amount of fluid.If a large leak is simulated, the results show that the brake pressure do not reach the expected value, but time response is similar.
Magnitude and position of the transfer function poles obtained for the simulated cases are summarized in Table 2. Larger leaks result in lower pole magnitudes located at higher frequencies.By looking at the transfer function obtained using the proposed method we were able to determine the presence of a leak on the pipe system.

Conclusions
The governing equations for a fluid transmission system are very complex and do not present an analytical solution that allows to represent its behavior in the time domain.Based on the frequency response of the system, we could obtain an approximated model with equal frequency response and obtain its transfer function.The approximated transfer function model allows to analyze the time domain response, and if needed, compared it with the response of a real model.This comparison, between the approximated model and the real model, could allow us to detect, via software, the presence of leaks.

Figure 2 .Figure 2 .
Figure 2. Effect of the Dissipative number into pressure wave magnitude

)Figure 3 .
Figure 3. Frequency response magnitude and phase plot for the system shown on 1.

PreprintsFigure 6 .
Figure 6.Step and impulse response of the system with and without leak

Table 1 .
Poles and Zeros of the brake system TF without leak

Table 2 .
Absolute values of Poles of TF per Leak value

Table 3 .
Magnitude of Poles of TF per Leak value