© P. Wang et al., published by IFP Energies nouvelles, 2018

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), reduced-order method (Behbahani-Nejad 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 long-term transient natural gas pipelines. However, in the implicit finite difference method, one system of large-scale 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 large-scale 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 large-scale 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. $$\frac{\partial U}{\partial t}+B\cdot \frac{\partial U}{\partial x}=F$$(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 cross-sectional 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): $$\frac{\partial U}{\partial t}+{\displaystyle \bar{B}}\cdot \frac{\partial U}{\partial x}+{\displaystyle \bar{G}}\cdot (U-{\displaystyle \bar{U}})={\displaystyle \bar{F}}+{\displaystyle \bar{S}}\cdot (U-{\displaystyle \bar{U}})$$(2) where the (i,j) elements of matrices G and S, respectively, are ${\left[G\right]}_{i,j}={\displaystyle \sum _{l=1}^n}{(\frac{\partial B}{\partial u_j})}_{i,l}\frac{\partial u_l}{\partial x}$ and ${[S]}_{i,j}=\frac{\partial F_i}{\partial u_j}$, u_i is the ith corresponding component of U, and n equals 2 in the hydrodynamic equations and equals 1 in the thermodynamic equation. ${\displaystyle \bar{B}}$,${\displaystyle \bar{G}}$,${\displaystyle \bar{F}}$,${\displaystyle \bar{S}}$ and ${\displaystyle \bar{U}}$ 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). $$\text{Give the pressure value }p=p(t)$$(3) $$\text{Give the flow rate value }m=m(t)$$(4) $$\text{Give the temperature value }T=T(t)$$(5)

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 $\left(i+\frac{1}{2}\right)$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: $$C{E}_i\cdot U_i^n+D{W}_i\cdot U_{i+1}^n=H_i$$(6) $${\displaystyle \bar{U}}=\frac{U_{i+1}^{n-1}+U_i^{n-1}}{2}\text{,}$$ $$C{E}_i=\frac{1}{2\mathrm{\Delta }{t}^n}I-\frac{1}{x_{i+1}-x_i}{\displaystyle \bar{B}}+\frac{1}{2}({\displaystyle \bar{G}}-{\displaystyle \bar{S}})\text{,}$$ $$D{W}_i=\frac{1}{2\mathrm{\Delta }{t}^n}I+\frac{1}{x_{i+1}-x_i}{\displaystyle \bar{B}}+\frac{1}{2}({\displaystyle \bar{G}}-{\displaystyle \bar{S}})\text{,}$$ $$H_i={\displaystyle \bar{F}}+({\displaystyle \bar{G}}-{\displaystyle \bar{S}}+\frac{1}{\mathrm{\Delta }{t}^n})\cdot {\displaystyle \bar{U}}\text{.}$$

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: $${\mathit{\text{UP}}}_i\cdot T_{i-1}^n+{\mathit{\text{CE}}}_i\cdot T_i^n+{\mathit{\text{DW}}}_i\cdot T_{i+1}^n=H_i\text{,}$$(7) where ${\mathit{\text{UP}}}_i=-\max \left(w_i^n,0\right)\frac{1}{x_{i+1}-x_i}{\mathit{\text{\hspace{1em}CE}}}_i=\left(\frac{1}{\mathrm{\Delta }t}-{\displaystyle \overline{S}}+\left|w_i^n\right|\frac{1}{x_{i+1}-x_i}\right)\text{,}$ $${\mathit{\text{DW}}}_i=-\max \left(-w_i^n,0\right)\frac{1}{x_{i+1}-x_i}{H}_i={\displaystyle \overline{F}}+\left(\frac{1}{\mathrm{\Delta }{t}^n}-{\displaystyle \overline{S}}\right)T_i^{n-1}\text{.}$$

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ⁿ. Second, the time step Δtⁿ⁺¹ 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ⁿ 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. $${ϵ}_t^n={\parallel \left(\frac{U^n-U^{n-1}}{\mathrm{\Delta }{t}^n}-\frac{U^{n-1}-U^{n-2}}{\mathrm{\Delta }{t}^{n-1}}\right)\frac{\mathrm{\Delta }{t}^n+\mathrm{\Delta }{t}^{n-1}}{2}\parallel }_2\text{,}$$(8) where ${ϵ}_t^n$ is the local error.

Compare the value of local error ${ϵ}_t^n$ to that of tolerable error TOL_t:

If ${ϵ}_t^n\le {\text{TOL}}_t$, the error in the nth time level is tolerable, and the next time step Δtⁿ⁺¹ 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; $$\mathrm{\Delta }{t}^{n+1}={\left(\frac{{\text{TOL}}_t}{{ϵ}_t^n}\right)}^{{\beta }_1}{\left(\frac{{\text{TOL}}_t}{{ϵ}_t^{n-1}}\right)}^{{\beta }_2}{\left(\frac{{ϵ}_t^n}{{ϵ}_t^{n-1}}\right)}^{{\beta }_3}\mathrm{\Delta }{t}^n\text{,}$$(9) where β₁, β₂ and β₃ are the parameters of the controller. In this paper, H211 controllers (Shampine, 2005) are adopted, β₁ = 0.25, β₂ = 0.25, and β₃ = − 0.25.

