A Novel Approach to Model a Gas Network

.


Introduction
The continuous uninterrupted supply of Natural Gas (NG) is crucial to today's economy. The balancing of gas supply from a wide range of sources with various end users can be challenging due to the unique and different behaviours of the end users, which in some cases span across a continent [1,2]. NG is one of the strategic primary sources of energy in the world, supplying approximately 24% of the world's primary energy in 2016 [3]. This increases to 30% in Ireland, where on occasion, 80% of peak power demand is provided by NG [4], and in Singapore, where generated electricity can be up to 95% [5] from NG. Transport of NG via pipeline is the most efficient and economic and has the lowest carbon footprint [6][7][8], and in many instances, the gas network is capable of supplying NG from gas fields directly to the end users (from Well to Wheels; WtW is a common expression for the life cycle of a product). A gas network includes both transmission and distribution pipelines, compressor stations to bring the transmission network to its operating pressure and City non-linear system describing the model, and by using matrix algebra, an alternative linear system is created. The solution of the alternative linear system can be used as a smooth approximation of the solutions of the non-linear system. This formulation is ready for software implementation and may also be used in atomic scale models as an alternative to existing empirical approaches with pair and cohesive potentials. An illustrative example, analysing a local region of a node, is presented to demonstrate the model performance.
In this study, two simplified gas networks were created, each based on the Irish Gas Network. Network 1 is A gas network with 13 nodes and 14 pipes, and a bigger network, Network 2, has 109 nodes and 112 pipes. Both networks are used for comparison of the novel approach with numerical solutions. Having a new reliable method for gas network modelling with fast calculation and precise results can be an option to replace numerical solutions, especially for large gas networks, such as those in continental Europe or North America, which contain several hundred nodes and pipelines and where large computation power is required. This new approach is a computational method that helps in producing faster results, and when the modelling of NG networks is coupled with other energy systems, for example the electrical network, faster and less computationally-demanding results are of key importance.

