NUMERICAL SIMULATION OF NATURAL GAS PIPELINE TRANSIENTS

Simulations of natural gas pipeline transients provide an insight into a pipeline capacity to deliver gas to consumers or to accumulate gas from source wells during various abnormal conditions and under variable consumption rates. This information is used for the control of gas pressure and for planning repairs in a timely manner. Therefore, a numerical model and a computer code have been developed for the simulation of natural gas transients in pipelines. The developed approach is validated by simulations of test cases from the open literature. Detailed analyses of both slow and fast gas flow transients are presented. Afterward, the code is applied to the simulation of transients in a long natural gas transmission pipeline. The simulated scenarios cover common operating conditions and abrupt disturbances. The simulations of the abnormal conditions show a significant accumulation capacity and inertia of the gas within the pipeline, which enables gas packing and consumers supply during the day time period. Since the numerical results are obtained under isothermal gas transient conditions, an analytical method for the evaluation of the difference between isothermal and nonisothermal predictions is derived. It is concluded that the nonisothermal transient effects can be neglected in engineering predictions of natural gas packing in long pipelines during several hours. The prescribed isothermal temperature should be a few degrees higher than the soil temperature due to the heat generation by friction on the pipelines wall and heat transfer from the gas to the surrounding soil.


Introduction 1.1 General
Over the past couple of centuries, fossil fuels, as primary energy sources, have been essential for global economic growth. During the industrial revolution in Europe in the 19 th century, coal played a key role in supporting technological progress in agriculture, manufacturing and transport. Since then, petroleum has superseded the position of coal, and is an essential factor in sustaining our very expensive and 'dangerous' lifestyle.
Nowadays, however, the need for the cleaner fuels usage with lower content of carbon, as well as proven sufficient reserves and more stable prices than in case of oil market prices lead to the strong increase of the natural gas usage. The exploitable reserves of the natural gas are enough for the consumption in longer future time period, the carbon emission during the combustion of natural gas is approximately half of the emission by coal combustion and its price is more stable than in case of oil [1].
The natural gas is used in various sectors of industry, both as a fuel and as a raw material.
As a fuel it is used in boilers and furnaces to generate steam, heat water or to provide heat for technological purposes. It is a raw material in petrochemical manufacturing, in polymer manufacturing and used to produce hydrogen, sulphur, carbon black, ammonia, and ethylene. In domestic sector natural gas is a fuel for district and individual building heating, for cooking and sanitary water preparation.
In contrast to petroleum or coal, natural gas can be used directly as а source of primary energy that causes less carbon dioxide and nitrogen oxide emissions (greenhouse gases). Besides substantially lower carbon dioxide emissions in comparison to usage of coal and oil, the combustion of natural gas leads to negligible emissions of sulfur dioxide, as well as lower nitrous oxide emissions. All these characteristics provide benefits such as elimination of acid rains, and reduced ozone layer depletion and effects of the greenhouse gasses in the atmosphere. In addition, it can be safely transported, stored and used [2].
Hence, the current position of natural gas as a primary non-renewable energy source (second to oil in OECD 1 countries) leads to the conclusion that the analysis, design and improvement of its processes, including transportation, play a significant role for both private and public sectors while offering a number of challenges to the scientific research community [1].

Natural gas and its transmission: history and present
The natural gas is known to humans since ancient times in the Middle East. At the beginning people had been aware of burning springs of natural gas. In Persia, Greece, or India, temples were built around these eternal flames for religious practices [2]. There is also historical evidence from the ancient times that people had started to harness natural gas springs with the aim of providing their living needs. Some 900 years B.C. the drilling of the ground was applied in China with the aim of obtaining springs of gas and that gas was used as a fuel for efficient provision of their living needs. Namely, the seawater was evaporated by natural gas combustion in order to obtain salt and drinkable water. In addition, by the first century, the Chinese had developed "an advanced techniques for tapping underground reservoirs of natural gas, which allowed them to drill wells as deep as 1,460 m in soft soil; they used metal drilling bits inserted through sections of hollowed-out bamboo pipes to reach the gas and bring it to the surface" [2].
The Romans were also aware about natural gas existence. It is supposed that Julius Caesar saw a "burning spring" near Grenoble in France. Also, there is evidence that religious temples in early Russia were built around burning sources of natural gas in the ground, which represented some kind of "eternal flames" [3].
The natural gas was discovered in Great Britain in 1659, but its commercial usage started more than a century later in 1790. A source of natural gas was discovered in Fredonia in United States US in 1821, as bubbles that rose to the surface from a creek. The first natural gas well in North America was dug by William Hart, who is called as "America's father of natural gas" [2]. He applied hollowed logs for the transport of gas from the well to a nearby building and the gas was burned for illumination. In 1865, the Fredonia Gas, Light, and Waterworks Company became the first natural gas company in the United States. The first transmission natural gas pipeline was built in 1872. It was some 40 km long and it supplied gas from the wells to the city of Rochester in New 1 Organization for Economic Cooperation and Development York. This pipeline was also built of hollowed logs. In 1885 Robert Bunsen developed so called "Bunsen burner", which enabled the usage of natural gas for heating and cooking, besides its use for lighting. Certainly, at the beginning of these commercial natural gas consumptions, an obstacle of its wider usage was the lack of pipeline infrastructure for natural gas transport and distribution. In addition, a need for facilities for gas storage was encountered [2].
Further development of technology related to natural gas usage led to the exploitation of a high-pressure gas deposit in central Indiana, which started in 1891. This gas was transported for consumption to Chicago in Illinois and a 192 km long pipeline was built for that purpose. Natural gas is also extracted together with oil from the oil wells, but during the early period of oil exploitation it was observed as burden. Hence, natural gas was leaked directly to the atmosphere at the oil fields or it was burnt and the flame illuminated the oil fields day and night. Oil companies realized that this is an unreasonable practice and they started to develop gas transmission pipelines and pipeline networks for gas distribution to the consumers in large cities. This activity was an additional source of profit for them. The technological progress after the World War II boosted the natural gas consumption, for example in pipeline manufacturing, metallurgy and welding. Gas transport companies started building and expanding their pipeline systems. The fast and steady growth of gas industry finally entailed the construction of various gas facilities, including processing and storage plants, as well as a number of sustainable projects around the world since the late 20 th century. In this way natural gas became an attractive alternative to electricity and coal [1].
Despite periodic economic and international crises, new oil and gas pipelines are being  Natural gas consumption rose by 96 billion cubic metres (bcm), or 3%, the fastest since 2010. This consumption growth was driven by China (31 bcm), the Middle East (28 bcm) and Europe (26 bcm)" [5].
"Global natural gas production increased by 131 bcm, or 4%, almost double the 10-year average growth rate. Russian growth was the largest at 46 bcm, followed by Iran (21 bcm). Gas trade expanded by 63 bcm, or 6.2%, with growth in LNG outpacing growth in pipeline trade. The increase in gas exports was driven largely by Australian and US LNG (up by 17 and 13 bcm respectively), and Russian pipeline exports (15 bcm)" [5].
"2017 was a bumper year for natural gas, with consumption (3.0%, 96 bcm) and production (4.0%, 131 bcm) both increasing at their fastest rates since the immediate aftermath of the financial crises. The growth in consumption was led by Asia, with particularly strong growth in China (15.1%, 31 bcm), supported by increases in the Middle East (Iran 6.8%, 13 bcm) and Europe. The growth in consumption was more than matched by increasing production, particularly in Russia (8.2%, 46 bcm), supported by Iran (10.5%, 21 bcm), Australia (18%, 17 bcm) and China (8.5%, 11 bcm)" [5]. Natural gas is foreseen as the fuel source with the highest increase in consumption in the near future. Huge projects of transmission pipelines are planned and conducted with the aim of transporting gas from distant gas fields with great reserves to industrial areas and big cities. Natural gas is transported through long distance pipelines by work of a series of compressor stations.

Natural gas origin and composition
"Natural gas exists in nature under pressure in rock reservoirs in the Earth's crust, either in conjunction with and dissolved in heavier hydrocarbons and water or by itself" [2]. It is exploited alone from the natural cavities or porous sediments or together with crude oil. "Natural gas has been formed by the degradation of organic matter accumulate in the past millions of years. Two mechanisms (biogenic and thermogenic) are responsible for this degradation" [2].
Natural gas is composed mainly of methane. Other ingredients are paraffinic hydrocarbons such as ethane, propane, and butane. Natural gas contains nitrogen as well as carbon dioxide and hydrogen sulfide [2]. A minor amount of argon, hydrogen, and helium may exist in it. Natural gas from geographically separated areas can have substantially different composition. Table (1.1) illustrates the typical composition of natural gas. Hydrocarbons C5+ can be also included and it can be separated as a light gasoline. Some toxic substances might be present in small quantities, such as benzene, toluene, and xylenes, as well as some acid contaminants like mercaptans R-SH, carbonyl sulfide (COS), and carbon disulfide (CS2). Mercury can also be present either as a metal in vapor phase or as an organometallic compound in liquid fractions [2] .
Typical composition of natural gas is presented in Table 1.1. It should be emphasise that the gas composition can vary substantially from the values presented in Table 1.1. Standard test methods were developed for the determination of the natural gas composition and description of these methods is available elsewhere [2].

