Regular Article
Adaptive implicit finite difference method for natural gas pipeline transient flow
^{1}
Beijing Key Laboratory of Pipeline Critical Technology and Equipment for Deepwater Oil and Gas Development, Beijing Institute of Petrochemical Technology,
102617
Beijing, PR China
^{2}
National Engineering Laboratory for Pipeline Safety, Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum,
102249
Beijing, PR China
^{*} Corresponding author: email: yubobox@vip.163.com
Received:
3
June
2017
Accepted:
20
March
2018
The implicit finite difference method is one of the most widely applied methods for transient natural gas simulation. However, this implicit method is associated with high computational cost. To improve the simulation efficiency of implicit finite difference method, an adaptive strategy is introduced into the simulation process. The proposed adaptive strategy consists of the adaptive time step strategy and the adaptive spatial grid strategy. And these two parts are implemented based on the local error technique and the multilevel grid technique respectively. The results illustrate that the proposed adaptive method can automatically and independently adjust the time step and the spatial grid system according to the gas flow state in the simulation process, and demonstrates a significant advantage in terms of computational accuracy and efficiency compared with the nonadaptive method.
© P. Wang et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The simulation of natural gas transient flow plays a significant role in risk assessment and safety management of natural gas pipeline network (Karimpour et al., 2014). It has been widely studied since the 1960s, and different numerical methods have been proposed for solving the governing equations of transient gas flow in a pipeline. These methods can be summarized as follows: the characteristics method (Gutiérrez et al., 2002), implicit finite difference method (Helgaker et al., 2014), implicit finite volume method (Liang et al., 2013; Wang et al., 2011; Zhang, 2016), finite element method (Ebrahimzadeh et al., 2012), equivalent circuit method (Wang et al., 2014), state space model method (Alamian et al., 2012), reducedorder method (BehbahaniNejad and Shekari, 2010). Among these methods, the implicit finite difference method is one of the most widely applied methods for transient natural gas simulation. Many industrial codes are developed by using finite difference for transient flow (Evje and Flåtten, 2005; Coquel et al., 2006; Flåtten and Munkejord, 2006; Andrianov et al., 2007). Especially, the famous pipeline simulation software, Stoner Pipeline Simulator software (SPS) and Realpipe, are both developed by the finite difference method (Advantica, 2007; Zheng et al., 2012). The main reason is that the time step of this method is not restricted by the spatial step (Wylie et al., 1971; Kiuchi, 1994), which is quite useful for simulating longterm transient natural gas pipelines. However, in the implicit finite difference method, one system of largescale nonlinear algebraic equations must be solved at each time level in the transient simulation, which impedes the computation efficiency to a certain degree (Wylie et al., 1971). Therefore, a series of studies have been devoted to improving the efficiency of the implicit finite difference method.
The governing equations of transient gas flow are nonlinear hyperbolic partial differential equations, which leads to nonlinear algebraic equations. In the early stage of natural gas pipeline simulation studies, the convective inertia term (nonlinear term) of the governing equations was ignored directly to improve the efficiency of simulation (Wylie et al., 1971; Kiuchi, 1994). However, this process could reduce the accuracy of simulation (Abbaspour and Chapman, 2008). Luskin (1979) linearized the nonlinear governing equations about the previous time step based on the Taylor expansion. Zheng et al. (2012) found that the simulation efficiency is improved by more than 5 times via the linearization process, while the calculation accuracy is almost unaffected. To further reduce the nonlinearity of the hydraulic and thermodynamic system for the governing equtions, the decoupled solution strategy was proposed (Barley, 2012; Helgaker and Ytrehus, 2012). By the decoupled solution strategy, the hydraulic and thermodynamic equations can be solved alternatively, and the simulation efficiency can be increased by 20%. Additionally, Wang et al. (2015) found that the form of hydraulic equations that takes the density and velocity as the solving variables is the most efficient and easiest form when solving hydraulic systems. And the simulation efficiency can be further improved by 1.5 times. After the above pretreatment of the governing equations, the nonlinear iteration process is no longer required, while the largescale linear discretized equations still need to be solved. More seriously, the coefficient matrix of discretized equations is a large sparse irregular matrix that is neither diagonally dominant nor symmetrical arranged (Wylie et al., 1971). The sparse matrix technique is recommended to efficiently solve these discretized equations (Wylie et al., 1971). However, this technique is implemented with such difficultly that it is only applied in cases that have heavy computational burdens, such as the simulation of largescale natural gas pipeline network (Zheng et al., 2012). In recent years, Madoliat et al. (2016) proposed an approach based on intelligent algorithms, such as Particle Swarm Optimization (PSO). It converges over 200 times faster than the traditional elimination algorithms. Wang et al. (2018) and Yu et al. (2017) proposed a fast simulation method based on the divide and conquer concept. Its simulation speed is 1.5 times higher than that of Stoner Pipeline Simulator (SPS) software.
On the basis of the above studies, the efficiency of the implicit finite difference method for the natural gas flow simulation has been improved considerably. However, a common feature of these studies is the use of fixed time step and spatial grid. The situation may occur that a small time step and a dense spatial grid are excessively used to simulate the slow transient flow, while it could otherwise be accurately simulated by a large time step and a sparse spatial grid. Thus, computing resources are unnecessarily wasted by using the fixed time step and spatial grid system. Therefore, a method that can adaptively choose the time step and spatial grid according to the gas flow state and allocate computing resources on demand would further improve the simulation efficiency. The adaptive simulation method is simply one type of the above methods (Tao, 2000) and has been widely used in the simulation of various engineering problems, for example, reservoir simulation (Jackson et al., 2015), multiphase flow simulation (Pivello et al., 2014), and flooding simulation (Ruponen, 2014). However, it is seldom adopted in natural gas pipeline simulations, and few studies have been performed on the method. Only two studies related to adaptive transient gas flow simulation have been retrieved, Tentis et al. (2003) and Liang et al. (2011). The numerical results of the adaptive coarse grid are in good agreement with those of the fixed finer uniform grid, while the fixed coarse uniform grid gave poor results. However, these two methods were both developed on the method of lines, which is a kind of explicit method. Additionally, the time step of the explicit method is restricted by the spatial step, and such a procedure would not occur in the widely used implicit method. Therefore, no adaptive simulation is built on the implicit method of natural gas pipeline.
To further improve the efficiency of the implicit finite difference method, this paper focuses on the development of an implicit adaptive finite difference method. The adaptive implicit finite difference method adopts two successful methods. One is the local error technique which is the basis theoretical for adaptive time step; the other is the multilevel grid technique which is the basis theoretical for adaptive spatial grid. Then the detailed implementation of the adaptive time step strategy and the adaptive spatial grid strategy are designed respectively based on flow characteristics of natural gas in pipeline. Through this study, the time step and space grid can be adjusted intelligently at the same time in the natural gas pipeline simulation, and the computing resources can be allocated on demand. Then the large complex natural gas pipeline network can be efficient simulated by using the adaptive implicit finite difference method.
The layout of this paper is as follows: first, the main process of natural gas pipeline simulation using the implicit finite difference method is introduced. Then, an adaptive strategy for the implicit finite difference method simulation process is presented. Last, numerical experiments are used to evaluate the performance of the proposed implicit adaptive finite difference method.
2. Implicit finite difference method
The implicit finite difference method is briefly introduced in this section. Major elements of this method are the governing equations and discretization. More detailed process of this method can be found in reference (Zheng et al., 2012; Wang et al., 2015, 2018).
2.1. Governing equations
The governing equations of a natural gas pipeline transient flow consist of the continuity equation, momentum equation and energy equation. The continuity equation and momentum equation are also called hydraulic equations, and the energy equation is also called the thermodynamic equation. These equations can be written in a general form (Sanaye and Mahmoudimehr, 2012; Zheng et al., 2012; Duan et al., 2013; Wang et al., 2015, 2018), as shown in Equation (1), and the parameters of the general form are given in Table 1. (1) where c_{v} is the specific heat capacity at constant volume, d is the pipe internal diameter, g is gravitational acceleration, m is mass flowrate, p is pressure, t is time, w is velocity, x is the spatial coordinate, A is the crosssectional area, D is the pipe outside diameter, K is the total heat transfer coefficient, T is temperature, T_{g} is ambient temperature, λ is friction, ρ is density, and θ is the inclination angle of the pipe.
Equation (1) is a nonlinear hyperbolic equation, and it can be linearized about the previous time step based on the Taylor expansion as below (Wang et al., 2015, 2018; Zheng et al., 2012): (2) where the (i,j) elements of matrices G and S, respectively, are and , u_{i} is the ith corresponding component of U, and n equals 2 in the hydrodynamic equations and equals 1 in the thermodynamic equation. ,,, and are calculated using the variables of the previous time step.
The boundary conditions of the natural gas pipeline simulation are always given as the pressure or flow rate at the supply and the demand, as well as the temperature at the supply, as Equations (3)–(5). (3) (4) (5)
Parameters in the general form.
2.2. Discretization
In this paper, the decoupled solution strategy (Barley, 2012; Helgaker and Ytrehus, 2012) is adopted to solve the hydraulic and thermodynamic system. So the hydraulic equations and thermodynamic equation are discretized individually. The pipeline is divided into N sections, and thus, there are N + 1 points. The ith section is the section between the ith point and the (i + 1) th point.
For the hydraulic equations, at the th point, that is the middle of the ith section, the time derivative term and convection term are respectively discretized by the forward difference scheme and central difference scheme (Kiuchi, 1994; Abbaspour and Chapman, 2008; Wang et al., 2015, 2018). The discretization of the hydraulic equation can be written as follows: (6)
For the thermodynamic equation, at the ith point, the time derivative term and convection term equation are respectively discretized by the forward difference scheme and upwind scheme (Keenan, 1996; Barley, 2012). The discretization of the thermodynamic equation can be written as follows: (7) where
After introducing the boundary conditions, the discretized equations can be solved by the decoupled solution strategy (Barley, 2012; Helgaker and Ytrehus, 2012). By the decoupled solution strategy, the discretized hydraulic equations are firstly solved base on the interpolated temperature variables, and then the discretized thermodynamic equation is solved base on the solved hydraulic variables.
3. Adaptive strategy
For the implicit finite difference method introduced above, the time step and spatial grid are generally fixed in the simulation procedures. To adaptively choose the time step and spatial grid, an adaptive strategy for the implicit finite difference method is proposed in this section. The proposed adaptive strategy consists of the adaptive time step strategy and the adaptive spatial grid strategy.
3.1. Adaptive time step strategy based on the local error technique
3.1.1. General adaptive time step strategy
The time step is automatically determined according to the local error technique (Shampine, 2005). The adaptive time step strategy involves two main steps. First, the local error of the nth time level is estimated, which is introduced by the time step Δt^{n}. Second, the time step Δt^{n+1} of the (n + 1)th time level is calculated according to the local error. The detailed implementation procedures of the adaptive time step strategy are as follows: Conduct the pipeline simulation at the nth time level using the known time step Δt^{n} through Equations (6) and (7);Estimate the local error using the numerical solution of the (n−2)th, (n−1)th and nth time levels. (8) where is the local error.
Compare the value of local error to that of tolerable error TOL_{t}:
If , the error in the nth time level is tolerable, and the next time step Δt^{n+1} is calculated by the controller, that is Equation (9). Then go to step (1) to conduct the pipeline simulation at the (n + 1)th time level; (9) where β_{1}, β_{2} and β_{3} are the parameters of the controller. In this paper, H211 controllers (Shampine, 2005) are adopted, β_{1} = 0.25, β_{2} = 0.25, and β_{3} = − 0.25.
If , the error in time nth level is intolerable. Then, the time step is reduced to , and return to step (1) to conduct the pipeline simulation at the nth time level.
3.1.2. Improved process
In the adaptive time step strategy, the local error estimation plays a significant role. Traditionally, the local error is estimated by Equation (8). However, due to the characteristics of natural gas pipeline simulation are complicated, if the twoorder truncation error shown in Equation (8) is exclusively adopted to control the time step, the sensibility of the simulation to the change of the flow state in the pipeline would be low, and the corresponding simulation accuracy would be low as well.
Therefore, the error estimation process is further improved in this paper. The firstorder derivative has a higher sensitivity than the secondorder derivative (Larsen and Hansen, 2014). Thus, a firstorder truncation error is introduced when estimating the local error, and estimation of the local error can be written as Equation (10), (10)
3.1.3. Notes
In addition, there are two notes about the bound of the time step and the tolerable error: To avoid excessive small or large time steps, the upper and lower bounds of the time step can be set up in practical application, that is, Δt ∈ [Δt_{min}, Δt_{max}].The tolerable error of the adaptive time step strategy is generally set up based on one's experience and has a great relation with the scale of the pipeline. Because the adaptive method is seldom used in natural gas pipeline simulations, there is no rule for the setting of the tolerable error. Therefore, referring to the adaptive simulation of a reservoir (Jackson et al., 2015), multiphase flow (Pivello et al., 2014) and flooding (Ruponen, 2014), the tolerable error of the adaptive time method of the natural gas pipeline simulation in this paper is set as , and .
3.2. Adaptive spatial grid strategy based on the multilevel grid technique
3.2.1. General adaptive spatial grid strategy
The spatial grid is automatically determined by the multilevel grid technique (Vasilyev and Bowman, 2000). The adaptive spatial grid strategy involves two parts: one is the initialization of the multilevel grid before the simulation process, and the other is the spatial grid adaptive control in the simulation process.
3.2.1.1. Initialization of the multilevel grid
The basis of the adaptive spatial grid is the multilevel grid system, as shown in Figure 1, denoted as the V set, and the spatial grid system (Vasilyev and Bowman, 2000) can be described by the following equation: (11) where is the ith point in the jth grid level, and j = 0 indicates the coarsest spatial grid level, while j = J represents the densest grid level. N^{j} + 1 stands for the number of grids at the jth grid level, and Δx^{0} represents the size of the spatial grid at the sparsest level. It is obvious that .
Fig. 1 Sketch map of the multilevel grid structure. 
3.2.1.2. The adaptive spatial grid strategy at the (n + 1)th time level
The adaptive spatial grid control at the (n + 1)th time level involves two main steps. First, the error of the discretized points at the nth time level are estimated, which are the deviation between the values obtained by numerical calculation and the values obtained by the interpolation method. Second, a new spatial grid system is yielded for the (n + 1)th time level based on the deviation. It is worth noting that the point of the sparsest grid level should always be used. To facilitate the illustration, the V1 set and V2 set are defined as the collection of points used and the collection of points not used in the process of the solution of the nth time level, respectively. It is quite evident that V1 ∩ V2 = Φ (empty set) and V1 ∪ V2 = V. The procedures of the adaptive spatial grid strategy at the (n + 1)th time level are implemented as follows: Conduct the simulation on the nth time level using the known spatial grid system, which is the V1 set, through Equations (6) and (7).Interpolate the values of all points in the V2 set.
After conducting the simulation on the nth time level, only the values of natural gas flow parameters (pressure, flow rate, temperature, etc.) of points in the V1 set are obtained. However, the values of all points in the V set are the premise to estimate the error of the points belong to V1. Thus, the values of all points in the V2 set should be interpolated from j = 1 to j = J through Equation (12), (12) where u stands for p, m, T, i = 1, 2, ⋯ N^{j}, and .
Estimate the error of all points in the V1 set.
The values of all points in the V1 set are interpolated from j = 1 to j = J through the following equation: (13) where u stands for p, m, T, i = 1, 2, ⋯ N^{j}, and .
Then, the error of all points in the V1 set are estimated through Equation (14). (14) where ϵ_{x} is the estimation error.
Compare the error to the tolerable error TOL_{x}: If , the point is not used on the (n + 1)th time level and is moved to the V2 set from the V1 set.If , the point is conserved in the V1 set, and the point near this point at the jth and (j + 1)th grid level should be added to the V1 set, denoted as , , and .
3.2.2. Improved process
According to the characteristics of natural gas pipeline simulation, the interpolation of the pressure and the adjustment process of the spatial grid system are improved here.
In view of this fact, the squared value of the pressure and the mileage show a linear relationship in natural gas steady flow (Li and Yao, 2009). In Equations (12) and (13), however, interpolation is performed with the pressure value rather than the squared value of pressure. This interpolation method does not match the characteristics of the pressure distribution along the pipeline, and the simulation accuracy is low, especially on the coarse grid system. In this paper, we present an improved interpolation of pressure as below (see Eq. 15). (15)
In the general adaptive spatial grid strategy, the spatial grid system needs to be renewed on every time step, either by being refined or coarsened. This process not only disturbs the stability of the simulation but also introduces extra calculation workload, especially for the smooth period when the flow state in the pipeline changes slightly. Therefore, a transition zone is set to ensure the stability of the simulation process. The improved procedure is as follows.
Compare the error to the tolerable error TOL_{x}: If , the point is not used on the (n +
1)th time level and is moved to the V2 set from the V1 set. If , the point remains in the V1 set.If , the point remains in the V1 set, and the neighboring points of at the jth and (j + 1)th grid level should be added to the V1 set, which denoted as , , and .
Here, α_{x} is an adjustable parameter whose value ranges between 0 and 1.
3.2.3. Notes
In addition, there are two notes about the spatial grid system and the tolerable error. In this paper, the spatial grid system is separated into two independent parts: hydraulic and thermodynamic spatial grid systems. The hydraulic spatial grid system is controlled by the error of the pressure and flow variables, and the thermodynamic spatial grid system is controlled by the error of the temperature variables.Similar to the tolerable error of the adaptive time step strategy, there is no rule for the setting of tolerable error. Referring to the adaptive simulation of a reservoir (Jackson et al., 2015), multiphase flow (Pivello et al., 2014) and flooding (Ruponen, 2014), the tolerable error for the adaptive spatial grid method of the natural gas pipeline simulation is set as , and .
4. Results and discussion
In this section, first, the effects of the above improved process of the adaptive method are validated, and then, the computational accuracy and efficiency of the adaptive method are analyzed.
The properties of the following cases are set as follows. The standard pressure and standard temperature are, respectively, 101.325 kPa and 20 °C. The BWRS equations (Benedict et al., 1940) and Colebrook formula (Colebrook, 1939) are employed as the state equation and resistance equation, respectively. The components of the studied natural gas are listed in Table 2.
The following transient flow cases are simulated by the nonadaptive method, adaptive method and SPS. In the nonadaptive method, several different fixed time steps (6 s ≤ Δt ≤ 600 s) and spatial mesh sizes (0.1 km ≤ Δx ≤ 2 km) are adopted. In the adaptive method, the maximum and minimum time step of the adaptive method are, respectively, 3600 s and 60 s, and the multilevel grid is the 7 spatial level, and the space step of the 0th level Δx^{0} is 6 km.
Main components of the studied natural gas.
4.1. Result of the improved process of the adaptive method
Effects of the improved process regarding the local error estimation, the pressure interpolation and the adjustment process of the grid system are validated here. The length, diameter, thickness and roughness of pipeline are 200 km, 700 mm, 8 mm, and 0.02286 mm, respectively. At t = 0 h, the pipeline is in a shutdown state with the pressure of 10 MPa. After t = 0 h, the inlet pressure of the pipeline is kept at 10 MPa, and the outlet is opened to hold a flow rate of 1.0 × 10^{6} Nm^{3}/h. At t = 144 h, the outlet valve opening is reduced abruptly, and the outlet flow rate becomes 0.8 × 10^{6} Nm^{3}/h. In the entire simulation period, the temperature of the pipeline is maintained at 15 °C. The comparisons of the three kinds of improved processes are shown in Figures 2–5, where SPS refers to the numerical results obtained by the software SPS.
Figure 2 illustrates the comparison between the pipeline outlet pressures obtained by the adaptive strategy where the local error is estimated by Equations (8) and (10). It can be easily found that the description of the pipeline flow state with Equation (8) lags behind the SPS solution. Especially when the step change at t = 144 h, there is a lag of 1 h, and Equation (8) cannot accurately characterize the change process of the flow rate. Meanwhile, the description of the pipeline flow state with Equation (10) agrees well with that from SPS. The comparisons demonstrate that the improved process of the local error estimation is highly effective.
Figure 3 shows the distribution of the relative error of pressure at t = 10 h and t = 140 h, obtained by the adaptive strategy where the pressure is interpolated using Equations (13) and (15). It should be noted that the relative error is obtained by comparing with the solution of SPS. It can be seen that the maximum relative error of pressure obtained by Equation (13) is 3.5% and 4.5% at t = 10 h and t = 140 h, respectively, while the maximum relative error of pressure obtained with Equation (15) decreases to less than 1.5%. Therefore, the improved process of the pressure interpolation has higher accuracy than the traditional process.
Figure 4 displays the spatial grid number versus time before and after a transition zone is set in the procedure of error judgment. Figure 5 shows the inlet flow obtained by the adaptive strategy before and after a transition zone is set in the procedure of error judgment. From Figure 4, we can observe that the fluctuation of the spatial grid number is weaker than before when the transition zone is set. Additionally, in the smooth period, the spatial grid number remains totally constant after setting a transition zone, which indicates better stability of the grid system. The spatial grid number decreases too, which indicates less computation. Moreover, it can be found from Figure 5 that the transition zone has little effect on the inlet flow rate. The above results indicate that the ‘transition zone’ setting in the error judgment is able to strengthen the stability of the simulation while the accuracy is not lowered. Therefore, the improved process of refreshing spatial grid system has better stability than the traditional process.
Fig. 2 Outlet pressure versus time. 
Fig. 3 Absolute error distribution of pressure along the pipeline. 
Fig. 4 Spatial grid number versus time. 
Fig. 5 Inlet flow rate versus time. 
4.2. The computational accuracy and efficiency of the proposed adaptive method
In this section, the computational accuracy and efficiency of the adaptive method are tested by the horizontal pipe case, whose length, diameter, wall thickness and inner wall roughness are 24 km, 323 mm, 8 mm, and 0.02286 mm. The ground temperature is assumed to be constant with a value of 15 °C, and the total heat transfer coefficient K = 0.5 W/(m^{2 }· K) is adopted.
The initial conditions of the flow are as follows: the gas in the pipe is static, and it has a constant temperature of 15 °C and a constant pressure of 3 MPa along the entire pipe. The boundary conditions are as shown in Figure 6. It is clearly seen that the outlet flow rate suddenly jumps to 1.0 × 10^{5} Nm^{3}/h at the beginning, 0.5 × 10^{5} Nm^{3}/h at 24 h and 0.1 × 10^{5} Nm^{3}/h at 48 h, while the inlet temperature suddenly changes to 30 °C at the beginning and 40 °C at 24 h. The pressure at the inlet remains 3 MPa during the entire 72 h.
Fig. 6 Boundary conditions. 
4.3. Computational accuracy of adaptive method
To analyze the computational accuracy of adaptive method, the inlet flow rate, outlet pressure, and outlet temperature obtained by the three methods are compared in Figures 7–9.
It is seen clearly in Figures 7–9 that the numerical results of the large time step and the uniform coarse grid gave poor results. For example, under the calculation condition of fixed time step Δt = 600 s and fixed spatial step Δx = 2 km, nonphysical solutions for the inlet flow rates are obtained, as shown in Figure 7, and the results of outlet pressure obviously lag behind the results obtained by SPS, as shown in Figure 8, and the outlet temperature at t = 3 h is more than 50 °C, which is clearly unphysical, as shown in Figure 9. However, the numerical results of the adaptive method are in good agreement with those of SPS. These analyses qualitatively illustrate that the adaptive method can effectively guarantee the accuracy of the numerical results and can well describe the transient phenomena in the natural gas pipeline.
To compare the computational accuracy more conveniently, the simulation process is separated into two parts. One is the severe changing period, including the time from 0–2 h, 24–26 h and 48–50 h. The other is the gentle changing period, including the time from 2–24 h, 24–48 h and 50–72 h. The relative error of the outlet pressure and the absolute error of the outlet temperature are obtained by the two methods in the severe periods as well as gentle periods and are given in the Tables 3–6. It is worth noting that the error is obtained by comparison with the solution calculated with SPS software.
As seen from Table 3, under the same spatial grid system, the average relative error of the outlet pressure in severe periods obtained by the adaptive time step strategy is smaller than that obtained by the nonadaptive method with a fixed time step Δt = 60 s and is larger than that obtained by the nonadaptive method with fixed time step Δt = 6 s. For example, under the calculation condition of fixed spatial step Δx = 0.2 km, the average relative error of the outlet pressure obtained by the adaptive time step strategy is 0.2779%, which is between the 0.0968% obtained under the calculation condition of fixed time step Δt = 6 s and the 0.4695% obtained under the calculation condition of Δt = 60 s. This conclusion can also be obtained by analyzing the average absolute error of the outlet temperature in Table 4. It can be seen from Table 5 that the average relative errors of the outlet pressure in gentle periods obtained by the adaptive spatial grid are 0.2419%–0.2985%, which is almost the same as the 0.2105%–0.2175% obtained by the nonadaptive method with fixed spatial step Δx = 0.2 km. Meanwhile, the average absolute error of the outlet temperature in gentle periods obtained by the adaptive spatial grid is almost the same as the error obtained by the nonadaptive method with fixed spatial step Δx = 0.2 km. Therefore, the accuracy of the result obtained by the adaptive method is almost the same as those by the nonadaptive method with fixed time step Δt = 60 s and fixed spatial step Δx = 0.2 km. These analyses quantitatively demonstrate that the result obtained by the adaptive method has high accuracy.
Fig. 7 The flow rate at the inlet versus time. 
Fig. 8 The pressure at the outlet versus time. 
Fig. 9 The temperature at the outlet versus time. 
Average relative error of pressure at the outlet in the severe periods (%).
Average absolute error of temperature at the outlet in the severe periods (°C).
Average relative error of pressure at the outlet in the gentle periods (%).
Average absolute error of temperature at the outlet in the gentle periods (°C).
4.4. Computational efficiency of adaptive method
To analyze the computational efficiency of the adaptive method, the CPU time of the nonadaptive method and adaptive methods are compared in Table 7. The computations are carried out for a transient flow of 72 h, and the configuration of the computer used is as follows: Intel Core i7 860 CPU, 2.8 GHZ, 3 GB memory.
It is known from the previous analysis that, in this case, the accuracy of the result obtained by the adaptive method is almost the same as those by the nonadaptive method with fixed time step Δt = 60 s and fixed spatial step Δx = 0.2 km. However, the CPU times of the adaptive method and the nonadaptive method with Δt = 60 s and Δx = 0.2 km are, respectively, 0.2656 s and 1.5938 s, and the ratio is 0.2656/1.5938 ≈ 1/6. This means that the computational efficiency of adaptive method is 6 times higher than that of the nonadaptive method for this case. Therefore, the adaptive method can improve the efficiency of the gas pipeline simulation.
To determine the reasons for the improvement of efficiency, as well as to analyze the performance of the proposed adaptive method, Figures 10 and 11 , respectively, show the time step size and the spatial grid numbers used in the adaptive method.
It is seen clearly in Figure 10 that the time steps are 1–10 s at approximately 0 h, 24 h and 48 h, which show severe changing periods, while the time steps are the maximum 3600 s during 4–20 h, 28–44 h and 52–68 h, which are gentle periods. It is seen clearly in Figure 11 that the hydraulic spatial grid numbers, respectively, reach the peak of approximately 200, 60 and 70 in the severe changing periods, while the hydraulic spatial grid numbers reach, respectively, only approximately 22, 22 and 7 during the gentle periods. The thermodynamic spatial grid numbers, respectively, reach peaks of approximately 200, 170 and 150 in the severe periods, while the thermodynamic spatial grid numbers, respectively, reach only approximately 120, 38 and 70 during the gentle changing periods. The above results show that the adaptive method can automatically adopt a small time step and a dense grid system when the flow state changes severely and adopt a large time step and a coarse grid system when the flow state changes gently. Therefore, the adaptive method can improve the computational efficiency while retaining considerable accuracy.
In addition, there are a total of 489 iterations in Figure 10, that is, the average time step of the adaptive method is 72 × 3600/489 ≈ 530 s. The average grid numbers of the hydraulic system and thermodynamic system in Figure 11 are 45 and 127, respectively, which indicates that the mean spatial steps of the hydraulic system and thermodynamic system in the adaptive method are 24 km/45 ≈ 0.5 km and 24/127 ≈ 0.2 km. The accuracy of the result obtained by the adaptive method is almost identical to the nonadaptive method with Δt = 60 s and Δx = 0.2 km. That is, the average time and spatial step of the adaptive method are larger than those of the adaptive method. This is the main reason why the computational efficiency of the nonadaptive method is 6 times higher than that of the nonadaptive method when retaining the same computation accuracy.
Compared to the conventional method, the adaptive method proposed in this paper can adaptively determine the time step and spatial grid according to the gas flow state and allocate computing resources on demand, which would further improve the simulation efficiency.
CPU time (s).
Fig. 10 Time step versus time. 
Fig. 11 Spatial grid number versus time. 
5. Conclusion
In this paper, the adaptive implicit finite difference method for gas pipeline simulation is proposed, which consists of two totally independent parts: the adaptive time step strategy and the adaptive spatial grid strategy. Then, a restart process of natural gas pipeline is simulated to evaluate the advantages of the method proposed. The following conclusions can be drawn:

