NUMERICAL SIMULATION OF THE TRANSIENT FLOW IN NATURAL GAS TRANSMISSION LINES USING A COMPUTATIONAL FLUID DYNAMIC METHOD

Until now, almost all tools that are used to analyz e gas transportation networks are based on steady f low conditions. Some design problems can be solved by t he assumption of a steady state. However, some circumstances force us to use transient analyses, w hich commonly include changes in customer demand, t he failure of a compressor station or gas pipeline rup tures. In this study, the dynamic behavior of natur l gas in transmission lines was investigated through Computa tion l Fluid Dynamic (CFD) simulations. The system was modeled using the equations of continuity, moti on (or momentum), energy and the ideal gas equation when considering turbulence effects in a two dimens ional cylindrical lattice. The coupled partial diff erential equations were discretized based on the finite volu me method. The accuracy of this method was verified by comparing the experimental field data with this app roach, showing errors of approximately 4 to 4.5%, which shows the precision of this approach. Finally , an important case study was used as an applicatio n for the model and the best possible operating solution was proposed for a compressor station failure durin g peak demand times when the pipelines were operating the maximum flow rate.


INTRODUCTION
The analysis of flows and pressure drops in pipe systems has been investigated by many different researchers and has usually been based on the assumption of steady state conditions (Chaczykowski et al., 2011;Nouri-Borujerdi and Ziaei-Rad, 2009;Osiadacz and Chaczykowski, 2001).Transient flows are generated in gas systems by a time dependency of the load at different distribution points and from adjustments made by system operators in response to the demands of the system.As larger loads are placed on existing systems, the ability to supply gas at contracted pressures during peak customer demand periods is reduced.To cope with this situation, a larger variation of the supply system input capacity must be provided, resulting in more severe transient flow effects and requiring that a careful watch be maintained to successfully operate the system.Dynamic effects are particularly important when calculating the effects of short term emergencies such as a temporary reduction in the supply resulting from equipment failure (Chaczykowski, 2010;Dorao and Fernandino, 2011;Majid et al., 2012).
There have been many studies on the flow of compressible gases in pipelines published in textbooks, papers and technical documents.Adiabatic flows and isothermal flows are two limiting cases that are often considered (Churchill, 1980;Crane, 1988;Crowl and Louvar, 2002;Farina, 1997;Levenspiel, 1998).Adiabatic flow conditions assume a flow through an insulated pipe.These conditions are usually valid for short pipelines where there is little heat transfer to or from the gas.Isothermal flow conditions assume that the flow through a pipe is maintained at a uniform temperature and these conditions are commonly assumed Science Publications AJAS when analyzing the flow of a gas in a non-insulated pipeline.Most natural gas pipelines are considered to be isothermal.Mekebel and Loraud (1983) experimentally investigated an unsteady compressible flow in a long gas transmission pipeline They concluded that the inclusion of heat transfer was necessary in the theoretical analysis, which was contrary to common assumptions that were made for isothermal or adiabatic flow systems.
The numerical simulation of unsteady internal compressible flows has been the goal of many researchers over the years and several algorithms have already been presented for developing such simulations.
The numerical methods that were previously used to solve systems of linear partial differential equations include the method of characteristics and a variety of explicit and implicit finite difference schemes (Osiadacz, 1996;Taylor et al., 1962;Wylie et al., 1971).Neglecting the inertia term in the momentum equation will result in a less accurate simulation result.
The dynamic behavior of the gas distribution pipeline network was simulated by MATLAB-Simulink.To simplify the computational process, some terms in the moment equation were assumed to be negligible through an estimation of their magnitudes, while maintaining the inclination term.The method of characteristics and the Crank-Nicolson method was applied to solve the model (Herran-Gonzalez et al., 2009).Mary et al. (2000) proposed an accurate secondorder algorithm for the simulation of unsteady, viscous and stratified compressible flows.The advantage of their method is its capability to address a broad range of subsonic Mach numbers, including nearly incompressible flows, with a single model that is based on the fully compressible Navier-Stokes equations.Gato and Henriques (2005) numerically modeled the dynamic behavior of high pressure natural gas flows in pipelines for a one dimensional compressible flow.The occurrence of pressure oscillations in natural gas pipelines as a result of the compression wave that originated from the rapid closure of downstream shutoff valves was investigated.
However, the advanced techniques that have been developed more recently in the area of computational fluid dynamics are now being investigated for solving the full conservation laws that govern the transport processes (Cadorin et al., 2010;Woldeyohannes and Majid, 2011;Zhou and Adewumi, 1995).In the work of Zhou and Adewumi (1995), a special type of hybrid Total Variation Diminishing (TVD) scheme with appropriate handling capability boundary conditions is developed, with the Roe scheme being used as the underlying algorithm.The formulation employed fixed point iterations rather than discrete time steps for developing the solution.Two recent reports in this area have been published by Cadorin et al. (2010) and Woldeyohannes and Majid (2011).The CFD commercial code ANSYS CFX 11.0 was used to evaluate the pressure drop through a pipe in a stream of high pressure gas that was characterized by a high Reynolds number in the study conducted by Cadorin et al. (2010).The flow was considered to be three-dimensional and the effect of turbulence was contemplated in equations for the first time based on the assumptions of isothermal flow and a steady mean limit.Woldeyohannes and Majid (2011) presented a simulation model for the analysis of a Transmission Pipeline Network System (TPNS) with detailed characteristics of the compressor stations using Visual C++ codes based on the Newton-Raphson solution technique to determine the pressure and the flow parameters of the network under various conditions.The model evaluated the energy consumption for various configurations to select the optimal TPNS and could assist in decisions regarding the design and operations of the pipeline.
Based on a review of the these papers, there is no comprehensive mathematical model using least restrictive hypotheses and a precise numerical technique that can be used to solve the model to be applicable in transient natural gas networks, which is the main purpose of this study and such as that developed in this study.Computational fluid dynamics are employed to solve the set of coupled partial differential equations that are developed from the conservation and turbulence equations without neglecting any terms.Finite volume and Multigrid methods have been used for discretizing the differential equations and achieving a convergence in the solution of algebraic equations.This scheme was found to be more stable, robust and accurate.Finally, one of the more practical applications of this simulation will be discussed.