Literature Review
Various methods have been published to model a gas network and solve the governing equations of gas flow through the pipelines, with the vast majority applying numerical methods due to the non-linearity of the equations. In the literature, a Steady-State Model (SSM) is frequently used where the pressure (the driving force) and flow rate at the inlet and outlet of the pipeline are assumed constant. The steady-state equations of gas flow are non-linear algebraic, and in most cases, for solving this type of equation, numerical solutions are employed. Some valuable studies of SSM have been published. The principles and fundamentals of graph theory can be used to model a network along with the empirical flow equations for low-, medium-, and high-pressure flow. Various different empirical flow equations have been proposed such as Lacey's equation, the Polyflo equation, the Panhandle equation, and the Weymouth equation (for high-pressure networks). Osiadacz [19] applied the Newton-Raphson method to deal with the non-linear equations. In this work, a step-by-step, Newton-Raphson method using a Jacobian matrix in graph theory was presented to model a horizontal gas pipe network. The essence of the work is the presentation of a computational method based on physical and mathematical laws for simplification of fundamental equations. While this work is one of the leading references in the modelling of gas networks, no specified example for a pipeline with inclination was presented. Many authors for simplification neglect the inclination term; however, Matko, Kwabena, and Iterran-Gonzales [22][23][24][25] are some authors who have considered inclination in the pipe flow equation, because of its importance for the validation of models with real-world data. Simulation of a natural gas network focusing on gas quality using a numerical-iteration method was carried out by Abeysekera [26]. An important point of this literature is the inclusion of gases with different density in the network. Varying gas composition will impact gas pipeline system characteristics such as pressure and flow rate. With a move to decarbonise the NG network, this is vital in the evaluation of the networks' ability to store synthetic natural gas or hydrogen. Szoplik [18], after describing three different gas networks, low, middle, and high pressure, indicated two methods to run the gas system model, which are the nodal method and loop method. In the nodal method, the flow rate is the dependent variable, and the pressure drop (differential pressures between two nodes of a pipe/branch) through each branch is the independent variable Q ij = f (dP ij ), while the loop method describes that pressure is a function of flow rate (Equation (1)). The article investigates the impact of ambient temperature on the flow rate.
Debebe's [27] paper focused on simulation of transmission pipelines in a gas network focusing on compressor stations. The equations are mainly based on the mass balance, flow principles, and compression characteristics of gas in pipes. The Newton-Raphson method has been applied to solve the equation using C++. The model evaluated the energy consumption for various configurations in order to optimize the transmission pipeline flow and pressure. There are a number of papers that have implemented the MATLAB-Simulink tool to simulate pipelines: Behbahani and Bagheri [28] and Herran-Gonzalez in [24] applied the Simulink toolbox in MATLAB to model the gas network. Behbahani and Bagheri [28] have manipulated the Simulink model of a gas system, and in order to verify the accuracy, the results have been compared with outputs from the finite difference scheme (FSM) and experimental test cases. The proposed model can predict the transient response of the dependent variable (pressure) with a good correlation with the nonlinear finite difference models. Herran [24] has implemented two models, the characteristic method and Crank-Nicholson, both using the MATLAB-Simulink library for solving the partial differential equations. The results for a simple triangle gas network (including two demand nodes and one supplier node) for both models produced similar results to the Osiadacz model [19]. Mohring [29] applied a numerical method to run the simplified loop model of a gas network. In this study, the inclination term has been neglected, and during simplification, the compressibility term was held constant. Matko [25] discussed the non-linear equations and presented a solution by linearizing the unsteady-state equations where density, viscosity, and temperature (in adiabatic flow) are constant, but in reality, the non-linear flow equation consists of a number of interdependent variables, e.g., flow through a pipeline gas flow rate relates the squares of the inlet and outlet pressures [30].
Yue-Wu [31] believed there were two important aspects in the modelling of gas systems; numerical solutions and optimization, where a numerical solution is applied to determine the actual behaviour of gas flow in the pipe. However, it is highlighted that the disadvantage of the numerical method is that it is not able to put constraints on the system to find minimal cost, which is necessary for optimization. Furthermore, due to the fact that sizing of pipe diameters depends on the users' demand, optimization could lead to significant operational advantages. The aim of optimization is to find the best set of control variables within a single optimization run. For some recent optimization methods useful for this we refer to [32][33][34].
Along with the equation parameters, the solving of non-linear equations is dependent on both the initial conditions and tolerance level chosen by the programmers, which can increase the computational time required and impact the convergence of the model (accuracy of the solution). Alternative methods of solving the pipe flow equations such as the one presented in this work are therefore of interest, as when a large-scale network is modelled, the numerical iteration method will take more time to solve the equations.
In recent years, there has been a significant development in using matrix theory to study networks related to engineering problems [35][36][37][38] and the solutions of non-linear algebraic systems [39][40][41][42][43][44][45]. The idea is to provide new techniques and methods ready for software implementation in order to solve non-linear equations, similar to those for modelling gas pipelines. This paper presents a novel method for solving the equations governing flow and pressure in a natural gas network using matrix theory to derive a unique theorem.