Demand for natural gas
The demand for natural gas has been steadily increasing over the last several years as shown in figure (1.2). The world consumption of natural gas in the year 2018 was 3.85 trillion cubic meters (Tm 3 ) (on the left vertical axis the consumption is presented in trillion cubic feets -TCF) [6,7]. The projected demand up to 2030 is also illustrated in the same figure.
It is difficult to predict the increase in natural gas demand in the future since it depends on several socioeconomic factors. Starting with worldwide energy demand, figure (1.3), the energy demand is expected to grow from 5.71×105 PJ (petajoule (PJ) = 1015 J) in 2010 to 9.54×105 PJ in 2050 for about 67% total increasing.   In figure (1.3), the energy demand increase is uneven across the world. In developing countries, the increase in demand is a lot higher (22% over last eight years), whereas, for industrialized countries the increase is slower (4% over last ten years). The shift in demand can have significant consequences on the demand for natural gas since transportation of gas is an important bottleneck in satisfying the demand for natural gas [8].
According to the U.S. Energy Information Administration (EIA) report (September 2019) [8] the world natural gas consumption demand will increase more than 60% from 2010 to 2050, from about 1.3×10 5 PJ to 2.1×10 5 PJ over forty years. "Natural gas use accelerates the most in countries outside of the Organization of Economic Cooperation and Development (OECD) to meet demand from increased industrial activity, natural gas-fired electricity generation, and transportation fueled by liquefied natural gas (LNG)." Natural gas consumption in non-OECD countries will grow from about 74×10 3 PJ in 2018 to around 126.5×10 3 PJ in 2050, a 71% increase. It is projected that the natural gas consumption during this time in the OECD countries will increase 17% between 2018 and 2050.
Also, the projected demand of natural gas is shown in figure (1.4). This is the most likely demand based on an assumption that fuel cell technology will not have significant contribution to transportation power. If fuel cell technology indeed becomes viable, the demand for natural gas can be even higher than predicted in figure. In figure (1.6), energy produced from the natural gas and its projected demand is show and compared with the energy provided from coal for the period of forty years. In 2010, the energy produced from coal was more than the energy from natural gas by about 3.6×103 PJ. For three years after, both fuels keep to increase. In 2013 coal shows decline in energy production and the production of natural gas continues to increase. The 2028 is the year where both coal and natural gas production is about 1.5×10 5 PJ. After this, the energy provided by natural gas will be more than the energy produced by coal. For about five years after 2028 the energy provided by coal is expected to remain constant and then start to slightly increase up to 2050. The energy produced from natural gas is expected to show a gradual increase and reach 2.1×10 5 PJ in 2050 that is about 0.2×10 5 PJ more than the energy provided by coal [8]. Hence, it is concluded that the outlook of the future natural gas roll in the primary energy mix in the World shows that its consumption, production and reserves will continue to increase for the foreseeable future [8]. The natural gas supply to the consumers is based on the following chain of technical and technological processes and activities [9].
 "Exploration: In this stage, the issue of how natural gas is found and how companies decide where to drill wells for it is addressed.
 Extraction: This stage deals with the drilling process, and how natural gas is brought from its underground reservoirs to the surface.
 Production: In this stage the processing of natural gas once is brought out from the underground takes place.
 Transport: The natural gas is transported from the processing plant to local distribution companies across a pipeline network in this stage.
 Storage: This stage accounts for the storage of natural gas.
 Distribution: In this stage, natural gas is delivered from the major pipelines to the end users.
 Marketing: This stage involves the buying/selling activity from the natural gas marketers." The reliable and efficient transportation of natural gas from production to consumption areas needs a developed transportation system. In the majority of cases the distance between the natural gas wells and consumers in industry or domestic sector is long over thousands of kilometres.
Therefore, long distance transmission pipelines are being built, accompanied with the development of complex distribution systems in urban and industrial areas with the aim of gas supply to final consumers. The supply of the natural gas is closely linked with its storage. The roll of the gas Year Energy from Natural Gas Energy from Coal storage is to adjust mainly constant gas extraction from the natural wells with variable seasonal or daily gas consumption by the final consumers. In winter periods natural gas consumption increases due to heating. Hence, the gas is accumulated in summer period usually in huge natural underground cavities, and discharged from them and consumed in winter period. On the daily level, during the reduced consumption periods, natural gas can be packed in long transmission pipelines, distribution networks and built gas storage facilities, while in later periods the gas is discharged from these storage units in order to cover peaks of increased consumption.
The whole transportation path, from the gas wells to the final consumers consists of three major types of pipelines: "the gathering system, the interstate pipeline system, and the distribution system.
The gathering system consists of low pressure, small diameter pipelines that transport raw natural gas from the wellhead to the processing plant. Natural gas from a particular well might have high sulfur and carbon dioxide contents (sour gas), a specialized sour gas gathering pipe must be installed. Sour gas is corrosive, thus its transportation from the wellhead to the sweetening plant must be done carefully" [9].

Transportation of natural gas
Natural gas is often found in places where there is no local market, such as in the many offshore fields or onshore fields in the deserts around the world. For natural gas to be available to the market it must be collected, processed, and transported.
Natural gas, as a result of the storage difficulties, needs to be transported immediately to its destination after production and processing from a reservoir. There are a number of options for transporting natural gas energy from oil and gas fields to market. These include pipelines, LNG (liquefied natural gas), MLG (medium conditioned liquefied gas), or CNG (compressed natural gas) [2].

Liquefied natural gas (LNG)
Liquefied natural gas (LNG) technology has proven to be effective over the last 30 years. In 2005, about 0.2 Tm 3 or 5.6% of natural gas was transported using LNG technology. By 2020, the worldwide demand for gas transported through LNG is expected to be 0.49 Tm 3 . The LNG was exported from eight countries (Indonesia, Malaysia, Algeria, Australia, Brunei, United Arab Emirates, United States, and Libya) and was imported by eight countries (United States, Japan, South Korea, Taiwan, Belgium, France, Spain, and Turkey).

Gas to liquid products
Gas to liquid (GTL) technology refers to the conversion of natural gas into synthetic hydrocarbon liquids, particularly middle distillates. By some estimates, 25.5 Tm 3 of natural gas are stranded too far from markets to be produced or transported profitably. This is sufficient to justify about 200 gas-to-liquid plants.
The technology of converting natural gas to liquids is not new. In the first step, natural gas is reformed and converted to hydrogen and CO. The mixture is called synthetic gas or syngas. This is the same process for converting natural gas to hydrogen, which can be used as a fuel in a fuel cell.
This step is the most expensive and consumes about 50% of the total GTL costs. In the second step, in a slurry reactor the syngas is blown over a catalyst at about 232°C and is converted to liquids. This is called Fischer-Tropsch synthesis. These liquids can be converted to other desirable products, such as synthetic fuels, using the cracking process.