Development of the Mathematical Model
The transient flows of gases in pipes are described by a two-dimensional approach.The basic equations that are used to describe the transient flows of gases in pipes are derived from the equation of motion (or momentum), the equation of continuity, the equation of energy and the ideal gas equation.We assume that the transmission line has a constant cross-sectional area and that the gas flow is highly turbulent.

Conservation of Mass: The Continuity Equation
Generally, the continuity equation is expressed in the form Equation 1 (Bird et al., 2006):

Newton's Second Law of Motion: The Momentum Equation
For gas flows in a pipe, the compressible Navier-Stokes equations can be written for the axis and the radial directions, respectively (Bird et al., 2006) Equation 2a and 2b: In these equations, u z and u r the axial and the radial velocity components, respectively.
The dimensionless viscous stress terms are defined for the tensor components as follows: In the above equations, µ eff is the effective viscosity (µ+µ t , where µ t is the turbulent viscosity, defined according to the turbulence model that was used) (Blazek, 2005).

Conservation of Energy
The basic form of the energy equation can be written as follows (Blazek, 2005;Ziaei-Rad and Nouri-Broujerdi, 2008) Equation 3: where, E is the total internal energy.The heat flux and the work performed by frictional forces in the energy equation are defined as follows: where, k eff is the effective conductivity (k+k t , where it is the turbulent thermal conductivity, defined according to the turbulence model that was used) (Blazek, 2005).

Equation of State
The equation of state for gases is related to the variables p, ρ and T. This equation is commonly expressed in the natural gas industry as follows Equation 4 (Kralik, 1988;Wylie and Streeter, 1978): The dynamic turbulent viscosity and the thermal conductivity are explained by the k-ε turbulence model (Blazek, 2005): , where C 0.09 The kinetic energies of the turbulence and the dissipation terms are computed by solving two transport Equation 5 and 6 (Launder and Spalding, 1974): where, C ε = 0.07 The right sides of Equations ( 5) and ( 6) contain the production and the destruction terms for ρk and ρε: where, c 1 = 0.129 and c 2 = 1.83.By definition, the following expressions for D and P in 2D can be used: The classical k-ε model is valid if the local Reynolds number is high.Therefore, it is not adequate to describe regions that are close to a solid wall.The idea is to use a two-layer approach by coupling K-ε the model to a oneequation model automatically, which enables us to compute the flow up to the wall with no empirical work required.However, this approach requires more computational resources as a finer mesh should be used.This method comprises introducing a local Reynolds number where, c 3 = 70 and c 4 = 0.41 (Crowl and Louvar, 2002).