Modelling of a Gas Network
A large amount of NG is transported long distances through the high-pressure pipeline network. Therefore, modelling/simulating of a gas transmission system with sufficient accuracy and fast calculations methods is beneficial in understanding how the system will behave. Evaluation of the feasibility of networks to cope with new customers, as well as the ability to determine the capacity of the current network, or to forecast the response of the system to any changes in supplying nodes (decrease/outage), as well as designing new gas network systems are areas where modelling can be used. The differential equations are used depending on the working pressure of the network and the operational conditions [19]. The major equations that must be satisfied in modelling gas flow in a pipeline are: According to forces' directions in a pipe, the governing equation has been extracted, which is shown in Figure 1 and from Newton's second law of motion-momentum Equation (3). The common equation for the steady-state flow of gas through a pipe is derived from Bernoulli's equation of fluid flow, including changing density of gas due to the pressure drop along the pipe in the direction of flow. The equation of mass conservation from inlet Point (1) to outlet Point (2) is shown in Equation (2).
where ρ is density and ω represents flow velocity.
∂ρ ∂t Newton's second law of motion (conservation of momentum) in Equation (4): where the following terms describe each portion of Equation (4): : inertia force (acting against the flow direction through the pipe) The variable pressure, p, and displacement, x, have been changed to p + dp and x + dx, respectively, which have been identified in Bernoulli's Equation (5): where "dh f " is the head losses, i and j are sending and receiving nodes, respectively, and "z" is pipe elevation. The first law of thermodynamics expresses the conservation of energy and for a system may be written as: where "Ω" is the heat added to the system, "W" is the work done by the system, and "E" is the change in energy of the system. Energy associated with the mass of the system is customarily separated into three parts: where "U" is the internal energy associated with molecular and atomic behaviour, " 1 2 mw 2 " the kinetic energy, and "mgz" the potential energy associated with the height. If the rate of work dW dt on the flow is zero, Equation (6) can be written as [46]: The final energy equation for gas flow in a pipeline can be expressed in two cases: firstly, an adiabatic flow (dΩ = 0); secondly, an isothermal flow. In adiabatic flow, the rate of heat change will be zero, and in the isothermal flow, temperature will remain constant (T = constant and dΩ = 0) [47]. The change of temperature within the gas due to heat conduction between the pipe and the ground (if it is constructed as an underground pipeline) is sufficiently slow to be neglected [48]. This means, due to the slow flowing gas, it has sufficient time to exchange heat with the ground and its temperature with the ground temperature (for underground pipes) and for aboveground pipes, with the surrounding temperature. Therefore, temperature changes can be neglected and assumed a constant temperature equal to the ground/environment temperature [19,23,24].
where "ρ n " and "Q" are the density at standard conditions and flow rate of the NG. "A" and "D" are the cross-sectional area and inner diameter of the pipe, and "p" in Equation (9) is the average pressure between two nodes. "η t " is the efficiency of pipe friction factor to convert theoretical friction to actual friction factor (Equation (10)).
Another significant and effective term in a pipe network equation is the compressibility factor, "Z". This is a correction factor in the state of gas equations for real gases, which accounts for the deviation from ideal gas behaviour. "Z" depends on pressure and temperature of gas and can be calculated using correlations such as Peng-Robinson, Redlich-Kwong, Soave-RK, the cubic equation, etc. [49]. In this research "Z" has been approximated by the PAPAY equation, which is applicable for high-pressure networks [50].
Due to slow changes of the variable in gas network pipelines, the focus of this study is a steady-state model. At steady state, outlet flow equals inlet flow, and the first term of Equation (9) is neglected. The equation then becomes: By integrating from Equation (12), the final steady-state equation is presented in Equation (13), where pressure drop depends on inertia force, frictional force of the pipe, and different elevations of the pipeline in a gas network. where: In the next section only the modelling of the Irish gas network with 13 nodes and 14 branches is explained step-by-step using the matrix approach, see Figure 2. Due to massive matrices associated with the larger Network 2, details are shown in Appendix D. Boundary conditions and network details can be seen in Table 1 and Table A1.
The following Theorem is provided.

Numerical Example Results
Both of the gas networks modelled include three supplier nodes of natural gas at a 70-barg pressure set point, which are shown in Figure 2 and Figure A1 (Nodes 1, 2, and 3). No constraint was placed on the flow of gas, which can be supplied from any of the supply nodes. Details and the boundary conditions of Network 1 and Network 2 are presented in Table 1 and Appendix D. The load at the demand nodes was weighted according to approximated town populations at the given node. These networks were then modelled as isothermal systems at 288.15 K using MATLAB (modelling the pipe flow equation using the Newton-Raphson method), SAINT (an energy system modelling software), and the new approach proposed in this paper. Results for Network 1 are presented in Tables 3 and 4, while those for Network 2 can be found in Appendix D. In both networks, there are three supply nodes, which represent the current sources of NG for the Island of Ireland Node 1, Moffat the Irish-Great Britain Interconnector, Node 2 the Corrib Gas field, and Node 3 the Kinsale Gas Field. Using the estimated gas consumption for each demand node, along with reported NG supplies [4], with the ZRT term representing the isothermal squared speed of sound in a pipeline, which was nearly 341 m/s in this case, the model was solved for pressure at each node and flow rate in each pipe. Using Equations (14) and (15), the parameters "a ij and b ij " for each pipe have been calculated, for which a ij relates to the friction factor, pipe length, ZRT, normal density of gas, and diameter, while b ij relates to the pipe angle, pipe length, and ZRT. The numbers shown in Table 2 for the smaller network and for the bigger network can be found in the Appendix D. In the pipe flow Equation (13), b ij is a secondary important term.
To supply the estimated overall demand of Network 1, the Moffat interconnectors supply nearly 53.5 m 3 /s, and gas flows from Nodes 2 and 3 are 38 m 3 /s and 88 m 3 /s, respectively. Based on the picture of Ireland's gas network (Figure 2), distributing this flow among the demand nodes resulted in an maximum pressure drop of 4 barg at Node 10 and and minimum pressure drop of 0.6 barg at Node 4 when the inlet pressure was 70 barg using SAINT. It is important to note that the maximum and minimum pressures did not occur at the minimum and maximum elevations, which was 130 m at Node 1 and 25 m at Node 11, respectively; therefore, the major source of pressure loss in the system was due to gas flow in the pipelines.  Tables 3 and 4, and the results of Network 2 are shown in Appendix D. The pipelines' and nodes' details including pipe diameter, length, and inclination are shown. In Figure 2, the designed network is illustrated.