Natural gas transportation via pipelines
Transportation of natural gas from gas fields with wells to consumers' areas is very important and crucial activity regarding reliability and economics of gas supply to the consumers.
Natural gas can be transported by different technological solutions, but the most economically acceptable method to transport large quantities of natural gas is by pipelines. This method of gas transport by pipelines has been boosted by metallurgical and welding techniques improvements.
Hence, there is a fast increase of pipeline networks deployment during the last decades all over the world, which enables economic gas transportation.
Pipelines can be installed both offshore and onshore, but there is an substantial difference in terms of security and construction prices. "Building pipeline systems under the sea is highly costly and technically demanding, a lot more than onshore" [9].
Transportation pipelines can be divided into three types: gathering pipelines, transmission pipelines, and distribution pipelines. Raw natural gas is transported from the production wells to the gas processing plant by gathering pipelines. Transmission pipelines transport natural gas from the gas processing plants towards storage facilities and distribution systems, while these distances can 13 be of the order of hundred or even thousands of kilometres. The transmission pipelines are under high pressure and the pressure is reduced at connections with the distribution network pipeline systems. Distribution pipeline systems can be found in communities and distribute natural gas to homes and businesses.
At the end, the natural gas pressure is further reduced in devices called regulators, which decrease the pressure to a level that is safe to enter homes or other facilities [10].
The main differences among these systems are the physical properties of the pipelines used, such as diameter, stiffness, material, etc., and the specifications of the maximum and minimum upstream and downstream pressures.
The major transportation of natural gas is carried through cross-border pipelines.
Throughout the world, major efforts are under way to increase the gathering, transmission, and distribution capacity in order to promote and support projected growth of natural gas demand [6].
The natural gas cross-border transmission pipeline infrastructure in the U.S. represents one of the largest and most complex mechanical systems in the world. This system of natural gas pipeline network is a highly integrated network that moves natural gas throughout the continental United States. More than 210 natural gas pipeline systems have about 490850 kilometres of interstate and intrastate transmission pipelines that link natural gas production areas and storage facilities with consumers. In 2017, this natural gas transportation network delivered about 0.708 trillion cubic meter (Tm 3 ) of natural gas to 75 million customers. These pipelines systems are driven by more than 1,400 compressor stations that maintain pressure on the natural gas pipeline network and assure continuous forward movement of supplies.
This system has been developed over the last 60 years, and is controlled at a very low level of sophistication [6]. Quite often, collected natural gas (raw gas) must be transported over a substantial distance in pipelines of different sizes. These pipelines vary in length between hundreds of meters to hundreds of kilometres, across undulating terrain commonly occurs because of the multicomponent nature of transmitted natural gas and its associated phase behaviour to the inevitable temperature and pressure changes that occur along the pipeline [2].

Mathematical flow modelling of gas pipelines
Optimal design of the gas transmission pipeline diameter for steady-state operational conditions is determined by the minimal overall exploitation costs. The major parts of these costs are investment cost in the pipeline and operational cost of fuel that energizes compressors. The greater pipe diameter means higher investment costs, but lower gas velocity, lower pressure drop and lower compressors' power and fuel consumption, and vice versa, the lower diameter reduces investments but increases fuel expenditures. But, the gas consumption has transient character and it also influence the optimal design of gas pipeline diameter. A design according to the maximum gas flow rate, without taking into account the possibility of gas accumulation, would lead to an uneconomical solution. Hence, there is a need for accurate prediction of operational parameters (mainly flow rates, velocities and pressure drops) of natural gas transmission pipelines both in steady-state and transient conditions. Such an engineering need has led to the development of various types of mathematical models for the prediction of gas transport. "Isothermal steady-state and transient pressure drop or flow rate calculation methods for single-phase dry gas pipelines are the most widely used and the most basic relationships in the engineering of gas delivery systems.
They also form the basis of other more complex transient flow calculations and network designs" [2].
There are several purposes that shape the demand of having precise and accurate pipeline mathematical flow models, mainly serving for pipeline balance, pressure monitoring, and deliverability. There are two types of mathematical models for lengthy pipeline flow; the steadystate and the transient models. The core difference between the two types of flow models lies in the equation of motion. In the transient flow models, there are terms that represent the change of transport parameters with time. When these terms are set to zero, the steady-state representation of the flow equation is obtained. Consequently, and due to this fundamental difference between the two types of models, the functionality of each of them differs. For purposes such as pressure monitoring and leak localization, in which the change of transport parameters with time is vital, transient flow models become a necessity. For other purposes such as pipeline design, sizing, line capacity estimations and line packing, where the changes of transport parameters with time are of no significance, steady-state models become ideal. Both of the two types are approximations of the actual conditions of the pipeline.
The prediction of the pressure drop due to gas friction on the inner pipeline wall is one of the most important tasks in the design of the gas transmission and network distribution systems. On the basis of this prediction the capacity and operational characteristics of compressor stations should be determined [11].
A mathematical modelling for the simulation of gas transport parameters is especially important for the transmission systems of large capacity due to its overall influence on the whole energy systems of regions and countries. Researchers have simulated and optimized gas pipeline networks and equipment for both steady-state and transient conditions with varying degrees of success. In Chapter 2 of this dissertation a literature review of previous research results are presented. "Historically, most of the efforts have been focused on steady-state flow conditions, but researchers have also identified the need for transient flow simulation" [12].
Fluid dynamicists and mechanical engineers are devising robust mathematical and numerical models to serve the gas transmission related purposes. Material engineers and scientists are developing advanced materials for the pipeline insulation and protection. Electric, control, and telecommunications engineers are developing sophisticated SCADA (Supervisory Control and Data Acquisition) systems in order to gain full control over the millions kilometres of gas pipelines worldwide.
A range of numerical schemes have been applied for the simulation of natural gas flow in pipelines, such as the method of characteristics, finite element methods, and explicit and implicit finite difference methods. The choice of the method is influenced more or less on the individual requirements of the system under investigation.

Problem background
In general, a mathematical model to simulate pipeline system operation, as well as the impact of design changes and equipment enhancements is urgently needed for this huge system to adjust the operation conditions. Simulation allows us to predict the behaviour of natural gas pipelines under different conditions. Such predictions can then be used to guide decisions regarding the design and operation of the real system. The control of natural gas pipelines system also requires simulation in order to obtain information about the pressure and flow rates at given points of the pipeline [13].
Natural gas driven by pressure is transported through pipeline for a hundreds or thousands of kilometres or miles (cross-border). As it flows over long distances through pipelines, energy and so pressure is lost due to both the friction of pipelines and heat transfer between the natural gas and its environment [14]. This lost pressure of the natural gas is added or recovered at the compressor stations which are installed along the route of the natural gas pipelines which consume un-neglected amount of money.
Many pipelines systems use online pressure and flow monitoring to detect leaks. In these systems, a computer algorithm compares actual pipeline operating conditions to calculated conditions. Discrepancies beyond a certain threshold are potential leaks. Numerical simulations have been used for solving these problems; these simulations are based on either transient or steadystate models of pipelines [14].

Problem statement
Simulations and analyses of natural gas pipeline transients provide an insight into flow parameters changes and a pipeline capacity to deliver gas to consumers or accumulate gas from the source wells under various abnormal conditions. This information is important in order to control gas pressure changes within acceptable minimum and maximum setpoints and to plan repairs in timely manner with the aim of sustaining gas accumulation and supply to consumers in cases of various disturbances.
The major concern of the present thesis is to devise a mathematical model and a proper solution algorithm for modelling compressible, frictional natural gas transient flow in complex pipelines to predict natural gas flow properties. The considered area of application of such model is natural gas flow in cross-border and network pipelines.
Therefore, a numerical model and a computer code have been developed for the simulation and analyses of natural gas transients, such as those typically found in high-pressure gas transmission pipelines. The model is based on the mass and momentum balance equations that describe one-dimensional, compressible, frictional natural gas transient flow, as well as on boundary conditions that enable simulation of gas flows in complex pipeline networks. The developed model is solved with the numerical procedure of the method of characteristics and implemented into the gas transient analysis (GTA) computer code.

Objective
The objective of this study is to develop a numerical model and a computer code for the simulation and analyses of one-dimensional, compressible, frictional natural gas transient flow in lengthy and shortened pipelines that are able to predict natural gas flow properties in normal and abnormal operation conditions under isothermal or non-isothermal conditions.

Thesis organization
The present thesis proposes a novel mathematical model and a computational algorithm to model the transmission of natural gas in (length/short) pipelines. A predictive numerical scheme is The natural gas transients in the long gas pipeline of the Western Libya Gas Project are presented and discussed in Chapter 6 for three scenarios related to (i) disruption of gas supply from the gas source wells to the pipeline inlet point and (ii) stopping of gas delivery to a thermal power plant and a terminal for further off-shore gas transport at transmission pipeline outlet points. The method for the evaluation of thermal effects during these gas accumulation and discharging transients are presented in this chapter, together with the comparison of pressure changes obtained with isothermal and non-isothermal model evaluations.
Chapter seven provides detailed conclusions for this thesis. The appendix section contains the derivation of formulas for mathematical modelling of natural gas one-dimensional unsteady compressible flow in pipelines based on the governing equations of one-dimensional, compressible, frictional natural gas transient flow in pipelines.

Major contribution
The major contributions of this thesis are as following: 1. Presenting chronological and technical reviews of the mathematical modelling of natural gas in pipelines. 5. Analysis of long transmission gas pipeline transients caused by disturbances of gas supply and delivery. The main aim of the performed simulations is to show the pipeline gas accumulation capacity, i.e. to accumulate gas during the disturbance at the gas delivery to consumers and to supply the consumers during disturbance in gas supply at the source.