Simulation Details
As mentioned in the introduction, the target of this study is to dynamically simulate the flow of natural gas in gas transmission lines to evaluate the sensitivity analysis of the effective parameters.
Therefore, the third Iranian Gas Trunkline (IGAT ΙΙΙ) that is used to transport natural gas from the south to the north of Iran for domestic consumption was examined.The composition of the gas in this trunk line and the main characteristics of this gas transmission system are shown in Fig. 1 and Table 1.
The pipelines between compressor stations six, seven and eight were selected for investigation.A grid refinement study was conducted to determine an adequate distribution and to ensure that the solution was not meshdependent.Finally, a mesh containing 873,760 and 1,103,583 cells were generated for the pipelines Ι and ΙΙ.

Numerical Techniques
CFD is a helpful tool for researchers and engineers to predict and analyze fluid flows and has been used since the early 1980s (Patankar, 1980).CFD is a numerical solution approach that is used to solve the coupled and highly nonlinear set of equations that characterize fluid flows and other engineering problems, such as chemical reactions.
As a first step in this numerical technique, the solution space or the domain is broken down into a large number of individual cells.Generating a good mesh for a particular flow situation is extremely important for CFD techniques.The equations are integrated over the volume of each cell in a discrete manner, which is also known as the Finite Volume Method.The introduction of a general variable φ to the conservative form of all fluid flow equations can be written as follows (Versteeg and Malalasekera, 2007) Equation 8: Equation ( 8) is the transport equation for the property φ.It clearly highlights the various transport processes: the rate of change term and the convective term on the left side and the diffusive term (Γ = diffusion coefficient) and the source term on the right side.To determine the common features among all these terms, we have to hide the terms that are not shared between the equations in the source terms.Applying the Gauss Divergence Theorem, Equation ( 8) can be rewritten as follows Equation 9: Science Publications

AJAS
In time-dependent problems, it is necessary to integrate using small increments of time ∆t (from t to t+∆t).The general forms of these integrals result in the following transport Equation 10: The integration of these equations employed an accurate second-order finite-volume method where the transport equations were discretized following a conservative control-volume approach, with the values of the dependent stored in the computational cell centers.Approximations of the diffusion and source terms were made using central differencing and a second-order upwind difference scheme was applied with respect to the convective and pressure fluxes.When an extended time period was used in this study (∆t = 1200s) an implicit scheme for unsteady state flows and the SIMPLER algorithm (Patankar, 1980;Patankar and Spalding, 1972) have been used to resolve the coupling between the velocity and the pressure.It is clear that the implicit method for any time period is unconditionally stable (Gerald and Wheatley, 2004) and for this reason, this method is recommended for solving transient fluid problems using CFD techniques.The values that are required for integration at the boundary of each cell are evaluated based on the neighboring cells and the boundary conditions that are supplied by the user.When the integrated equations at each cell are collected, they form a matrix of algebraic equations that can be solved using iterative numerical techniques.The Gauss-Seidel and the Incomplete LU decomposition (ILU) iterative methods (Axelsson, 1996) have been employed for scalar and coupled systems, respectively, because their application in computer programs is straightforward.However, in complicated problems, the convergence of the solutions is very slow.The Multigrid method (Chung, 2010) was utilized to overcome this problem and accelerate the convergence.
The total number of iterations to obtain the final converged results ranged from 650 to 700.The numerical computation was considered to be converged when the residual sum across all computational nodes at the nth iteration was less than or equal to 1.0e-05 for continuity and the velocities in the r and z directions.The convergence criterion for the energy equation was 1.0e-08, while the convergence criterion for k and ε was 1.0e-04).

Boundary Conditions
The boundary conditions for gas flow in a 2D pipeline are shown in Fig. 2.

Steady state condition
• Fluid boundary • Inlet: P, u, T are definite and ρ is defined according to the state equation • Outlet: The general condition of the fluid, which is commonly applied in the finite volume method, is as follows (Ferziger and Peric, 2002):

Initial Conditions for Unsteady State
The parameters P, u, T and ρ must be specified according to their initial values at t = 0 Therefore, the steady state condition data are also the initial conditions of the unsteady state.

RESULTS
When the CFD technique is used, a comparison between the simulation data that was produced from this amethod with the experimental (actual) data or from other successful projects is necessary to verify the accuracy of the simulated data.To illustrate the practicality of using this numerical simulation method on a gas transmission network, we compare the computational results with the real field data measured on-site from the Iran National Gas Company. Figure 3 shows the change in the gas flow rate in the pipeline over time.Figure 4 demonstrates the outlet pressure change as a function of time as the flow rate changes for both the simulated data from CFD and the information that was obtained from the National Gas Company According to the above figures, the field data and the predicted results show a satisfactory correlation, which implies that the CFD method is significantly accurate in predicting and analyzing the natural gas flow parameters in transmission pipelines, with an error of only 4 to 4.5%.