the proposed adaptive method has high accuracy. In the case of this paper, the accuracy of the result obtained by the adaptive method is almost the same with the nonadaptive method with fixed time step Δt = 60 s and fixed spatial step Δx = 0.2 km, as well as SPS software;

the proposed adaptive method can further improve the simulation efficiency. The main reason is that this method can adjust the time step and the spatial grid system according to the flow state. That is, the adaptive method can allocate the calculation resources on demand to save time consumption. As a result, to reach a certain accuracy, the computation speed of the adaptive method is 6 times faster than that of the nonadaptive method in the given case.
Acknowledgments
The study is supported by the Project of Construction of Innovative Teams and Teacher Career Development for Universities and Colleges Under Beijing Municipality (No. IDHT20170507), the Program of Great Wall Scholar (No. CIT&TCD20180313), and the Beijing Talents Project.
References
 Abbaspour M., Chapman K.S. (2008) Nonisothermal transient flow in natural gas pipeline, Int. J. Appl. Mech. 5, 3, 031018. [CrossRef] [Google Scholar]
 Advantica I. (2007) Stoner Pipeline Simulator (SPS) Help and Reference, 9.600 Bent Creek Blvd. [Google Scholar]
 Alamian R., BehbahaniNejad M., Ghanbarzadeh A. (2012) A state space model for transient flow simulation in natural gas pipelines, J. Nat. Gas Sci. Eng. 9, 51–59. [CrossRef] [Google Scholar]
 Andrianov N., Coquel F., Postel M., Tran Q.H. (2007) A relaxation multiresolution scheme for accelerating realistic two‐phase flows calculations in pipelines, Int. J. Numer. Meth. Fluids 54, 2, 207–236. [CrossRef] [Google Scholar]
 Barley J. (2012) Thermal decoupling: An investigation, PSIG Annual Meeting, Pipeline Simulation Interest Group, Santa Fe, New Mexico. [Google Scholar]
 BehbahaniNejad M., Shekari Y. (2010) The accuracy and efficiency of a reducedorder model for transient flow analysis in gas pipelines, J. Petrol. Sci. Eng. 73, 1, 13–19. [CrossRef] [Google Scholar]
 Benedict M., Webb G.B., Rubin L.C. (1940) An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures I. Methane, ethane, propane and n‐butane, J. Chem. Phys. 8, 4, 334–345. [CrossRef] [Google Scholar]
 Colebrook C.F. (1939) Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws, J. Inst. Civil Eng. 11, 4, 133–156. [CrossRef] [Google Scholar]
 Coquel F., Postel M., Poussineau N., Tran Q.H. (2006) Multiresolution technique and explicitimplicit scheme for multicomponent flows, J. Numer. Math. 14, 3, 187–216. [CrossRef] [Google Scholar]
 Duan J.M., Wang W., Zhang Y., Zheng L.J. (2013) Energy equation derivation of the oilgas flow in pipelines, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 68, 2, 341–353. [CrossRef] [Google Scholar]
 Ebrahimzadeh E., Shahrak M.N., Bazooyar B. (2012) Simulation of transient gas flow using the orthogonal collocation method, Chem. Eng. Res. Des. 90, 11, 1701–1710. [CrossRef] [Google Scholar]
 Evje S., Flåtten T. (2005) Weakly implicit numerical schemes for a twofluid model, SIAM J. Sci. Comput. 26, 5, 1449–1484. [CrossRef] [Google Scholar]
 Flåtten T., Munkejord S.T. (2006) The approximate Riemann solver of Roe applied to a driftflux twophase flow model, ESAIM: Math. Model. Numer. Anal. 40, 4, 735–764. [CrossRef] [Google Scholar]
 Gutiérrez J.A., Moreno P., Naredo J.L., Gutiérrez J.C. (2002) Fast transients analysis of nonuniform transmission lines through the method of characteristics, Int. J. Electr. Power & Energy Syst. 24, 9, 781–788. [CrossRef] [Google Scholar]
 Helgaker J.F., Müller B., Ytrehus T. (2014) Transient flow in natural gas pipelines using implicit finite difference schemes, J. Offshore Mech. Arct. Eng. 136, 3, 031701. [CrossRef] [Google Scholar]
 Helgaker J.F., Ytrehus T. (2012) Coupling between continuity/momentum and energy equation in 1D gas flow, Energy Procedia 26, 82–89. [CrossRef] [Google Scholar]
 Karimpour K., Zarghami R., Moosavian M.A., Bahmanyar H. (2014) New fuzzy model for risk assessment based on different types of consequences, Oil Gas Sci. Technol. 71, 17. [Google Scholar]
 Keenan P.T. (1996) Collation and upwinding for thermal flow in pipelines: the linearized case, Int. J. Numer. Meth. Fl. 22, 9, 835–849. [CrossRef] [Google Scholar]
 Kiuchi T. (1994) An implicit method for transient gas flows in pipe networks, Int. J. Heat Fluid Fl. 15, 5, 378–383. [Google Scholar]
 Jackson M., Percival J., Mostaghimi P., Tollit B., Pavlidis D., Pain C., Gomes J., Elsheikh A.H., Salinas P., Muggeridge A., Blunt M. (2015) Reservoir modeling for flow simulation by use of surfaces, adaptive unstructured meshes, and an overlappingcontrolvolume finiteelement method, SPE Reserv. Eval. Eng. 18, 2, 115–132. [CrossRef] [Google Scholar]
 Larsen P.M., Hansen N.E. (2014) Computer aided design in control and engineering systems: advanced tools for modern technology, Pergamon Press, New York. [Google Scholar]
 Li Y.X., Yao G.Z. (2009) Design and operation of gas pipeline, China University of Petroleum Press, Beijing (in Chinese). [Google Scholar]
 Liang Y.M., Zheng G., Liu H. (2011) Dynamic simulation of natural gas network based on adaptive step, Comput. Eng. Appl. 47, 7, 233–235 (in Chinese). [Google Scholar]
 Liang J.F., Guo K., Huangpu L.X., Li N., Wang G.P. (2013) Operation analysis of city gas pipelines by the finite volume method, Nat. Gas Ind. 33, 104–108 (in Chinese). [Google Scholar]
 Luskin M. (1979) An approximation procedure for nonsymmetric, nonlinear hyperbolic systems with integral boundary conditions, SIAM J. Numer. Anal. 16, 1, 145–164. [Google Scholar]
 Madoliat R., Khanmirza E., Moetamedzadeh H.R. (2016) Transient simulation of gas pipeline networks using intelligent methods, J. Nat. Gas Sci. Eng. 29, 517–529. [CrossRef] [Google Scholar]
 Pivello M.R., Villar M.M., Serfaty R., Romab A.M., SilveiraNetoa A. (2014) A fully adaptive front tracking method for the simulation of two phase flows, Int. J. Multiphas. Flow 58, 72–82. [CrossRef] [Google Scholar]
 Ruponen P. (2014) Adaptive time step in simulation of progressive flooding, Ocean Eng. 78, 35–44. [CrossRef] [Google Scholar]
 Sanaye S., Mahmoudimehr J. (2012) Technical assessment of isothermal and nonisothermal modelings of natural gas pipeline operational conditions, Oil Gas Sci. Technol. 67, 3, 435–449. [Google Scholar]
 Shampine L.F. (2005) Error estimation and control for ODEs, J. Sci. Comput. 25, 1, 3–16. [CrossRef] [Google Scholar]
 Tao W.Q. (2000) Advances in computational heat transfer, Science Press, Beijing (in Chinese). [Google Scholar]
 Tentis E., Margaris D., Papanikas D. (2003) Transient gas flow simulation using an adaptive method of lines, C.R. Mec. 331, 7, 481–487. [CrossRef] [Google Scholar]
 Vasilyev O.V., Bowman C. (2000) Secondgeneration wavelet collocation method for the solution of partial differential equations, J. Comput. Phys. 165, 2, 660–693. [CrossRef] [Google Scholar]
 Wang H., Liu X.J., Zhou W.G. (2011) Transient flow simulation of municipal gas pipelines and networks using semi implicit finite volume method, Procedia Eng. 12, 217–223. [CrossRef] [Google Scholar]
 Wang J.R., Wang T., Wang J.Z. (2014) Application of π equivalent circuit in mathematic modeling and simulation of gas pipeline, Appl. Mech. Mater. 496, 943–946. [CrossRef] [Google Scholar]
 Wang P., Yu B., Deng Y.J., Zhao Y. (2015) Comparison study on the accuracy and efficiency of the four forms of hydraulic equations of a natural gas pipeline based on linearized solution, J. Nat. Gas Sci. Eng. 22, 235–244. [CrossRef] [Google Scholar]
 Wang P., Yu B., Han D., Sun D., Xiang Y. (2018) Fast method for the hydraulic simulation of natural gas pipeline networks based on the divideandconquer approach. J. Nat. Gas Sci. Eng. 50, 55–63. [CrossRef] [Google Scholar]
 Wylie E.B., Stoner M.A., Streeter V.L. (1971) Network: System transient calculations by implicit method, Soc. Petrol. Eng. J. 11, 04, 356–362. [Google Scholar]
 Yu B., Wang P., Wang L.Y., Xiang Y. (2017) A simulation method for natural gas pipeline networks based on the divideandconquer concept, Oil Gas Storage Transp. 36, 1, 75–84 (in Chinese). [Google Scholar]
 Zhang L. (2016) Simulation of the transient flow in a natural gas compression system using a highorder upwind scheme considering the realgas behaviors, J. Nat. Gas Sci. Eng. 28, 479–490. [CrossRef] [Google Scholar]
 Zheng J.G., Chen G.Q., Song F., AiMu Y., Zhao J.L. (2012) Research on simulation model and solving technology of large scale gas pipe network, J. Syst. Simul. 14, 133 (in Chinese). [Google Scholar]
All Tables
All Figures
Fig. 1 Sketch map of the multilevel grid structure. 

In the text 
Fig. 2 Outlet pressure versus time. 

In the text 
Fig. 3 Absolute error distribution of pressure along the pipeline. 

In the text 
Fig. 4 Spatial grid number versus time. 

In the text 
Fig. 5 Inlet flow rate versus time. 

In the text 
Fig. 6 Boundary conditions. 

In the text 
Fig. 7 The flow rate at the inlet versus time. 

In the text 
Fig. 8 The pressure at the outlet versus time. 

In the text 
Fig. 9 The temperature at the outlet versus time. 

In the text 
Fig. 10 Time step versus time. 

In the text 
Fig. 11 Spatial grid number versus time. 

In the text 