Introduction
Natural gas pipelines systems are becoming more complex as the use of this energy source increases. Mathematical modelling is one of the most important tools used in both design and operation of natural gas pipelines. Plentiful efforts have been spent and continue to be spent on steady-state and transient mathematical models.
In recent decades, world consumption of natural gas has grown and it is now one of the most commonly used primary energy source worldwide, accounting for 24% of total primary energy supply, behind coal and oil (33.5% and 27%) respectively Simulations and analyses of natural gas pipeline transients provide an insight into flow parameters changes and a pipeline capacity to deliver gas to consumers or accumulate gas from the source wells under various abnormal conditions. This information is important in order to control gas pressure changes within acceptable minimum and maximum setpoints and to plan repairs in timely manner with the aim of sustaining gas accumulation and supply to consumers in cases of various disturbances. Therefore, a numerical models and a computer codes have been developed for the simulation and analyses of natural gas transients, such as those typically found in high-pressure gas transmission pipelines. It is noteworthy to mention that many operation conditions problems can be solved by transient flow modelling.
The most notable efforts dealing with the problem of mathematical modelling of natural gas pipelines are reported in this chapter.

Steady-state models
Steady-state gas flows have been studied in numerous research papers and theses. Steadystate flows of natural gas through pipelines are simulated in order to investigate and analyse the behaviour of gas flow in both operational and design conditions. In some researches, the investigators as Stoner [16], Mohitpour et al [17], Costa et al [18], etc., developed models which describe isothermal gas flow, and some others applied models that analyse non-isothermal gas flow as what Borujerdi and Rad [19] have done -they analysed the gas flow in pipelines subjected to wall friction and heat transfer. A few researchers presented comparisons between isothermal and non-isothermal pipelines gas flow models, as it was done by Alghlam [20]. On the other hand, some models are developed analytically such as those solved by Cameron [21], Szoplik [22], and Zhou and Adewumi [23] and the others are solved numerically as done by Mohitpour et al. [17] and some other investigators by using a various of numerical methods such as the method of characteristics, the implicit finite difference method, the explicit finite difference method, the finite elements method; the choice depends upon the particular requirements of the system under studying.
In addition, in this subsection, in order to provide some reading structure, it would be useful to find some other common characteristics among these research papers and thesis, such as solving the gas flow models with taking into account or neglecting the term of kinetic energy in the energy equation, studying the gas flow in pipelines systems as a single-phase flow and two-phase flow, and investigate the non-isothermal gas flow with and without consideration of gas wall friction, etc. fluids flow through the pipes were described by some researcher as Ouyang and Aziz [24], Rhoads [25] and Schroeder [26].
Abbaspour establishes general flow equations of simple form as a principle of the pressure loss calculations due to friction, elevation and kinetic energy [13]. Stoner had developed a new methodology for getting a steady-state solution of an integrated gas system model comprising of pipelines, compressors, control valves and storage fields. He utilized Newton-Raphson method for solving nonlinear algebraic equations [16,27].
Berard et al. developed a simulation based on computer software to perform a steady-state gas transmission network utilizing the Newton-Raphson method for solving nonlinear equations.
The simulated computer software has several features that facilitate efficient, accurate simulation of large nodal systems, including 1) optimal number of nodes, 2) implicit compressor fuel gas consumption calculation, 3) the ability to prorate gas volumes entering the network system, and 4) gas temperature distribution calculation [28].
Hoeven and Gasunie [29 ] used a linearization method to describe some mathematical aspects of gas network simulation. Tian and Adewumi [30]  Mohitpour et al. [17] addressed the significance of a dynamic simulation on pipeline transmission systems design and optimisation. The authors demonstrate in this paper that steadystate simulations are enough to optimize a pipeline when supply / demand conditions are relatively stable. In general, steady-state simulations should give a reasonable degree of confidence to the designer when the system is not subject to radical changes in mass flow rates on operating conditions. The mass flow rate varies in practice, so the most practical and general simulation is one that allows for transient behaviour. Borujerdi and Rad [19] analysed the gas flow in high pressure buried pipelines treated with wall friction and heat transfer. The governing equations for one-dimensional compressible pipe flow are derived and solved numerically. This examines the effects of friction, heat transfer from the pipeline wall and inlet temperature on several parameters such as gas pressure, temperature and mass flow rate. By using some previous numerical experiments and available experimental data, the numerical scheme and numerical solution was verified.
Zhou and Adewumi [37] presented a mathematical model describing steady-state gas flow in pipeline. The model was reduced to a second-order ODE (Ordinary Differential equations) system of first order initial-value problem with gas pressure and temperature as the two dependent variables. The fourth-order Runge-Kutta method was used to solve this ODE system.

Transient models
Some researchers developed or presented natural gas models for prediction of pressure transients, temperature transients, leakages, etc. Numerical simulations of gas pipeline network of transient flow were conducted by Osiadacz [38] to predict the gas pressure and flow rate distribution and time change within the network, with the aim of minimization of compressors' operational costs. Osiadacz used the theory of hierarchical systems to explain the dynamic optimisation of high-pressure gas networks. The author explains that mathematically the transient optimization is more complicated than the steady state simulation, but the advantage of using a dynamic simulation is that the operator can achieve greater savings. He further states that it is of great importance to be able to optimize large-scale systems represented by partial differential equations as fast as possible in order to achieve real time optimization.
Gas pressure and flow rate changes in the long transmission pipeline network under the partial reduction of the gas supply at the inlet point were predicted by Pambour et al. [39]. They applied the approach with the isothermal gas flow model, also in an investigation of long transmission pipeline transients. According to these authors, a prediction of the influence of thermal effects on the gas flow would require a good knowledge of the thermal resistance of the ground and the distribution of ground temperature, which is typically difficult to estimate. Moreover, due to the slow dynamics in transport pipelines (with gas velocity lower than 15 m/s) the flowing gas typically has sufficient time to exchange heat with the ground and adapt its temperature to ground temperature. Thus, it is reasonable to neglect the temperature changes and assume a constant temperature equal to the ground temperature, as it is done by many authors in the literature.
Also, Mohitpour et al. [17] showed that the natural gas transmission pipelines should be designed by taking into account the transient pressure changes and compressible gas accumulation in the volume of pipelines, which are caused by the daily changes of gas consumption by consumers. A calculation of the transmission pipeline diameter on the basis of an average daily gas consumption and the minimum pressure setpoint would lead to an under design of the supply capacity. A calculation of the pipeline diameter on the basis of the maximum gas consumption, but without taking into account the transient accumulation of the gas within the transmission pipelines would lead to an over design of the pipeline and increased construction costs. Zuo et al. [40] investigated gas pipeline transients as a support to the prediction of setpoints for the action of automatic line-break control valves closure, which are used to prevent the gas release in the event of a pipeline rupture accident.
On the other hand, for the temperature profile of buried gas pipelines prediction, M. Edalat As same as the steady-state flow models of natural gas pipelines, some transient models are based on isothermal flow, and the rest on non-isothermal flow. But a few of these models made a comparison between the two types of flow as was done by Osiadacz and Chaczykowski [11]. They Chebyshev (RKC) methods to solve ordinary differential equations resulting from the line approach applied to parabolic-type partial differential equations. Lewandowski [72] presented an application of an object-oriented methodology to model a network for the transmission of natural gas. For organized modeling and sensitivity analysis of dynamic systems, this approach was applied using a library of C++ classes. The model of a gas pipeline network can be formulated as a directed graph.
Each arc of this graph represents a segment of the pipeline and has associated a partial differential equation which describes the gas flow through this segment. Graph nodes corresponding to gas pipeline nodes may be categorized as: source nodes, sink nodes, passive nodes, and active nodes.   [30] determine the flow of natural gas through a pipeline system Isoth. flow Steady-state Deriving analytical equation based on mass and momentum ballance Costa et al [18] provided a steady-state gas pipeline simulation Isoth. flow Steady-state Model based on flow equation associated with the energy equation Borujerdi and Rad [19] analysed the gas flow in high pressure buried pipelines subjected to wall friction and heat transfer.

Concluding remarks
The above literature review leads to the following conclusions:  Most of the researchers focused on isothermal conditions where they neglected the effect of heat transfer from and to the gas flow in pipelines.
 Most of researchers have studied the natural gas problems in terms of steady-state and transient flow of one-phase of transported natural gas, also, in one-dimensional flow.
 Most of researchers who focused on the non-isothermal flow conditions did not take into account the heat generation due to the friction between pipe wall and the gas flows in the pipeline.  Researchers have developed numerical schemes for a flow dynamics of natural gas pipelines using different methods such as the implicit finite difference method, the explicit finite difference method, the finite elements method, and the method of characteristics.