DISCUSSION
Hereon, an important case study related to the main goal of this study will be studied.During peak demand times in winter, the pipelines are operated at their maximum allowable flow rates.Therefore, compressor stations must operate with their maximum compression ratio to achieve desirable pressures at the end of the line.During these critical days, if one of the compressor stations fails to work for any reason, what changes have to be introduced into the pipeline system?What can be done to solve this problem?These questions will be discussed in the following paragraphs.
It is assumed that natural gas is flowing in the pipeline at the maximum operating flow rate (90 MMScMD).Compressor stations 6 and 7 compress the gas to a maximum operating condition of 72.74 bar or 1055 psi in the input line.The input gas temperature is 50°C.The natural gas pressure gradient along the length of the pipe under these conditions is shown in Fig. 5.The gas pressure is reduced to 57.13 bar at station 7 and it is compressed to 72.72 bar.Eventually, the pressure at the end of the line (station 8) is reduced to 47.81 bar, which is extremely close to the critical pressure (or the minimum allowable pressure) of 47.07 bar.
If compressor station 7 suddenly fails to work under these conditions, the outlet pressure changes in line ΙΙ as a function of time as illustrated in Fig. 6.
As shown, the compressor failed at 2 o'clock.A sudden pulse was observed simultaneously at the end of the line and the output pressure was reduced slowly to a steady value of 22 bar.
In gas transmission lines, when the pressure at the end of the line is dropped to less than 47 bar, the gas flow is cut from that compressor station onward.
One of the general solutions for this issue is to reduce the gas flow rate in the line.The gas pressure drop would decrease and the outlet gas pressure would therefore increase because of the reduction in the flow rate.

AJAS
The changes in the outlet pressure with flow rate at steady state can be seen in Fig. 7.If the flow rate is decreased to a value of 65 MMScMD, the outlet pressure is increased to 48.63 bars, regaining a proper pressure in the pipeline.
According to Fig. 8 when the outlet pressure was reduced to a steady value of 25.08 bar, there was a reduction in the flow rate after 22 h.The reduction in the flow reduction was conducted stepwise such that the outlet pressure increased to a steady value of 48.63 bar in h 36.
In this study, the natural gas flow rate in the pipeline was cut for approximately 25 h between hours 3 and 28.Using natural gas storage facilities would be very helpful during these situations to alleviate the gas production shortages.

CONCLUSION
In this study, the natural gas flow was modeled based on transient conditions in transmission lines.Mass transfer, momentum and energy equations were used for modeling the flow.Turbulence flow effects were applied in equations based on the two dimensional geometry of the system.For solving the set of six PDE coupled equations, Computational Fluid Dynamics (CFD) and the finite volume method were used to discretize the equations.The results that were obtained from this method were compared with the experimental data that was produced by National Gas Company.The errors between these comparisons were approximately 4 to 4.5 percent, which showed that the data generated from this model was accurate.
The results provide interesting insights into the transient behavior in transmission pipelines under operational scenarios and also indicate that a decrease in the flow rate leads to an increase in the outlet pressure at the pipe.In other words, in order to decrease the pressure loss, one should reduce the flow rate of the natural gas.
point that is closest to the wall and y is the distance between the current point and this point.To compute y + at values that are less than 200, is used according the following transport Equation 7(Apsley, 2007): iss = ρk 3/2 /I ε and the dynamic viscosity of turbulence is based on the relationship t c kI are two length scales that include the damping effects in the near wall regions and are defined as follows:

Fig. 1 .
Fig. 1.The structure of the gas transportation system

AJASFig. 2 .
Fig. 2. The flow domain and the boundary conditions

Fig. 3 .
Fig. 3.The changes in the pipeline flow rate with time out aiming considered mass flow rate.where, n is the normal outward vector of the outlet surface.• Solid boundary • No slip condition u = u w = 0 • Constant temperature on the wall T = T w • Symmetric boundary condition ∂φ/∂n = 0 Fig. 4. A comparison of the pressure outlet changes between the field data and the predicted results

Fig. 7 .
Fig. 5.The gas pressure gradient along the length of the pipeline during periods of peak demand

Table 1 .
Gas composition and main characteristics of IGAT-ΙΙΙ Compressibility factor was calculated for every discretization section of the pipeline using SGERG 88 (ISO, 2006) equation