If ${ϵ}_t^n>{\text{TOL}}_t$, the error in time nth level is intolerable. Then, the time step is reduced to $\mathrm{\Delta }{t}^n=\frac{\mathrm{\Delta }{t}^n}{2}$, 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 two-order 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 first-order derivative has a higher sensitivity than the second-order derivative (Larsen and Hansen, 2014). Thus, a first-order truncation error is introduced when estimating the local error, and estimation of the local error can be written as Equation (10), $${ϵ}_t^n=\max \left\{{\parallel \left(\frac{U^n-U^{n-1}}{\mathrm{\Delta }{t}^n}-\frac{U^{n-1}-U^{n-2}}{\mathrm{\Delta }{t}^{n-1}}\right)\frac{\mathrm{\Delta }{t}^n+\mathrm{\Delta }{t}^{n-1}}{2}\parallel }_2{\parallel U^n-U^{n-1}\parallel }_2\right\}\text{.}$$(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 $\text{TO}{\text{L}}_t^p={\parallel p\parallel }_2\times {10}^{-3}$, $\text{TO}{\text{L}}_t^m={\parallel m\parallel }_2\times {10}^{-1}$ and $\text{TO}{\text{L}}_t^T={\parallel T-273.15\parallel }_2\times {10}^{-2}$.

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: $$\text{V}:x_i^j=\frac{\mathrm{\Delta }{x}_0}{2^j}(i-1),i=\mathrm{1,2},\cdots ,N_j+1\text{,}$$(11) where $x_i^j$ 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⁰ represents the size of the spatial grid at the sparsest level. It is obvious that $x_{i/2+1/2}^{j-1}=x_i^j=x_{2i-1}^{j+1}$.

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), $$u\left(x_i^j\right)=\frac{9}{16}\left(u\left(x_{i-1}^j\right)+u\left(x_{i+1}^j\right)\right)-\frac{1}{16}\left(u\left(x_{i/2-1/2}^{j-1}\right)+u\left(x_{i/2+3/2}^{j-1}\right)\right)\text{,}$$(12) where u stands for p, m, T, i = 1, 2, ⋯ N^j, and $x_i^j\in \text{V}2$.

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: $$u^{*}\left(x_i^j\right)=\frac{9}{16}\left(u\left(x_{i-1}^j\right)+u\left(x_{i+1}^j\right)\right)-\frac{1}{16}\left(u\left(x_{i/2-1/2}^{j-1}\right)+u\left(x_{i/2+3/2}^{j-1}\right)\right)\text{,}$$(13) where u stands for p, m, T, i = 1, 2, ⋯ N^j, and $x_i^j\in \text{V}1$.

Then, the error of all points in the V1 set are estimated through Equation (14). $${ϵ}_x(x_i^j)=\left|\frac{u^{*}(x_i^j)-u(x_i^j)}{u(x_i^j)}\right|\text{,}$$(14) where ϵ_x is the estimation error.

Compare the error ${ϵ}_x(x_i^j)$ to the tolerable error TOL_x: If ${ϵ}_x(x_i^j)\le {\text{TOL}}_x$, the point $x_i^j$ is not used on the (n + 1)th time level and is moved to the V2 set from the V1 set.If ${ϵ}_x(x_i^j)>{\text{TOL}}_x$, the point $x_i^j$ 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 $x_{i-1}^j$, $x_{i+1}^j$, $x_{2i-2}^{j+1}$ and $x_{2i}^{j+1}$.

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). $$p^{\ast }\left(x_i^j\right)=\sqrt{\frac{9}{16}\left(p^2\left(x_{i-1}^j\right)+p^2\left(x_{i+1}^j\right)\right)-\frac{1}{16}\left(p^2\left(x_{i/2-1/2}^{j-1}\right)+p^2\left(x_{i/2+3/2}^{j-1}\right)\right)}\text{.}$$(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 ${ϵ}_x\left(x_i^j\right)$ to the tolerable error TOL_x: If ${ϵ}_x\left(x_i^j\right)\le {\alpha }_x\cdot {\text{TOL}}_x$, the point $\left(x_i^j\right)$ is not used on the (n + 

1)th time level and is moved to the V2 set from the V1 set. If ${\alpha }_x\cdot {\text{TOL}}_x\le {ϵ}_x\left(x_i^j\right)\le {\text{TOL}}_x$, the point $\left(x_i^j\right)$ remains in the V1 set.If ${ϵ}_x\left(x_i^j\right)>{\text{TOL}}_x$, the point $\left(x_i^j\right)$ remains in the V1 set, and the neighboring points of $\left(x_i^j\right)$ at the jth and (j + 1)th grid level should be added to the V1 set, which denoted as $\left(x_{i-1}^j\right)$, $\left(x_{i+1}^j\right)$, $\left(x_{2i-2}^{j+1}\right)$ and $\left(x_{2i}^{j+1}\right)$.

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 ${\text{TOL}}_x\left(p(x_i^j)\right)=\left|p(x_i^j)\right|\times {10}^{-4}$, ${\text{TOL}}_x\left(m(x_i^j)\right)=\left|m(x_i^j)\right|\times {10}^{-4}$ and ${\text{TOL}}_x\left(T(x_i^j)\right)=\left|T(x_i^j)-273.15\right|\times {10}^{-4}$.

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 non-adaptive method, adaptive method and SPS. In the non-adaptive 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⁰ is 6 km.

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 shut-down 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⁶ Nm³/h. At t = 144 h, the outlet valve opening is reduced abruptly, and the outlet flow rate becomes 0.8 × 10⁶ Nm³/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.

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²· 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⁵ Nm³/h at the beginning, 0.5 × 10⁵ Nm³/h at 24 h and 0.1 × 10⁵ Nm³/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.

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, non-physical 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 non-adaptive method with a fixed time step Δt = 60 s and is larger than that obtained by the non-adaptive 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 non-adaptive 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 non-adaptive 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 non-adaptive 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.

4.4. Computational efficiency of adaptive method

To analyze the computational efficiency of the adaptive method, the CPU time of the non-adaptive 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 non-adaptive 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 non-adaptive 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 non-adaptive 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 non-adaptive 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 non-adaptive method is 6 times higher than that of the non-adaptive 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.

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 non-adaptive 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 non-adaptive 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

All Tables

All Figures