The presented literature survey shows that there is limited information about the natural gas transmission pipelines behaviour during operational disturbances. Some of these disturbances are likely to happen during the gas pipeline exploitation, such as: (a) the stoppage of gas delivery from gas source (wells) to consumers or storage, or (b) disruption of gas consumption while the gas input from the source is available. In such cases, it is worth to know the accumulation capacity of the long transmission pipelines or a time period during which the gas pipeline accumulation capacity can satisfy consumers' needs without disruption. Regarding a need for the insight into these transient operational characteristics under likely disturbances, the topic of research of the present PhD thesis is stated.
In addition, the literature survey shows that there is a lack of a simple method for the prediction of uncertainty that is introduced by the isothermal gas flow assumption into transient gas pipeline simulations. Hence, a derivation of an original analytical method is a topic of research in the presented thesis.
The motivation of the present research is to investigate the capacity of the long natural gas transmission pipelines to deliver gas to consumers in a case of abrupt disturbance of gas supply at the inlet point. A time period is determined from the trip of the gas supply to the instance of reaching a low pressure level at the delivery points at consumers. Also, an accumulation capacity of the transmission pipeline is evaluated in cases of cease of gas delivery to consumers under sustained gas pressure at the pipeline inlet point. The results should support the operation procedures and guidelines in cases of abnormal condition operations. In order to numerically simulate the gas transmission pipeline transients, the code GTA (Gas Transient Analysis) is developed, based on the model of one-dimensional, compressible and transient natural gas flow. The model mass and momentum balance governing equations are solved with the method of characteristics, which has the potential to produce the most accurate results (Wulff, [91]). Its high accuracy originates from the fact that it reduces partial differential equations to ordinary differential equations, as well as being the only method that accurately tracks the propagation of discontinuities in first-order derivatives. The characteristic coordinates are Lagrangian coordinates for such discontinuities. An analytical method for the evaluation of the difference between isothermal and non-isothermal transient pressure predictions and non-isothermal temperature change is derived. It supports the application of the isothermal simulations by the GTA code of transient gas accumulation and discharging of the long transmission pipeline within time periods of several hours. The motivation for the evaluation of the influence of thermal effects on the pressure changes in gas transmission pipelines was also initiated by the ambiguity of previously published results.

Introduction
The In this chapter, physical properties and flow dynamic parameters of natural gas, which is treated as mixture of non-ideal gases are discussed. Then the mathematical modelling of natural gas flow in pipeline systems is presented. Finally, numerical methods for computing of these natural gas flows are described.

Gas properties
In this section the properties of natural gas that influence gas flow through a pipeline are discussed. The relationship of pressure, volume, and temperature of a natural gas is presented and how the gas properties such as density, viscosity, and compressibility change with a variation of temperature and pressure.

Density of Gas
Density (ρ) is the ratio of the mass (m) of gas and the volume (V) that the gas occupies.  Density is expressed in kg/m 3 in SI units.

Specific Gravity
Specific gravity (G) is a measure of how heavy the gas is compared to air at a particular temperature. Sometimes it is called gravity or relative density.
where, ρ g : density of gas.
ρ air : density of air.
It is noted that ρ air is the density of dry air at the temperature of 20 o C and the pressure of 101.325 kPa. In terms of the molecular weight (M), gravity of gas can be calculated as following: where, M g : molecular weight of gas.
M air : molecular weight of air. Because natural gas is formed of a mixture of several gasses (methane, ethane, etc.), molecular weight M g in equation (3-3) is referred to as the gas mixture apparent molecular weight.
where, M i : molecular weight of natural gas component i, g/mol. y i : mole fraction of natural gas component i, %

Viscosity
The viscosity of fluid (gas or liquid) represents its resistance to flow. It depends on fluid temperature and pressure. Since natural gas is a mixture of pure non-ideal gases such as methane and ethane, the following formula key rule is used to calculate the viscosity from the viscosities of component gases: Where, µ i is a dynamic viscosity of natural gas component i (kg/ms), M i is a molecular weight of natural gas component i (g/mol), and y i is a mole fraction of natural gas component.
A related quantity to the dynamic viscosity µ is the kinematic viscosity (ʋ):

Ideal gas law
The ideal gas law sometimes referred to as the perfect gas equation, states that the pressure, volume, and temperature of the gas are related as following: where, p stands for pressure, T represents temperature, R is the ideal gas constant (8.314 J/mol K), and n is a number of moles which can be calculated as:

Real gas properties
The ideal gas equation presented in section (3.2.4) can be applied when dealing with real gases, and get adequately accurate results when the pressure levels are similar or close to the atmospheric pressure. For most real gases, the ideal gas equation will not be appropriate if the pressure values are considerably higher. To achieve reasonably accurate results, ideal gas equation should be modified.
It is necessary to define two terms which are called critical temperature and critical pressure.
A real gas critical temperature is defined as the temperature above which a gas cannot be compressed to form a liquid, whatever the pressure. The critical pressure is known as the minimum pressure required for the compression of gas into a liquid at the critical temperature [92].
Real gases can be treated with a modified form of the ideal gas law described in section (3.2.4), if the modifying factor, known as the compressibility factor z is included. This factor is also called the deviation factor. It is dimensionless number less than 1 and varies with gas temperature, pressure, and gas composition.
Including the compressibility factor z, the ideal gas equation gets the following form: The ratio of the gas temperature (T) to its critical temperature (T c ) is called the reduced temperature and is defined as: The reduced pressure is the ratio of gas pressure (p) to its critical pressure (p c ) and is given by:

Natural gas composition and pseudo-critical properties
In reality, natural gas is a mixture of several gaseous components. The critical temperature and critical pressure can be found for each pure component that constitutes this mixture of gases.
However, the critical values of temperature and pressure of the gas mixture, which are called respectively the pseudo-reduced temperature (T pr ) and pseudo-reduced pressure (p pr ) need to be calculated as follows [92]: where, T pc and p pc represent pseudo-critical temperature and pseudo-critical pressure. These quantities are determined in an analogous way to one used to calculate the molecular weight.
Therefore, the apparent molecular weight is defined in equation (3-4) as following: In an analogous fashion, Kay's rule can be used as following to calculate the average pseudo-critical temperature (T pc ) and pseudo-critical pressure (p pc ) of the gas mixture: For the given mole fractions (y i ) of gas components.
In equations (3-14) and (3-15) T ci and p ci represent the critical temperature and critical pressure of a pure component i within the gas mixture.
For the case that the composition of gas mixture is not exactly known, i e. the mole fractions of the various components in the natural gas mixture are not available, the pseudo-critical properties of the gas mixture can be computed if the specific gravity (G) of gas is known in the following approximate way [2]:

Compressibility factor
As introduced in section 3.2.5, the compressibility factor is a measure of how similar real gas is to the ideal gas. The compressibility factor z is defined as the ratio of the volume of gas at a given pressure and temperature to the volume of the gas would occupy at the same temperature and pressure if it were an ideal gas. The factor z is a dimensionless number close to 1 and its value depends on the gas gravity, gas temperature, gas pressure, and the critical gas properties.
Generalized plots showing the variation of z with pseudo reduced temperature (T pr ) and pseudo reduced pressure (p pr ) can be used for most gases for calculating the compressibility factor, as shown in Besides using the chart, the compressibility factor z can also be computed. The methods for the calculation of the compressibility factor z are presented in the following. 1 where, C 1 = 5260. For C 1 value derivation see Appendix A-3. Further, p ave and T ave represent the average pressure and temperature at any location on the pipeline. Therefore, for two points along the pipeline at pressure p 1 and p 2 the average pressure is (p 1 + p 2 )/2 and average temperature (T 1 +T 2 )/2. For more accurate evaluation of average pressure the following formula can be used.
In addition, at any particular point along the pipeline; the compressibility factor is determined as follows: where, p is a gauge pressure of gas in kPa and T is in K.

Flow regimes
In high-pressure gas transmission lines with moderate to high flow rates, two types of flow regimes are normally observed, which are turbulent flow and laminar flow. A determination of whether a given flow in pipe is laminar or turbulent is necessary, since the two different flow regimes often need different methods to analyse the flow behaviour.  The regime of flow is defined by the Reynolds number, which is a dimensionless expression, which represents the ratio between the momentum forces of the flow to the viscous forces of the fluid: where, Re: Reynolds number.
D: inner diameter of pipe.
u: gas flow velocity.
Reynolds number is used to characterize the type of flow in a pipe, such as laminar, transitional, or turbulent flow. It is also used to calculate the friction factor in the pipe flow. For Reynolds numbers less than 2100 the flow in pipes is normally laminar or stable. Turbulent flow in pipes occurs when the Reynolds number is greater than 4000. For the so-called the transition region (2100 < Re < 4000) the flow may be either laminar or turbulent, depending upon factors like the entrance conditions into the pipe and the roughness of the pipe surface. In general transition region conditions should be avoided in designing piping systems. In natural gas transmission the Reynolds number is much greater than 4000 [92]. Therefore, transport of natural gas in a pipeline is typically turbulent flow.

