Numerical modeling of transients in gas pipeline

A set of equations governing an isothermal compressible fluid flow is analytically and numerically analyzed. The obtained equations are written in characteristic form and resolved by a predictorcorrector lambda scheme for the interior mesh points. The method of characteristics (MOC) is used for the boundaries. Advantages of explicit form of these schemes and the flexibility of the MOC are used for an isothermal fast transient gas flow in short pipeline. The results, obtained for a simple practical application agree with those of other methods.


INTRODUCTION
Gas transients equations in pipelines can be linear or generally non linear.They also may be parabolic or hyperbolic of the first or second order.As a rule, simple models are an alternative which presents a reasonable compromise between the description accuracy and the cost of solution.These models are obtained by neglecting some terms in the basic set of equations, as a result of a quantitative estimation of the particular elements of the equations for given operating conditions of the pipeline.
Several relatively new numerical schemes were tested to integrate equations of conservation.These include Godunov and TVD schemes (Leveque and Yee, 1990).These second order schemes have the advantage that shock waves problems and other discontinuities can be treated with relatively good accuracy.
In their studies, many authors have considered fast gas transients employing numerical techniques as method of characteristics (MOC) (Kameswara and Eswaran, 1993;Greyvenstein, 2001) and finite differences (Greyvenstein, 2001;Gato and Henriques, 2005), with a relatively good agreement each other.In the last decades Behbahani-Nejad and Bagheri (2010a) have used transfer functions of a single pipeline in order to develop a mathematical Simulik library.Obtained results are satisfactory with the classical methods.
With a reduced order modeling approach, Behbahani-Nejad and Shekari (2010b) have compared their results with the conventional numerical techniques for a simple gas transient example.A good agreement was observed.By the use of time space least square spectral method with a technique based on hierarchical interpolations in space and time, on numerical examples, including fast gas transients, particularly in the case of severe conditions flowing, Dorao and Fernandino (2011) have successfully handled the problem of strong shock wave.
Simulating gas transients in pipes, Ebrahimzadeha et al. ( 2012) have used an orthogonal collocated method technique to solve the corresponding governing equations.
Its performances are verified and tested for two practical examples corresponding to isothermal and non isothermal cases.
In this work an, old idea, relative to a physically meaningful simples schema, have been developed by Moretti (1979), Zannetti and Colasurdo (1981) and Gabutti (1983).They are known as the lambda schemes.An important advantage of these schemes is related to the concept of non reflecting boundary condition :a characteristics form of the boundary conditions equations is applied in order to avoid the use of the improperly reflecting technique.An explicit form of them is adapted and applied to an isothermal fast transient gas flow in short pipeline.

THEORETICAL MODELING
If we consider an isothermal flow in pipeline with variable crosssectional area in which one dimensional continuity equation is: and momentum equation is given by : Writing equation of state for natural gas as: and taking into account the isothermal conditions, the acoustic wave speed becomes: In the absence of field data, steady state variable distributions constitute the initial conditions.These steady state initial conditions are obtained by the use of an appropriate analytical equation (Zhou and Adewumi (1996) for Z=1: . This equation, which is implicit in  , is well suited for iterative method to determine density or pressure distribution.Kessal et al. 83 In earlier work, considering above equation set, two applications of fast and slow fluid flows have been developed by using two explicit finite-difference schemes (Kessal, 2000).Then, in the same way, a one-dimensional lambda scheme is proposed to study the first case.

NUMERICAL SCHEME
In order to analyze the gas transients phenomena in short pipelines the characteristic method is used to convert the initial partial differential equation set (1) and (2) into ordinary differential equations (Lister, 1960).The physical interpretation is that the waves travel with the speed "a", given by the relation (4), propagating in this way the effect of the initial boundary conditions.Then the transformation of Equations ( 1) and ( 2) yields: These equations are associated respectively with the following characteristics directions: Using Equations ( 8) and ( 9), Equations ( 6) and ( 7) can be written as: Adding Equations ( 10) and ( 11) and simplifying yields to: Subtracting Equation (10) from Equation ( 11) and simplifying, yields to: Application of an explicit lambda scheme for the above equation set needs the following transformations.
Equations ( 12) and ( 13) are in so-called lambda form.Note that the spatial derivatives are marked with subscripts + andto indicate the characteristic directions along which these derivatives are approximated.Predicted values of * i V and * i P can be obtained by substitution of finite-difference approximations for the time derivatives into Equations ( 12) and ( 13): where: Then, the predictor-corrector scheme applied to our equation set yields to the following procedure (Gabutti, 1983).The spatial derivatives in Equations ( 14) and ( 15) are approximated as follow: Replacing the above finite-difference approximations in Equations ( 14) and (15) yields the predicted values of * i V and * i P .

Part 2
The predicted values of the time derivatives 12) and ( 13) are calculated by using the following finite difference approximations: x