The Novel Approach Results and Discussion
By applying Theorem 1 (see Appendix C), the following outputs are shown in Tables 3 and 4 for the results generated for Network 1. The outputs' tables of the smaller network are presented in this section; however, all outputs for the bigger gas network are shown in Appendix D. The two variables of pressure and flow rate are the main outputs from the three models and used for the benchmarking of the proposed method. As can be observed from Tables 3 and 4, there was a slight deviation between the results generated using the MATLAB and SAINT model for Network 1. Applying initial approximations for the non-linear part (F(X)) of linearized Equation (16) provided the outputs of the novel method of pressure and flow rate, which can be compared with the results from the SAINT and MATLAB models and the novel approach. As an example, the pressures at Node 5 from different solution methods were 68.32, 68.28, and 67.93 barg from SAINT, MATLAB, and the novel approach, respectively. The comparison between flow rates from the numerical methods and the novel method is shown in Table 4. As a result, The explored flow rates of Pipe Number 5 were 16.38, 16.73, and 16.22 sm 3 /s using SAINT, MATLAB, and the novel approach, respectively.
The comparison of outputs between the new method calculations and numerical solution using MATLAB for Network 1 showed 0.32% deviation of pressures and 1.38% deviation between flow rates in Table 5 and 0.3% and 2.95% average deviations of pressure and flow rate from the outputs of SAINT.
Furthermore, for the bigger network with 109 nodes, the average deviations of pressure and flow rate were 0.07% and 1.23%, respectively, using MATLAB, and the average deviations between the novel approach and SAINT software for pressure and flow rates were 0.77% and 3.48%. The comparisons are shown in Appendix D, Table A2. Along with accuracy, the time taken for each of the simulations was measured for both modelled networks. The proposed new method was significantly faster, taking 0.31 s, while 1 s was taken for the solution in SAINT and 2.80 for the MATLAB solution of the bigger network. This time is even less for the smaller network, with 0.058 s, shown in the Table 6. This method can be used for much bigger networks at continental scales with fast calculations.

Conclusions
Currently, there are several numerical methods along with commercially-available software to solve the non-linear equations describing the compressible gas flow through the gas network. In this paper, a novel, unique mathematical method has been developed for simulating a steady-state non-linear gas network to quantify the main characteristics, pressure and flow rate through the transmission pipeline system.
Using two approximate models of the Irish natural gas transmission network as a case study the novel method was compared and benchmarked against a numerical modelling solution using the Newton-Raphson method carried out in MATLAB and the commercial software SAINT. The comparison of the results from all solutions showed that the new model is much faster than the alternative methods, and the average deviation between the methods was less than 0.31% (0.3% and 0.32%) for pressure and less than 2.1% (2.95% and 1.38%) for flow rate. Based on these key findings, the new method can be applied as a reliable fast method to model various conditions in a gas network even when there is a large number of nodes and branches. The method also lends its self ideally to the modelling of integrated energy systems as it is fast, simple, and can be easily integrated with linear models of electrical networks. The networks presented here only contain boundary conditions associated with pressure for supply nodes and load for demand nodes; in reality, there may be other real-world constraints that impact gas supply or demand. These constraints need to be included when modelling real-world systems, so that in the future, the results from the novel method can be validated with accurate data from a real gas network.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix C. Matrices and Tools for the Numerical Example
Substituting the data in the model as described in the matrix equations in Theorem 1, we have: With the above result, we arrive at the results in Table 5.  Figure A1. The 109-node map of Ireland's gas network.