Friction factor calculation
When gas flows in a pipeline, friction occurs between the flow stream and pipeline walls and causes pressure losses. This pressure loss is computed by introducing friction factor. The friction factor is a dimensionless parameter depending on the Reynolds number of flow and roughness of pipe walls. In engineering literature, there are two formulation of friction factor; Darcy friction factor and Fanning friction factor. The relationship between the both factors is given by: where, f d is a Darcy friction factor, and f f is a fanning friction factor.
Darcy friction factor is more general and will be used in this study. For the sake of simplicity, the Darcy friction factor hereafter will be denoted by the symbol f. Re f  (3-23) The friction factor for turbulent flow is a function of the Reynolds number and relative roughness of pipe walls (defined as the ratio of absolute wall roughness e and inside pipe diameter D). This dependence is graphically presented by Moody diagram in Figure 3.3. As shown in Figure   3.3, turbulent flow in pipes (Re > 4000) is divided into three zones; turbulent flow in smooth pipes, turbulent flow in rough pipes, and transition flow between smooth pipes and rough pipes.
The friction factor f only depends on Reynolds number for turbulent flow in smooth pipes.
For fully rough pipes, f is more dependent on relative roughness of pipe walls (e/D). The value of friction factor depends on both the roughness of pipe wall, and Reynolds number in the transition zone. As shown above the roughness plays an important role in determination of friction factor.
For that reason, in Table 3.3 typical values of absolute pipe roughness (e) are given. There are many correlations for calculating the friction factor. The most widely used ones for evaluation of friction factor in the gas flow in pipelines are presented below. Re In order to calculate the friction factor f from equation (3-24) one must use a trial and error approach. Re

 American Gas Association (AGA) correlation
American Gas Association (AGA) correlation is derived as a result of a study which dealt with determination of the transmission factor for gas pipelines. The transmission factor F is related to the friction factor f in the following way: The transmission factor F is determined using the method of two separate equations. First, F is calculated for the zone of turbulent flow in rough pipe. Next, F is determined for the zone of turbulent flow in smooth pipe. Finally, the smaller of the two values of the transmission factor is used for pressure drop calculation.
Based on these investigations, AGA suggests using the following formula for F for the fully turbulent region, based on relative roughness e/D and independent on the Reynolds number.
where, V is the volume flow rate of the natural gas, and E is pipeline efficiency.
Summary of various correlations for friction factor used in the gas pipeline industry is presented in Table 3.4. In the present thesis, for a code, it is more complicated to write all flow equations in their different forms; so, friction factor is calculated individually then substituted as an input in a code for flow calculations.

Velocity of natural gas in pipeline
Unlike a liquid pipeline, the natural gas velocity depends upon the pressure and, hence, will vary along the pipeline even if the pipeline diameter is constant, that is due to the change in compressibility of gas. In addition, if the flow is non-isothermal, the gas velocity is affected by the variation of gas flow temperature, because of its impact on the natural gas compressibility.
The highest gas velocity will be where the pressure is least and that is at the downstream end. On the opposite, the lowest value of velocity of gas will be at the upstream end, where the pressure is the highest.
Mathematically, the calculation of the velocity of the one-dimensional, compressible, frictional natural gas transient flow could be done numerically by the combination of the mass balance and momentum balance. The derivation of pressure and velocity of the natural gas transient flow in pipe will be described in detail in Chapter 4.

Equation Application
Colebrook

Erosional velocity
The velocity of natural gas flows in a pipeline is directly related to the pressure. The gas velocity increases as the flow pressure decreases. With the velocity increase, the vibration and noise occur. Another problematic issue is that higher velocities cause erosion of the pipeline during long period of time. If the gas velocity exceeds the erosional velocity calculated for the pipeline, the erosion of the wall is increased to rates that can significantly reduce the life of the pipeline.
Therefore, it is necessary to control gas velocity in natural gas transmission lines to prevent it from rising above this limit. The upper limit of the gas velocity is usually calculated approximately from the following equation [92]: where, C 2 is an empirical constant (C 2 = 122 for continuous service as per API 14E 2 )

Heat transfer consideration of gas flow in pipeline
Generally, in some applications, where pipelines are relatively short and at low pressure, an isothermal (i.e. constant temperature) assumption for the gas flow is fairly sufficient. There are certain characteristics of lengthy pipelines (e.g. cross-border pipelines) that make the implementation of an isothermal flow model inadequate. The majority of these pipelines transport massive amount of gas every day, which requires the line to be at high pressure values all along its route. The energy loss due to pressure drop is mostly caused by friction, this lost energy is transformed into heat that is dissipated in the ground. In some cases, when the pipeline routes from north to south or from east to west and vice versa, the climatic changes along the year create relatively large difference in soil temperature, which can pump the heat out from the gas reducing its pressure. For all these reasons, it is useful to include in some studied cases a heat transfer model that takes into consideration the heat transfers between gas and its surrounding (Osiadacz and Chaczykowski, 2001) [11].
Natural gas temperature in a pipeline is affected by the conductive and convective transfer of heat in a radial direction, by the accumulation of heat in the surrounding soil, and by the Joule-Thomson effect. Modelling the flow of natural gas in pipelines requires consideration of the physical processes that govern the flow. In this section, the physical laws governing the processes that take place during natural gas transportation are applied in the derivation of mathematical expressions to model natural gas flows.

Governing equations
The flow of natural gas in pipelines is governed by the time-dependent continuity and momentum for isothermal flow and continuity, momentum, and energy equations for nonisothermal flow and an equation of state for homogenous, geometrically one-dimensional flow. By solving these equations, the behaviour of gas parameters can be obtained along the pipe network [2]. Some of investigators developed the basic equations for one-dimensional unsteady compressible flow that include the effects of wall friction and heat transfer.  The equation of state is differentiated by time t and by spatial coordinates x The mass balance equation (3-35) is transformed by the introduction of derivatives (3-38) and (3-39) and the following form is obtained with the pressure material derivatives where the speed of sound is expressed as Determination of the speed of sound is presented in Appendix A-2.
Equations (3-40) and (3-36) present a set of two partial differential equations of the hyperbolic type as follows where (3-44) In this system of equations dependent variables are the pressure and velocity of fluid, and the independent variables are the time and space coordinate. In order to solve the above system of equations it is necessary to specify the appropriate initial and boundary conditions. The initial conditions are defined with flow parameters of the fluid at the initial time prior to disturbance.
Boundary conditions are defined on the basis of the state of the fluid at the inlet and outlet of the The analytical solution of this system cannot be obtained, so a numerical method is applied to determine the particular integral.

Solution methods
Various numerical schemes have been developed for the solving of the mass, momentum and energy balance equations for one-dimensional transient pipeline flow, such as the method of characteristics, the finite elements method, the explicit finite deference method, and the implicit finite difference method. The choice depends partly upon the particular requirement of the system under investigation.

Method of characteristics
The method of characteristics is a technique for solving hyperbolic partial differential equations (PDE). Typically the method applies to first-order equations, although it is valid for any hyperbolic-type PDEs. The method involves the determination of special curves, called characteristics curves, along which the PDE becomes a family of ordinary differential equations (ODE). Therefore, it can be used to transform the partial differential of the continuity, momentum and energy equations into ordinary differential equations [27]. The resulting characteristics equations are solved numerically either on a grid of characteristics or on a rectangular coordinate grid. This method has the potential to produce the most accurate results (Wulff, 1987) [91]. Its high accuracy originates from the fact that it reduces partial differential equations to ordinary differential equations, as well as being the only method that accurately tracks the propagation of discontinuities in first-order derivatives. The method of characteristics was also used for the simulation of natural

Finite element method
This method can handle some boundary conditions better than finite difference methods. On the other hand, the method has not been commonly used for gas transient flow modelling because computing time and the storage requirement are high. The element size, shape, and distribution are relatively flexible, so that nonuniform internal distribution of nodal points is possible. This method was compared with the application of the finite difference for the simulation of gas pipeline transients by Osiadacz and Yedroudj (1989) [47].

Explicit finite difference methods
There are several explicit methods of finite difference such as first-order and second-order approximations. A first-order approximation is typically not sufficiently accurate to model gas transients in a pipeline and therefore attention is focused on second-order methods [45]. The main drawback of the second-order approximation is that these techniques require a greater amount of computer time and are therefore not ideal for examining large systems or analysing unsteady flows over long periods of time.