Corrector part
By considering the following finite difference approximations and using V * and P * instead of V and P in Equations ( 12) and ( 13 x Finally, the values of P and V at the unknown time level are determined from the following equations: It can be noticed that above discretization is possible only at nodes 3, 4,…, N-1 as it is shown by part 2 of the predictor scheme.Then, a special treatment is needed at points near the boundaries.For this, a one sided finite-difference approximation can be used at nodes 2 and N-1. Note that in Gabutti (1983) paper, two points finite-difference approximations was used at the nodes adjacent to the boundaries if three points were not available in the desired directions.In this study computational time interval was selected so that the Courant stability condition was satisfied at all nodes of the mesh (Streeter and Wylie, 1969).If necessary, the time interval can be reduced in some cases.

INITIAL AND BOUNDARY CONDITIONS
Initial and boundary conditions to the previous Equations ( 1) and ( 2) must be specified in order to obtain the appreciable solution for this differential equation set.Initial conditions of these systems are required to resolve initial pressure and velocity as a function of the position x along the pipeline.In this study they are given by the relation (8) for the pressure distribution.Boundary conditions must also be specified to obtain a unique solution.They depend on the considered cases.It is proposed in this work to treat numerically the two boundary conditions by the characteristics method.Then, integrating Equations ( 6) and ( 7) along the negative and the positive characteristics lines Equations ( 8) and (9) (Figure 1) yields to the following finite-difference equations: Applying these equations to limit conditions yields (Figure 2): A computational procedure to obtain P or V is necessary with the

RESULTS AND DISCUSSION
An example concerning a gas transient in a relatively short pipeline with an impulse supply of gas mass flux at the inlet, has been simulated using the previous predictor-corrector scheme.This example is taken from Zhou and Adewumi (1996) in which solutions have been obtained, respectively, by using the method of characteristics (MOC), and, a first-order three-point explicit Godunov scheme and a source free secondorder five-point TVD scheme.A pipeline 91.44 m long, 0.609 m interior diameter and having initially a static pressure of 4136.8 kPa (with initial velocity V i =0) with a shut downstream extremity.At time zero, upstream inflow begins to increase linearly and reaches 196 Kg/s at 0.145 s, then decreases linearly to zero again at 0.29 s, and then remains constant.The downstream end is closed.
For the simulation of the above fast transient problem, the predictor-corrector lambda scheme adopts with the characteristics method, at the boundaries, the same t imposed by the stability Criteria (28).
For numerical simulation of this example, the previous predictor-corrector scheme adopts a grid size x=0.9144m, t=0.811x 10 -3 s and C.F.L=0.312.These values can be compared with those used in the Godunov and TVD 0,0 0,2 0,4 0,6 0,8  Schemes (19) respectively: x=0.9144m, t=0.09144x 10 -3 s and C.F.L=0.0348;x=0.9144m, t=0.9 x 10 -3 s and C.F.L=0.348.Some numerical results are illustrated by the following figures.It shown a time evolution rate of the gas (Figure 3) at midpoint of the gas line (x=0.5*L),where our results are compared with those obtained by Zhou and Adewumi (1996).Good agreement between them can be observed.The gas flow rate fronts are completely solved within the first 0.8 s.The behaviour of gas flow rate evolution at the midpoint of the pipeline is the result of reflected pressure impulse at the upstream end of the pipe.A comparison of the predicted pressure (by the present model) at the inlet of the pipeline (x/L=0) with the reported data as shown in   reported data is obtained.It can be noticed that duration of the pressure pulse of the peak is the same.Pressure wave is maintained and captured during the first 0.80 s.
The pressure wave fronts are reproduced during 2.4 s, without any loss of accuracy (Figure 5).Slight differences with an other methods can be due to some simplifications introduced by Zhou and Adewumi (1996), that is, the value of the friction losses coefficient.A comparison, as regards the pressure at the outlet point of the pipeline (x/L=1), between obtained numerical results and the reported data is shown in Figure 6.Good agreement can be noticed.At the outlet of the pipeline, the pressure wave fronts are completely resolved within the first 0.8 s.
In order to check the numerical method described herein, the pressure evolution has been calculated for time close to the end of the transient phenomenon.It is shown in Figure 7 that agreement is satisfactory.This indicates that the method described in the previous  sections is reliable.Also computed results show that the employed numerical simulation has not produced any undesirable effect.The evolution of the pressure, corresponding to successive times relative to the pipeline filling, can be observed in Figure 8. Depending on the speed of propagation of sound in gas, then we note that the wave front (of pressure) is reflected from the downstream end.In Figure 9, we can observe the correspondence with Figures 2 and 4, namely the time change in pressure at the inlet and outlet of the pipeline.The incident and reflected pressure waves are explicitly shown in this figure.The longitudinal evolution of the flow rate from the inlet of the pipe is shown in Figure 10.Nevertheless, the speed reflected by the downstream end takes negative values to the upstream end, thereby causing a depression waves which propagate to the upstream.
Finally, the propagation phenomenon of the pressure and the flow velocity disturbances, initiated at upstream point that the pipeline, is well reproduced by the two   forms of difference schemes, that is, a second order scheme for the interior points and a method of characteristics for the extremities.The gas transients are contained in a very comprehensive and convenient.

Figure 1 .
Figure 1.x, t Grid for method of characteristics.

Figure 3 .
Fig 3. Flow rate variation at the midpoint of the pipeline

Fig 4 .
Fig 4. Pressure response at the inlet of the pipe.

Figure 4 .
Figure 4. Pressure response at the inlet of the pipe.

Figure 5 .
Figure 5. Pressure response at the inlet of the pipeline.

Figure 5 .Figure 6 .
Figure 5. Pressure response at the inlet of the pipeline.

Figure 7 .
Figure 7. Pressure response at the inlet of the pipe.

Figure 9 .
Figure 9. 3D Gas transients between the inlet and outlet

Figure 9 .
Figure 9. 3D Gas transients between the inlet and outlet.