Detecting Optimal Leak Locations Using Homotopy Analysis Method for Isothermal Hydrogen-Natural Gas Mixture in an Inclined Pipeline

: The aim of this article is to use the Homotopy Analysis Method (HAM) to pinpoint the optimal location of leakage in an inclined pipeline containing hydrogen-natural gas mixture by obtaining quick and accurate analytical solutions for nonlinear transportation equations. The homotopy analysis method utilizes a simple and powerful technique to adjust and control the convergence region of the infinite series solution using auxiliary parameters. The auxiliary parameters provide a convenient way of controlling the convergent region of series solutions. Numerical solutions obtained by HAM indicate that the approach is highly accurate, computationally very attractive and easy to implement. The solutions obtained with HAM have been shown to be in good agreement with those obtained using the method of characteristics (MOC) and the reduced order modelling (ROM) technique.


Introduction
One of the strategies to reduce gas transportation costs is the use of natural gas pipeline networks by petroleum companies [1]. These networks are capable of supplying gas in long distances under high pressure and through compression stations [2]. Changes in pipeline pressure are a function of gas velocity, valve closure time, and arrangement of the closing valve [3].
When the valve is closed at the end of the pipeline, there is the possibility of the occurrence of maximum pressure, which can be decreased in short times during its closure. It is of utmost importance to control factors affecting transient pressure, such as initial pressure and mass ratio. This is because the damage caused by this pressure is not evident shortly after the event [4][5][6].
Several studies have been conducted on transient flow in the mixtures of hydrogen and natural gas with the use of isothermal flow and horizontal pipelines, which is not the case in reality [2,[7][8][9]. Furthermore, another study has made an attempt to study the flow of these mixtures under high pressure through inclined pipelines [10]. In most pipelines working under high pressure, there are slow and fast fluid transients. As gas properties are not constant, a one-dimensional and non-isothermal gas flow model should be presented to simulate these transients [2].
The reason for proposing hydrogen and natural gas mixtures is their transportation through the same pipelines for the purpose of cost reduction. This is while the existing lines are just designed for natural gas, whose properties are significantly different from that of hydrogen [11,12]. The solution to this problem has been the mixture of the both with a great deal of care and attention, as hydrogen is a reactive gas with high pressure that can cause leakage [13,14]. This problem is of great importance Figure 1 shows an inclined pipeline, which has a reservoir at the top and a valve at its bottom. The governing equations consist of three partial differential equations that are all coupled, non-linear and hyperbolic. The non-isothermal flow in the pipeline, a homogenous mixture of hydrogen and natural gas, was considered to be one-dimensional that is compressible and covers transient condition [7].

Governing Equation
The governing equations for the transport of hydrogen/natural gas mixture in an inclined pipeline from the principle of conserving mass and momentum are given by the following, with u = Q/A, A = πD 2 /4, ρ is density, u is the gas velocity, P is the pressure, e is the gas internal energy per unit mass, D is the diameter of the pipeline, f is the coefficient of friction, g is the gravitational force and θ is an angle between the friction force and the direction. Boundary conditions of this equations depend on the types of closure and the valve operational time. The boundary conditions at the initial point x = 0 and at the end point x = L, respectively are given by, where ρ 0 and u 0 are defined as density and gas velocity at the inlet pipeline, respectively and ρ L and u L are defined as density and gas velocity at the outlet pipeline, respectively. The initial conditions that are assumed to be in a steady state condition at t = 0 are [7], The commonly used equation of state for perfect gas is as follows: where, R: is the specific gas constant. T: is temperature.
The equation of state for the compressible flow, where there is a celerity pressure wave, is: The following equations are also achieved from ideal gas relation, where, C v : is the specific heat at constant volume. C p : is the specific heat at constant pressure. R: is the specific gas constant. P: is pressure.
γ: is the flow process index.

Hydrogen-Natural Gas Mixture Equation
The mass ratio and the density of hydrogen-natural gas mixture are defined as, Where m g , m h , V g and V h are defined as the mass of natural gas and hydrogen and volume of natural gas and hydrogen, respectively. Therefore, the expression of the average density of the gas mixture is given by, The celerity pressure wave for compressible flow is defined as, where the subscript s is defined the constant entropy condition. The derivative of Equation (11) with respect to P, and substituting into Equation (12), then the celerity pressure wave yields [7], The properties of hydrogen and natural gas used in the calculations are shown in the Table 1. For the simulation, the parameters are assumed as Table 2.

Homotopy Analysis Method
A brief description of the standard homotopy analysis method (HAM) presented by [28][29][30][31][32]. This will be followed by a description of the algorithm of the homotopy analysis method (HAM). First, we consider the following differential equation, where N are nonlinear operators, x and t denotes the independent variable, u(x, t) are unknown functions, and G (x, t) are known analytic functions. For G (x, t) = 0, Equation (14) reduces to the homogeneous equation. By means of generalizing the traditional homotopy method, Liao [28] constructed the so-called zero-order deformation equation, where p ∈ [0, 1] is an embedding parameter,h are nonzero auxiliary functions, L is an auxiliary linear operator, u 0 (x, t) are initial guesses of u(x, t), H (x, t) denotes a nonzero auxiliary function and Ψ(x, t; q) are unknown functions. It is important to note that one has great freedom to choose auxiliary objects such ash and L in HAM. Obviously, when q = 0 and q = 1, Equation (15) becomes, Thus, as q increases from 0 to 1, the solution Ψ(x, t; q) varies from the initial guesses u 0 (x, t) to the solutions u(x, t). Expanding Ψ(x, t; q) in Taylor series with respect to q, one has where, If the auxiliary linear operator, the initial guesses, the auxiliary parametersh, and the auxiliary functions are so properly chosen, then series Equation (17) converges at q = 1, and one has, which must be one of the solutions of the original nonlinear equations, as proved by Liao [28]. Ash = −1 and H (x, t) = 1, Equation (15) becomes, which is used mostly in the homotopy perturbation method. Define the vectors, Differentiate the zeroth-order deformation Equation (14) m-times with respect to q and then dividing them by m! and finally setting q = 0, we get the following mth-order deformation equation, where, with, It should be noted that the linear Equation (22), which has linear boundary conditions, governs u m (x, t) for m ≤ 1 [33]. Boundary conditions stem from the main problem, the solution for which can be provided by Matlab, Maple, or Mathematica. The requirement for the limit of Equation (17) is that it should meet the conditions of the main equation N u(x, t) = 0 when it is convergent at q = 1. It is noteworthy that drawing "h-curves" or "curves for convergence-control parameter" aim to find a proper convergence-control parameterh, a convergent series solution, or and accelerate its convergence rate. It is such that these curves with unknown quantities are drawn againsth to approximately find the convergence region, though they are just graphical. This is because it is not possible to find which h 0 ∈ R h provides the fastest convergent series (see Liao [28,34] for further reading). Another note to be made is that a unique solution is achieved when Equation (14) accepts a unique solution; otherwise, many possible solutions will be obtained from HAM.

Solving the Steady State Equations by High-Order Deformation HAM
We define the vectors, Differentiating Equations (5) and (6) m times with respect to the embedding parameter q and then setting q = 0 and finally dividing them by m!, we have the so-called mth-order deformation equations, with the initial conditions, where, with the celerity pressure wave c i defined as follows, with the following linear operators, with the property that, which implies that, Now, the solution of the mth-order deformation Equations (5) and (6) becomes, which can be easily solved by a symbolic computation software such as Matlab, Maple, and Mathematica. Therefore, we will have P(x) and u(x) as follows, Furthermore, to construct the zeroth-order deformation equations we can define the nonlinear with the celerity pressure wave c defined as follows,

Solving Isothermal Flow of Hydrogen-Natural Gas Mixture by HAM
We define the vectors, Differentiating Equations (1) and (2) m times with respect to the embedding parameter q and then setting q = 0 and finally dividing them by m!, we have the so-called mth-order deformation equations, with the initial and boundary conditions as follows, where, with the celerity pressure wave c i defined as follows, with the following linear operators, with the property that, which implies that, Now, the solution of the mth-order deformation Equations (1) and (2) becomes, which can be easily solved by a symbolic computation software such as Matlab, Maple, and Mathematica. Therefore, we will have P(x, t) and u(x, t) as follows, Furthermore, to construct the zeroth-order deformation equations we can define the nonlinear operators N 1 Ψ 1 (x, t; q) and N 2 Ψ 2 (x, t; q) as follows, with the celerity pressure wave c defined as follows,

Results and Discussion
For solving the Equations (5) and (6) by suing the homotopy analysis method according the Equations (25)-(36) we can have, Equation (47) is a approximation solution for pressure P to the problem Equations (25)- (36) in terms of the convergence parametersh and order m = 12 with H (x) = 1. To find the valid region ofh, theh-curves given by the 12th-order HAM approximation at different values of x are drawn in Figure 2; this figure shows the interval ofh in which the value of P 12 is constant at certain x, and M; we chose the horizontal line parallel to x-axis (h) as a valid region which provides us with a simple way to adjust and control the convergence region.    (1) and (2) with the homotopy analysis method (Equations (37)-(46)) using the following initial approximations, we guessed the initial approximations Equations (48) and (49) using the initial and boundary conditions (for x = 0, L results will be as Equations (3) and (4)). Therefore, using the homotopy analysis method for solving the Equations (1) and (2) with the initial approximation Equations (48) and (49) can obtain the following results, therefore, pressure P(x, t) is as follows, Equation (50) is a approximation solution for pressure P(x, t) to the problem Equations (1) and (2) in terms of the convergence parametersh with H (x) = 1. To find the valid region ofh, theh-curves given by the 12th-order HAM approximation at different values of t and x = 0 are drawn in Figure 4; this figure shows the interval ofh in which the value of P 12 (0, t) is constant at certain t, and M; we chose the horizontal line parallel to t-axis (h) as a valid region which provides us with a simple way to adjust and control the convergence region.  (2) at different values of t and x = 0.

Leak Detection Using Homotopy Analysis Method
Because of a small orifice between the high-pressure pipeline and the environment, the orifice of leak can be simulated leaning on the flow rate. The discharged flow from the orifice can be computed by the following Equation [7], where A l is the leak orifice area with radius r l , P l is the pressure of gas mixture at the leak position and ρ l is the density of gas mixture at the leak position respectively, C d is a discharge coefficient and X L is the distance of leak from the reservoir. Analyzing transient pressure wave for hydrogen/natural gas mixtures is based on transmission and reflection properties of pressure wave effected by a downstream valves sudden closure. When the initial pressure wave reaches the leak, it will produce a reflection as it arrives back at the downstream end section. Then, the difference in time between the initial transient wave and the reflected wave is measured and the leakage position in the pipeline is computed by, where X L is defined as the distance between the leak and upstream end section, ∆t l is the difference of time between the initial transient wave and reflected wave and c ∆t l is defined as the transient celerity wave at time. The transient pressure of mixture of natural gas and hydrogen with a mass ratio of φ = 0.5 is shown in Figure 7 in case of isothermal flow and leak location at X L = L/3 with diverse angles. The homotopy analysis method from order 12 andh = −0.5 has been used. Red line is for θ = 0 and black line is for θ = 15.

Results and Discussion
The celerity wave distribution is presented in Figure 8 as a function of time. In this case, the valve of the horizontal pipeline containing different mass ratios of a mixture of gas and hydrogen is abruptly closed when the leakage is at X L = L/3. The values of celerity wave of the leak point for various mass rations are 819.20 ms −1 , 964.60 ms −1 and 1086.60 ms −1 for φ = 0.25, 0.5 and 0.75, respectively.
As shown in Figures 6 and 7, the occurrence of the leakage is possible when ∆t l is equal to 0.808 s. Equation (53) can be used to calculate the leak location of the mixture of natural gas and hydrogen in case of an isothermal flow in a horizontal pipeline as follows: As seen earlier, there are various mass ratios of the mixture and various angles of the pipeline each with a specific leak location at X L = L/3, the values of which are presented in Table 3. It can be inferred that the leak location is not a function of pipe angle, it is rather a function of the mass ratio of the natural gas and hydrogen mixture. Therefore, mass ratio is of utmost importance here.
The real location of leak is 200 m, when the leak location is at X L = L/3. The leak location calculations by Subani et al. and HAM turned out to be 211.10 m and 210.30 m, respectively. It is a mixture of natural gas and hydrogen with a mass ratio of 0.5. When the mass ration is increased to 0.75, the leakage location is less than 200 m. Therefore, when the mass ratio is decreased, the location is greater than 200 m. This is contrary to the calculations since the calculated value is less than 200 m when the mass ratio is 0.5. This is an indication of the dependence of leak location of mass ratio of the mixture considered. As Elaoud et al., (2010) state, the most important part in early determination of a leak close to the reservoir or compressor is the bottom of the pipeline. Figure 9 shows the leak location with respect to the gas mixture (φ). As can be seen from this figure, there is a steep slope for the values φ ∈ [0, 0.25] and φ ∈ [0.75, 1], but for values φ ∈ [0.25, 0.75] there is a mild slope. Table 3. Leak location for the hydrogen-natural gas mixture for isothermal flow at leakage X L = L/3.

Gas Mixture
Pipeline      In real (physical) pipelines, noise is expected to affect measurements [35,36]. The possible effects of noisy signals on the performance of the proposed method are Brownian motion or Wiener process or White noise, as the physical model of the stochastic procedure, as an indexed collection random variables. A Wiener process (notation W = (W t ) t≥0 ) is named in the honor of Prof. Norbert Wiener; other name is the Brownian motion (notation B = (B t ) t≥0 ). Wiener process is Gaussian process. As any Gaussian process, Wiener process is completely described by its expectation and correlation functions. A Brownian motion, also called a Wiener process, is obtained as the integral of a white noise signal as follows, The effects of noisy signals on the effectiveness of the proposed method and possible effects of noisy signals on the performance of leak locations will be proposed in the future works, by introducing white noise in the simulations.
For accurate pinpointing, we can use the zero-gradient control (ZGC) method which we have discussed in our recently published paper [6] about optimal mixture and controlling the pressure. In our next manuscript with title "Detecting Optimal Leak Locations using Delta Method and Zero Gradient Control for Non-isothermal Hydrogen/Natural Gas Mixture in an Inclined Pipeline" we used the delta method (DM) and zero gradient control (ZGC) method for detecting optimal leak locations. In our future works we will mixed the proposed methods with Artificial intelligence, Neural Network and Deep Learning [37] to predict and estimate the optimal mixture parameter for achieving more accurate pinpointing.

Conclusions
The homotopy analysis method used to solve the flow equations of hydrogen natural gas mixture in an inclined pipeline. To validate the approximation series for pressure compared with the Subani et al. The leak locations were detected using the homotopy analysis method for horizontal pipeline (θ = 0 • ) and inclined pipeline (θ = 15 • ) for gas mixture φ = 0, 0.25, 0.5, 0.27, 1. Using the homotopy analysis method the celerity wave at leak point of the pipeline are 819.20 ms −1 , 964.60 ms −1 and 1086.60 ms −1 for φ = 0.25, 0.5 and 0.75, respectively.
In an inclined pipeline θ = 15 • the leak location for gas mixture φ = 0.5 using the Subani et al. method (ROM) and homotopy analysis method respectively are 211.1 m and 210.3 m. Because of the real leak location is supposed at 200 m when the leak is located at X L = L/3, the result of HAM method is more accurate than ROM method. As can be seen from Figure 9, with increases the gas mixture φ from 0 to 1 the leak location decreases and there is a steep slope for φ ∈ [0, 0.25] ∪ [0.75, 1], and a mild slope for φ ∈ [0.25, 0.75].
The proposed HAM method is employed without using linearization, discretization, or transformation. It may be concluded that the HAM is very powerful and efficient in finding the analytical solutions for a wide class of gas transportation equations in a pipeline.