Implicit finite difference methods
The main advantage of using an implicit method over the explicit method is that some sort

Central difference method
In this method, the partial derivatives are approximated for sections of the pipeline rather than node points. It was used by Wiley at al.
[53] to solve for the transient isothermal flow field gas pipeline network. For the non-linear equations, the Newton-Raphson method was used, and sparse matrix algebra reduced the solution time for the simultaneous equations. Although this method requires a large amount of computer storage to handle the coefficient matrix and lengthy execution times, these major disadvantages can be overcome by using a sparse matrix method.

Crank-Nicolson method
This method is a central difference solution of high-order accuracy. It was utilized by Heath and Blunt (1969) [46] to solve the conservation of mass and momentum equations for slow transients in isothermal gas flow. The main advantage of this method is that it does always give a stable solution according to the Neumann stability analysis of large time step for nonlinear problems.

Fully implicit method
Whereas the explicit finite difference methods are forward difference methods, the fully implicit method is a backward method. This method mostly is unconditionally stable. It is very robust for the gas pipeline industry because of relatively slow transient. The implicit method garantees stability for a large time step, but requires a numerical method such as the Newton-Raphson method to solve a set of nonlinear simultaneous equations at each time step.

MODEL FORMULATION AND SOLUTION ALGORITHM
In this chapter the procedure for the computation of natural gas properties and flow field variables (pressure and velocity) is presented. The newly proposed model is based on the method of characteristics.

Application of the method of characteristics for the simulation of natural gas pipeline transients
The transient one-dimensional natural gas flow in pipelines is described with the mass and momentum balance equations. These equations are partial differential equations of the hyperbolic type. In this research, the method of characteristics is used for the numerical solution of the system of partial differential equations of hyperbolic type. This method can solve the system of two quasilinear partial differential equations (3-42) to (3-43), with the two dependent variables (pressure and velocity) and the two independent variables (time and space coordinate).
The method of characteristics converts the quasilinear system of partial differential equations (3-42) and (3-43) into a system of differential equations with the total differential, wherein the family of curves is determined in the space- The dependant variables are marked as a general function and their total derivatives are . x By substituting corresponding equation (4-2) for each dependant flow parameter into (4-1) the following equation is obtained where coefficients  1 and  2 are determined from the condition that the expressions in equation (4-3) that multiply the partial derivatives of dependant variables p and u with respect to time t are zero. Hence, a system of two linear homogeneous equations is obtained of dependant variables and a system of ordinary differential equations is obtained as follows By substituting corresponding equation (4)(5)(6) in equations (4-9), and (4-10) it is obtained   By approximating the total differentials in equations (4-11) and (4-12) with finite differences along the characteristic directions, the following system of algebraic equations is obtained, Solving of algebraic equations (4)(5)(6)(7)(8)(9)(10)(11)(12)(13) and (4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14) provides expressions for the calculation of p D and Equations ( Solving equations (4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21) to  by u R , c R and p R it is obtained Solving the equations (4-30) to (4-33) by u L , c L and p L is obtained   Calculation of all dependent variables (p, u), as well as specific volume v, in the nodes of R and L enables the prediction of pressure and velocity in the node D according to equations (4-15) and (4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16). Specific volume of fluid in the node D is determined from the equation of state, and the local speed of sound in the fluid can be determined from the appropriate theoretical expressions or an empirical correlation for the speed of sound.
In the presented method, the spatial step, i.e. the distance between nodes, is constant. The time steps is determined by Courant criterion that provides the stability of numerical solutions min , 1, 2, ... , , wherein the minimum time step, for a given value of spatial step ∆ , is determined by the maximum value of the sum of the speed of sound and the absolute value of the fluid velocity.

Boundary conditions
Boundary conditions are defined for the pipe inlet and outlet. In case of the pipe inlet the C + characteristic path in (Fig. 4.1) and corresponding characteristic equation condition inside a pipe network is a junction of two or more pipes (Fig. 4.2) and it is derived as follows. The pipes that transport fluid towards the junction node D in (Fig. 4.2) are denoted with J i , while the pipes that transport fluid from the node D are denoted with I j . The characteristic equations in J i pipes can be written for C + paths from point R to node D ( Fig. 4.1) as follows , , , , 1, 2, ...
The characteristic equations in I j pipes can be written for Cpaths from point L to node D ( Fig. 1) ,

Flowchart of the calculation process
The calculation flowchart of the GTA code is developed in a way to enable defining the pipe network and appropriate boundary conditions by input parameters. The flowchart is shown in  The transient was simulated with the GTA code that is developed in this presented thesis.
The pipeline was discretized with 371 nodes. Further grid refinement by increasing the number of nodes provided practically the same results. The maximum gas flow rate at the pipeline outlet is reached after 13 hours, while the minimum pressure is reached after 15 hours. Again, the delay of minimum pressure occurrence after the maximum flow rate at the pipeline outlet is attributed to the gas accumulation in the pipeline and inertia of the gas mass along the pipeline. After 18.7 hours the gas flow rate at the pipeline outlet remains constant (Fig. 5.3). A certain discrepancy between measured and calculated data is shown.
The measured maximum pressure is higher for approximately 0.1 MPa than the calculated values after 8 hours (Fig. 5.4)

Case 3
The ability of the GTA code to predict transients in gas pipeline networks is validated by a simulation of transient in the gas network shown in figure 5.5.
Dimensions of three pipelines that form the network are presented in Table 5.1. The gas specific gravity is 0.6, the operational temperature is 278 K, and the friction factor is considered to be constant and equal to 0.003.

Conclusion remarks
The GTA code is validated by computer simulations of transient cases reported in the literature.
The simulated cases include transients caused by the variable gas consumption and boundary pressure pulses.
It is shown that the calculation procedure is numerically stable and the good agreement is obtained between the GTA code results and the previous published results. The presented model derivation and analysis of validation results show that the applied method is relatively easily implemented in the computer code, the calculation procedure is robust and the reliable simulations are obtained for both slow and fast gas pipeline transients.

TRANSIENT BEHAVIOR OF A LONG TRANSMISSION GAS PIPELINE
The GTA code was used to analyse the behaviour of natural gas transient flow in long transmission pipeline. Real natural gas transmission pipeline in Libya was taken for studying, and some scenarios were assumed and simulated to predict the gas flow parameters to investigate its behaviour. Also here, the presented results are published in reference [96].

Project
The developed GTA code was applied to the analysis of transients in the onshore gas    In the later period between 1 and 8 hours both calculated and measured values show the pressure decrease. In the period between 8 and 9 hours the gas flow rate decreases at the outlet in the Mellitah Complex (Fig.6.3) and both measured and calculated values show the pressure increase ( Fig.6.4) due to this gas flow rate change. In the period between 10 and 11 hours the gas flow rate at the outlet in the Mellitah Complex increases (Fig.6.3) and this leads to the pressure decrease as shown in Fig.6.4 by both measured and calculated values. The maximum difference between these values is lower than 0.02 MPa and the calculated pressure transient behaviour is in the complete agreement with measured behaviour in the period when the influence of the uncertainty of the initial condition is diminished.
Further work was directed towards investigation of gas pipeline transport capacity in transients caused by a trip of gas source at Wafa Desert Plant and by a trip of gas delivery at the TPP and the Mellitah Complex.

Scenario 1
The trip of the gas supply in Wafa Desert Plant is assumed. As presented in figure 6.5 the gas supply in Wafa Desert Plant is constant for 2 hours and then suddenly stops. The gas delivery in the Mellitah Complex and to the TPP is kept constant at the initial level that corresponds to the nominal operation. These flow rates are specified boundary conditions for this simulation. Calculated pressure values are shown in figure 6.6. The pressure at the main pipeline inlet at Wafa Desert Plant is denoted as (Pipe 1in) as in figure 6.6. The pressure in the junction D (Fig.6.1) equals the pressure at the outlet of Pipe 1 and inlets of Pipe 2 and Pipe 3, as presented in figure 6.6 (Pipe1out = Pipe2in = Pipe3in). The pressures at the outlet in the Mellitah Complex and at the outlet in the TPP are denoted as (Pipe 2out and Pipe 3out) in figure 6.6. All these pressure values are constant for the first 2 hours till the gas supply trip. decreasing flow rate at the inlet of pipeline 2 (denoted as (Pipe2in) as in Fig.6.7) still exists in the long period of 11 hours after the outlet flow stoppage due to the pressure increase and the gas accumulation.
where the product of specific heat capacity at constant pressure p c and temperature T represents enthalpy and k is the heat transfer coefficient from the gas to the surrounding soil at temperature s T . The first term on the right hand side of equation (6-1) represents the heat generation due to the friction between the flowing gas and the inner pipeline wall. It is noted that the heat generation due to the wall friction is of the order of MW in long transmission gas pipelines. In case of the Western The second term is the heat transfer rate from the gas stream to the surrounding. Differential equation (6-1) is solved analytically by applying the following relations and assumptions: a) The product of density and velocity u  is constant under a steady-state condition.
b) The heat transfer coefficient is determined by the heat conduction from the pipeline outer surface through the soil.
The heat transfer rate per unit length of the buried pipeline is calculated as [100].
, where x is the depth from the ground surface to the centerline of the buried pipeline. The soil temperature in the massive of the ground is T s , T is the gas temperature in the pipeline and  is the soil thermal conductivity. The relation between surface and linear heat flux is and introducing equation , it follows c) The gas velocity changes along the pipeline due to the pressure, temperature and consecutive density change.
d) The soil temperature depends on the ground surface temperature change, which is determined by the seasonal and day-night period changes, and on the soil conductivity. The soil conductivity changes along the pipeline, especially in cases of hundreds of kilometres long pipelines. The precise information about the soil characteristic is usually not available. Further, the soil temperature at some distance from the ground surface changes slowly with time and usually it can be assumed constant during a 24 hours day period [101]. According to the above presented analyses, the parameters (u), u, c p , f, k and T s are approximated fairly well with constant values. Therefore, Eq. (6-1) is solved analytically in the following form Friction factor and compressibility factor were calculated respectively with Colebrook-White equation and California Natural Gas Association (CNGA) method [14].
The thermal effect is evaluated first by the introduction of the above defined parameters and the value of the heat conduction coefficient = 0.64 W/(mK) into equation   The calculation of natural gas temperature changes along the pipeline with equations  and  according to equation  are presented in figure 6.13 for the heat conduction coefficient values = 0.64 W/(mK) and 1.28 W/(mK). As shown, the gas temperature decreases within the first hundred kilometres even with the low value of the heat conduction coefficient related to the dry send. For higher values of , which are most common, the natural gas temperature decrease at the pipeline inlet part will be more intensive. The difference between gas temperatures calculated with The temperature change in figure (6.13) is presented till 370 km since at that distance the gas mass flow rate is reduced in the main transmission pipeline due to the branch towards the TPP, which leads to the further decrease of difference between gas and soil temperatures. It is also noted that after 200 km the gas temperature is practically constant in case of = 1.28 W/(mK).
According to these results, the conclusion can be derived that the difference in the pressure change calculation with an isothermal and a non-isothermal model is small, as follows. Since the natural gas density change is about 5% with the temperature change form 315 K to 299 K (the density change with temperature is related to (315/299 = 1.05), the difference in the pressure drop calculated with the non-isothermal and isothermal models is lower than 5% (assuming that the adopted isothermal gas temperature is between the maximum value of 315 K and the minimum value of 299 K, and according to the well-known Darcy relation that ). This   Energy Balance Total enthalpy H is expressed as The total enthalpy H in Eq. (6-8) is replaced with the product of fluid mass M and specific enthalpy h and the specific enthalpy is expressed as the product of the specific heat capacity c p and temperature T. After derivation of the left hand side term in Eq. (6-8), by taking into account the mass balance Eq. (6-7) and by assuming that the change of the specific heat capacity by pressure and temperature is negligible (c p =const. in the range of pressure and temperature change during the analysed transients) the following expression is obtained for the temperature change.
Differentiation of equation (6-9) gives, and from equation (6-14) is derived Substation of equations (6-7) and (6-13) into equation (6)(7)(8)(9)(10)(11)(12)(13)(14)(15) gives The heat ratėis determined as The first term on the right hand side represents the heat generation due to the gas friction on the pipeline wall, while the second term is the heat transfer rate from the gas stream to the surrounding and the linear heat transfer rate is determined with Eq. (6-2).
As an assumption, in case of the trip of gas delivery the gas flow rate at the pipeline outlet is ; hence, equations (6-13) and (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16) The compressibility factor z is function of p and T and its derivative is Dividing the above equation with d and assuming the isothermal conditions, i.e. 0 dT  , the following expression is obtained Introduction of Eq. (6-24) into Eq. (6-22) leads to By expressing the ideal gas law in the form T zR p g    and differentiation by temperature for isobaric conditions leads to Finally, Eq. (6-19) is written as Substitution of Eq. (6-29) into Eq. (6-18) leads to 2 2 dT a b T dt   (6-32) Parameters on the right hand side of Eqs. (6-30,6-31) and (6-33,6-34) are taken as constant for a certain range of gas pressure and temperature change and gas packing with the constant gas inlet mass flow rate and temperature. This assumption enables an analytical solving of the differential equations (6-29) and (6-32). Equation (6-32) is solved as Substitution of Eq. (6-35) into Eq.  leads to the following solution In order to evaluate the temperature and pressure change during the gas packing of Scenario 3, the stated balance equations (6-7) and ( The same is concluded for the case when the gas inflow is stopped (̇= 0 ) and its delivery to the consumers is continued with unchanged flow rate, such as in Scenario 1 applied to the Western Libya Gas Project in Section 6.1. In this case the temperature and pressure change differential equations (6)(7)(8)(9)(10)(11)(12)(13) and (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16) have the form During the pipeline discharge transient the gas outlet temperature T out is not constant.
Therefore, in the above equations the outlet temperature is approximated with the mean gas temperature. It should be noted that this approximation leads to an even more conservative approach to the estimation of the temperature and pressure change. Namely, the term The differences of mean temperature and pressure changes under non-isothermal and isothermal conditions, in case of the pipeline empting Scenario 1 in Section 6.1, when the gas supply at the Wafa plant is stopped and the gas delivery is continued with the value of the initial flow rate, are respectively -0.8 K and 0.17 MPa (this pressure difference is 11% of the calculated pressure change under non-isothermal conditions).

Conclusions
The shows that the required delivery flow rates at the TPP and at the Mellitah Complex are maintained for a period even longer than 24 hours, due to the gas accumulation in the large volume of the gas transmission pipeline.
-Scenario 2 shows that even though the delivery flow rate in the Mellitah Complex is stopped after 1 hour, the decreasing flow rate at the inlet of the long pipeline towards Mellitah Complex still exists in the long period of 11 hours after the outlet flow stoppage due to the corresponding pressure increase and the gas accumulation. After one hour of the trip of gas delivery in the Mellitah Complex, the pressure along the several hundred kilometers long pipeline towards the Mellitah Complex is practically constant. The pressure increase within the whole pipeline system indicates gas accumulation.
-Scenario 3 shows that although the whole gas delivery to consumers (Mellitah Complex and TPP) is stopped, there is still a gas inflow at the main pipeline inlet point for a period about 11 hours, due to the gas accumulation and corresponding pressure increase. Because of its relatively short length of 5 km, the flow rate in the branch towards the TPP practically immediately stops. The flow rate sustains for a time period longer than 11 hours in the main transmission pipeline due to its long length of 525 km. -The gas temperature in steady-state condition is determined by the heat generation due to the gas friction on the pipeline's wall and by the heat transfer from the pipeline to the surrounding ambient, as presented by equation . The heat generation by friction in the long transmission pipelines is of the order of MW and according to equation (6-37) there is a difference of the gas and soil temperatures by a few Celsius degrees. This difference should be taken into account also in case of isothermal calculations (in Subsection 6.2.1 the adopted soil temperature is 295 K, the gas inlet temperature is 315 K and the calculated mean gas temperature in steady-state operation is 301.5 K). Hence, the assumption that in the long pipeline the gas temperature is equal to the surrounding soil temperature is not adequate. The gas temperature is a few degrees higher and, as explained, it is determined by the heat generation by friction and its transfer from the gas to the surrounding soil.
The developed GTA code and presented results are a support to planning and specification of operational and repair procedures and guidelines in cases of abnormal conditions.

 Continuity Equation
In figure A1 below, the control volume of the continuity is shown, where the conservation of mass can be written as follows:

Figure A1
Control volume of continuity equation For the control volume illustrated in figure A2, and using the following force component summation, the momentum equation can be written as below:

  
Dividing by dx and the cross section area A leads to where:̇is the heat transfer per unit volume, W/m 3 , e is the internal energy per unit mass in J/kg and Q is the heat transfer in W. The separation of the second term leads to

From momentum equation (A-4) and multiplying by u it is obtained
where w is a work of frictional force per unit length of pipe The above equation is written as On the other hand, it is known that = , so, we can consider the following equation:  o Pressure history in the gas pipeline of the Western Libya Gas Project during the gas supply trip (Fig. 6