Simulation of Natural Gas Pipeline Networks Based on Roughness Optimization Algorithm and Global Mesh Refinement

Abstract

ABSTRACT

Natural gas pipeline network simulation technology is the fundamental technology of system capacity analysis, pipeline design, operation planning and optimization as well as emergency decision‐making for the whole life cycle of a given pipeline network system. There has been an increased demand for the computation efficiency and numeric accuracy of pipeline simulation with the increase of the total mileage of China's pipeline network as well as the centralize mode of one control center operating the entire system. This paper proposes a new numerical simulation method for natural gas pipeline systems based on roughness optimization and global mesh refinement. The numerical model of the natural gas pipeline system, which consists of governing equation of the pipeline fluid flow and characteristic equations of equipment, is firstly obtained by using an implicit finite difference method for discretization. The roughness identification problem is then transformed into an optimization problem by minimizing the error between measured and simulated values. The GA‐based algorithm is applied thereafter. Finally, a two‐step nonlinear iterative algorithm is proposed, which uses the coarse mesh to obtain the initial solution and the refined mesh to solve the problem to achieve accuracy and efficiency performance. The proposed method was verified by three industrial pipeline network examples. It is found that the average relative errors between the simulated and the measured data of the three cases are reduced by 3.87%, 5.06%, and 6.0%, respectively. The computational costs under 24‐h transient simulation conditions were reduced by 39%, 56%, and 65%, respectively. These numeric results show that the developed method has the advantages of stability, computation efficiency, and convergence, which provide a technical basis for the subsequent simulation of the national pipeline network.

Introduction

China's natural gas pipeline network system is an important energy transmission infrastructure [1, 2] which is composed of various elements such as gas sources, pipelines, customers, and storage facilities, forming a large-scale, structurally complex, and integrated hydrodynamic system. Based on natural gas flow dynamics in a pipeline network system, simulation technology is an efficient way for the comprehensive life cycle management and operation [3, 4]. This technology provides a quantitative assessment of the mass and heat transfer phenomena of natural gas in the pipeline infrastructure, accurately predicts the operation status of gas sources, customers, pipelines and equipment of the given pipeline network system.

With the rapid growth of the pipeline network, the demands on the computation cost and accuracy of pipeline simulation have been raised for high efficiency, energy saving and optimized operation of the pipeline network. Many studies have focused on the technology to accelerate simulation efficiency, including the model preprocessing method [5], Adaptive technology [6], discretization solution [7], and parallel computation [8, 9]. Wang et al. [10] proposed a natural gas pipeline network modeling method based on the concept of “divide-and conquer.” The method first solves all the connection points of the elements in the pipe network at one time and decouples the pipe network into several elements that can be solved independently. Fleming and Wang [11] uses large sparse matrices for parallel computation, which is about 10 times faster than serial computation.

Challenges have emerged due to the difficulty in accurately measuring the surface roughness of the pipelines with the presence of dirt, sediments and internal and external corrosion. Therefore, the calculation of the pipeline molecular resistance coefficient is unsettled which results in errors of natural gas pipeline simulations [12]. Abdolahi et al. [13] analyzed the effects of dominant factors of natural gas pipeline simulations using measured data from Iranian natural gas trunk lines. It was reported that the deviation between the simulated values and measured data can be minimized by adjusting the pipe roughness and heat transfer coefficient. Adopting the idea of a surrogate model, Yin et al. [14] simulated the transient model with a data-driven modeling and applied it for parameter identifications of the transient process. Leung et al. [15] used a error-in-variable model (EVM) in parameter estimation techniques to determine the surface condition of natural gas pipeline walls with a steady-state data set. However, it is difficult to be applied to the simulation of the conditions of pipeline transient operation.

Despite the implementation of correction at each time step, the oversimplification inherent in the mechanistic model diminishes the fidelity of the computational outcomes. Aimed at enhancing the computational efficiency and precision of natural gas pipeline network simulations, the present paper proposes a new method based on roughness optimization strategy and global mesh refinement. The methodology is delineated in Section 2, where a one-dimensional mathematical model is articulated, employing an implicit finite difference discretization scheme to model gas flow dynamics. A comprehensive implementation of the roughness optimization algorithm and the global mesh refinement process is provided. Subsequently, in Section 3, the efficacy of the proposed method is evaluated through the analysis of three distinct case studies. Concluding remarks are presented in Section 4, where the focus for future research is highlighted.

Methods

Mathematical Model

Pipeline Model

Considering the transient flow of natural gas in pipelines, a one-dimensional mathematical model consisting of the continuity equation, momentum equation, energy equation, and state equations is established based on the principle of the conservation of mass, momentum, and energy [16–18]. 1$\frac{\partial M}{\partial x}+A\frac{\partial \rho }{\partial t}=0$ 2$\frac{\partial }{\partial x}\left(A^{2}p+\frac{M^2}{\rho }\right)+A\frac{\partial M}{\partial t}+A^{2}g\rho \sin \theta +\frac{\lambda }{2}\frac{M^2}{D\rho }=0$ 3$\frac{\partial }{\partial x}\left[\left(h+\frac{M^2}{2A^2{\rho }^2}\right)M\right]+A\frac{\partial }{\partial t}\left[\left(h-\frac{p}{\rho }+\frac{M^2}{2A^2{\rho }^2}\right)\rho \right]+\frac{4AK\left(T-T_0\right)}{D}+Mg\sin \theta =0$ 4$\rho =\rho \left(p,T\right)$ 5$h=h(\rho ,T),$where M, ρ, x, t, A, p, g, D, θ, λ, h, T, K, and T₀, are gas mass flow (kg/s), density (kg/m³), position (m), time (s), cross-sectional area of the pipeline (m²), pressure (Pa), gravity acceleration (m²/s), inner diameter of a pipeline(m), the angle between pipeline and horizontal plane (rad), the friction coefficient, gas enthalpy (J/kg), temperature (K), total heat transfer coefficient of a pipeline (W/(m²·K)), and ambient temperature (K), respectively.

The accurate calculation of the friction coefficient in Equation (2) is the key to obtain the pipeline pressure. Colebrook's formula is recommended across a broad range of Reynolds numbers and roughness scales, and it is expressed as follows 6$\frac{1}{\sqrt{\mathrm{\lambda }}}=-2\lg \left(\frac{k}{3.7D}+\frac{2.51}{Re\sqrt{\mathrm{\lambda }}}\right),$where k is the pipe roughness (mm); Re is Reynolds number.

The BWRS equation is widely known for its accuracy and it is used to calculate properties of the natural gas physical as shown in Equation (4) [19].

Characteristic Equation of Nonpipe Components

Nonpipe components refer to the equipment, such as the compressors, valves, storage tanks, and other facilities that have an important impact on the fluid liquidity in the pipelines [20]. Taking a centrifugal compressor as an example, the compressor model including the pressure head equation, efficiency equation, pressure ratio equation, power equation, and temperature rise equation is established in Equations (7–11) [21]. 7$H=(1n{n}^2)\left(\begin{array}{l}{a}_1{a}_2{a}_3\\ a_4{a}_5{a}_6\\ a_7{a}_8{a}_9\end{array}\right)\left(\begin{array}{l}1\\ Q_{vol,i}\\ Q_{vol,i}^2\end{array}\right)$ 8$\eta =(1n{n}^2)\left(\begin{array}{l}{b}_1{b}_2{b}_3\\ b_4{b}_5{b}_6\\ b_7{b}_8{b}_9\end{array}\right)\left(\begin{array}{l}1\\ Q_{vol,i}\\ Q_{vol,i}^2\end{array}\right)$ 9$\epsilon =\frac{p_o}{p_i}$ 10$P_{power}=M\frac{H}{\eta }$ 11$T_o=T_i{\epsilon }^{\frac{k_g-1}{k_g\eta }}\frac{z_i}{z\left(p_o,T_x\right)}$where H, n, Q_vol,i, η, ε, p_i, p_o, P_power, T_o, T_i, T_x, k_g,and z_i are the adiabatic pressure head (kJ/kg), the compressor speed (rpm/min), volume flow rate at the compressor inlet (m³/s), adiabatic efficiency, compressor pressure ratio, inlet pressure (Pa), outlet pressure (Pa), shaft power (kW), compressor outlet temperature (K), compressor inlet temperature (K), outlet temperature after the ideal gas compression (K), natural gas insulation index and the compression factor, respectively. a₁…a₉, b₁…b₉ are the coefficients.

Junction Equation

In natural gas pipeline networks, junctions typically refer to the nodes connecting different pipeline segments or equipment. Junctions are the most common structure, which play an essential role in natural gas transportation, distribution, and scheduling. At each junction, the mass conservation equation is built based on the fact that the mass flow into and out of any node at any time must be equal to each other [21]. The energy conservation equation is based on the principle of the isenthalpic mixing principle. 12$\sum _{in=1}^{N_{in}}{M}_{in}=\sum _{out=1}^{N_{out}}{M}_{out}$ 13$\sum _{in=1}^{N_{in}}{M}_{in}{h}_{in}=\sum _{out=1}^{N_{out}}{M}_{out}{h}_{out}$where N, in, out are the number of nodes, the import node number and export node number, respectively.

The Boundary Conditions

The boundary conditions such as flow, pressure and temperature must be combined with the governing equations of pipe elements and nontube components:

The boundary conditions such as flow, pressure and temperature must be combined with the governing equations of pipe elements and nontube components: 14$\begin{array}{l}M=M\left(t\right)\\ P=p\left(t\right)\\ T=T\left(t\right)\end{array}.$

Discrete Method

The implicit center method is adopted for discretization of the variables in continuity and momentum equations. 15$\varphi =\frac{{\varphi }_{i+1}^j+{\varphi }_{i+1}^{j+1}+{\varphi }_i^j+{\varphi }_i^{j+1}}{4}$ 16$\frac{\partial \varphi }{\partial x}=\frac{{\varphi }_{i+1}^j+{\varphi }_{i+1}^{j+1}-{\varphi }_i^j-{\varphi }_i^{j+1}}{2\Delta x}$ 17$\frac{\partial \varphi }{\partial t}=\frac{{\varphi }_i^{j+1}+{\varphi }_{i+1}^{j+1}-{\varphi }_{i+1}^j-{\varphi }_{i+1}^j}{2\Delta t},$where ϕ represents M and p.

For T in the energy equation, a first-order windward scheme is used: 18$\left\{\begin{array}{l}\frac{\partial T}{\partial x}=\frac{T_i^{j+1}-T_{i-1}^{j+1}}{\Delta x},v_i>0\\ \frac{\partial T}{\partial x}=\frac{T_{i+1}^{j+1}-T_i^{j+1}}{\Delta x},v_i<0\end{array}\right..$

And the partial derivative of T with respect to time are approximated as follows 19$\frac{\partial T}{\partial t}=\frac{T_i^{j+1}-T_i^j}{\Delta t},$where the subscript i refers to the i-th computational node of each pipeline, and the superscript j refers to the j-th time layer.

Roughness Optimization

The roughness of the pipe is a key factor in determining the friction coefficient in Equation (6). However, the roughness of different parts of a pipeline undergoes different levels of changes as its usage time grows [21]. To address this issue, a new method is proposed to correct the segment roughness of the pipeline by using SCADA data along the pipeline (Figure 1).

[IMAGE OMITTED. SEE PDF]

Initially, taking a pipe as an example, each pipe is partitioned into a grid with N pipe segments. This segmentation assumes that the roughness of each individual segment is independent of others. Consequently, the roughness of each pipe is characterized by a distinct correction vector containing all segments, which is as follows 20$\mathbf{k}=\{{\mathit{k}}_1,{\mathit{k}}_2,\dots ,{\mathit{k}}_i,\dots ,{\mathit{k}}_N\},$where k represents the corrected vector of roughness.

The status vector at each segment is (M_i, p_i, T_i), where i ranges from 0 to N, constituting the status vector of each pipe, that is, X ={M_i, p_i, T_i}| i = 0,1,2,…, N. specially, i = 0 and i = N correspond to the pipeline's inlet and outlet boundaries, respectively. The input status vector of the pipe is defined as U = {p₀, T₀, p_N}, representing inlet pressure, inlet temperature, outlet pressure of each pipe. And the output of status vector of the pipe is defined as Y = {M₀, M_N, T_N}, representing inlet mass flow, outlet mass flow, outlet temperature of each pipe. The output status vectors of the pipeline system is obtained by solve the equations of pipes, equipment and junctions, which is defined as the status equation. And the status vector between the two adjoined time layers of j and j + 1 is expressed as, 21$\begin{array}{cc}\text{Status}\text{vector}& {\mathit{X}}^j={\left\{M,p_1,T_1,p_2,T_2,\dots ,p_i,T_i,\dots ,p_{N-1},T_{N-1},p_N,T_N\right\}}^j\\ \text{Corrected}\text{vector}& \mathbf{k}=\left\{{\mathit{k}}_1,{\mathit{k}}_2,\dots ,{\mathit{k}}_i,\dots ,{\mathit{k}}_N\right\}\\ \text{Input}\text{vector}& {\mathit{U}}^j={\left\{p_0,T_0,p_N\right\}}^j\\ \text{Output}\text{vector}& {\mathit{Y}}^j={\left\{M_0,M_N,T_N\right\}}^j\\ \text{Status}\text{equation}& {\mathit{X}}^{j+1}=\mathit{F}\left(\mathbf{k},{\mathit{U}}^j\right)\end{array}$

From Equation (21), the output vector Y is implicitly dependent on the input vector U and the corrected vector k. Furthermore, variables in the output vector Y are quantifiable variables. And Z = {M_N, M₀, T_N} is defined as the measurable vector. By minimizing the error of the measurable vector Z and output vector Y, a roughness optimization strategy is proposed as Equation (22). 22$\begin{array}{cc}\text{Objective function}:& \min \sum _{j=0}^m\{\parallel {\mathit{Y}}^{j+1}-{\mathit{Z}}^{j+1}\parallel \}\\ \text{Constraints}:& {\mathit{Y}}^{j+1}=\mathit{F}(\mathit{k},{\mathit{U}}^j,{\mathit{U}}^{j+1})\\ & {\mathit{U}}^j={\{M,p_1,T_1,p_2,T_2,\dots ,p_i,T_i,\dots ,p_{N-1},T_{N-1},p_N,T_N\}}^j\\ & {\mathit{U}}^{j+1}={\{p_0,T_0,p_N\}}^{j+1}\\ & {\mathit{Z}}^{j+1}={\{M_0,M_N,T_N\}}^{j+1}\\ \text{Optimising variables}:& \mathbf{k}=\{{\mathit{k}}_1,{\mathit{k}}_2,\dots ,{\mathit{k}}_i,\dots ,{\mathit{k}}_N\}\end{array}.$

The genetic algorithm method is applied to solve the optimal parameter search problem for roughness identification. The genetic algorithm is a maximum value problem and the objective function of Equation (22) is taken as the inverse of the fitness function. The coding of the correction variables is coded as real numbers. The minimum value is taken as 0 and the maximum value is taken as 2ε, ε is the pipe design roughness value. The roughness of some pipes is reduced by gas erosion. Oil deposits on the inner wall surface of some pipes lead to an increase in flow resistance, characterized as an increase in pipe roughness. The initial population size is 5 N and N is the number of correction variables. The selection operator uses the tournament operator. The crossover operator uses two-point crossover operator and the crossover rate is set to 0.8. The variation operator uses uniform variation operator and the variation rate is set to 0.03.

Global Mesh Refinement

A global mesh refinement is introduced in the solution procedure to speed up the computation. At the start of the simulation, each pipe can be discretized into coarse grids, followed by several iterations to obtain the initial solution. Once the initial solution vector is close enough to the optimal solution, a refinement grid is automatically reconstructed for precise computation. The criteria to start the global mesh refinement can be defined as: 23$\mathbf{F}\left(\mathbf{X}\right)<\xi \mathbf{F}\left({\mathbf{X}}^0\right),$where ξ is a smaller value, F(X⁰) is the residual value of the objective function with coarse meshes. The variables obtained in the coarse mesh simulation are served as the initial values of the fine mesh simulation. Linear interpolation is used to assign values to the fine grid nodes from coarse grid nodes to handle the inconsistency of spatial step sizes between coarse and fine grids.

Solution

To address the nonlinear partial differential equations presented in Section 2.1, the discretization techniques outlined in Section 2.2 are employed. This approach necessitates the transformation of the problem into an unconstrained optimization problem. The formulation of the objective function is as follows: 24$\mathbf{F}\left(\mathbf{X}\right)=\sum _{i=1}^N{f}_i\left(\mathbf{X}\right),$where f_i(X), i = 1, 2,…, N represents the residual of the i-th discrete equation, which includes equations of pipelines, equipment, junctions, and boundary conditions, X represents the variables to be solved. For a transient simulation, it is required that F (X^j, X^j + 1) = 0, where F, X^j + 1 represent the residual vector and the point vector at the j + 1 time step, respectively. Equation (23) can be transfered into the following expression, 25$\mathbf{J}\left({\mathbf{X}}^{j+1}\right)\Delta \mathbf{X}=-\mathbf{F}\left({\mathbf{X}}^{j+1}\right)\mathbf{,}$where J (X) represents the Jacobian matrix of Equation (23), $\mathbf{J}=\frac{\partial \mathbf{F}}{\partial \mathbf{X}}$. And the variables are updated as follows: 26${\mathbf{X}}^{j+1}={\mathbf{X}}^j+\mathrm{\beta }\Delta \mathbf{X}\mathbf{,}$where β is the step size along the search direction, and X can be solved by LU decomposition calculation formula (24). Updating X^j+1 until the convergence conditions for the objective function are satisfied, the solution vector at time layer j + 1 is obtained. Proceed with the aforementioned iterative process until the end of the simulation. Coupled with roughness optimization and global mesh refinement, the solution process of the new method can be summarized as shown in Figure 2.

[IMAGE OMITTED. SEE PDF]

Results and Discussions

To verify the validity and applicability of this method, three representative cases are conducted for tests. The computers used for the tests are equipped with Intel Core i7-7700HQ (2.80 GHz, quad-core processor), and 12GB of DDR2400 RAM.

Case 1: A Single Pipeline

Basic Data

The HM-YMS pipeline has a total length of 139.2 km, with one gas source and one user, a pipe with a diameter of 1016 mm, a wall thickness of 16.9 mm and a design pressure of 12 MPa. The pipe is made of X80 steel and has been in service for 12 years. The inlet and outlet of the pipeline are arranged with sensors for flow rate, temperature, and pressure. 60 measured data points, including the flow rate, temperature and pressure of the inlet and outlet of the pipeline, are selected and recorded during 1 h of actual operation, including the flow, temperature and pressure data of the inlet and outlet of the pipeline.

Discussion of the Results

Before implementing the genetic algorithm in solving the corrected pipe wall roughness, the HM-YMS pipeline is divided into 20 segments, and each segment takes (10 μm, 50 μm) as the search range of the pipe wall roughness. Then, the corrected values of the segmented roughness of the HM-YMS pipeline are solved, as shown in Table 1. Segment 12–14 of the pipeline is the low-lying segment of the whole pipeline, and the pipeline roughness is obviously lower than other positions. The reason for this phenomenon is the scouring effect of the pipe wall by the change of gas flow direction.

Table 1 Section roughness correction values of HN-YMS pipeline.

Segment ID	Roughness (μm)	Segment ID	Roughness (μm)	Segment ID	Roughness (μm)	Segment ID	Roughness (μm)
Section 1	25	Section 6	23.38	Section 11	24.002	Section 16	24.045
Section 2	23.5	Section 7	25.08	Section 12	20.613	Section 17	24.92
Section 3	23.5	Section 8	24.79	Section 13	20.808	Section 18	24.67
Section 4	23.8	Section 9	23.67	Section 14	19.897	Section 19	23.42
Section 5	23.85	Section 10	23.84	Section 15	23.968	Section 20	25.16

Figure 3 shows the comparison between the pressures at the pipe inlet with the measured value. The results show that the simulated pressures with the roughness optimization algorithm at pipeline inlet are consistent with the measured values, and the average relative error is 0.67%, which is 3.87% lower than that without the roughness optimization algorithm.

[IMAGE OMITTED. SEE PDF]

A transient simulation case is set up to compare computational speed with and without global grid optimization. The transient simulation is conducted by keeping the flow rate of the gas source changing with time. The simulation time is 24 h with a fixed timestep of 180 s. Table 2 shows the computational time consumption by the two computational methods. The results show that the computational speed is increased by 30% using the global mesh refinement.

Table 2 Comparison of the calculation times of the two methods.

Method name	Without global mesh refinement	With global mesh refinement
Method setting	Δx = 0.1 km	ΔxCoarse = 5 km	ΔxFine = 0.1 km
Time (s)	10	2	5
Total time (s)	10	7

Case 2: A Simple Pipeline Network

Basic Data

Case 1 verifies the feasibility in a single natural gas pipeline, but complex natural gas pipeline networks usually lack internal flow meters and only have flow meters at the pipeline inlet and outlet. Therefore, a pipeline network in the Southern Pipeline Network is used to analyze the applicability of the method under the condition of a lack of internal flow sensors.

The pipeline has a total length of 500 km and consists of three gas sources, four distribution points, 25 pipelines, and two pressure stations (Figure 4). To supplement the mass flow rate values of the pipeline network without flow sensors, the simulated value of the flow rate of the network is used as the real value, which is then utilized for the correction of the pipeline roughness.

Results and Discussion

[IMAGE OMITTED. SEE PDF]

Figure 5 illustrates the comparative analysis of the compressor outlet pressure, both with and without roughness optimization. It is indicated that a 1.11% enhancement is achieved in the accuracy of the compressor outlet pressure. Figure 6 presents a comparative study of the pressure distribution along the branch pipeline. An average improvement of 5.06% is revealed in the simulation accuracy of the pressure profile with roughness optimization.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Two scenarios of simulation durations of 1 h (Scenario 1) and 24 h (Scenario 2). The comparative analysis of computational time for simulations, with and without the implementation of global mesh refinement, is presented. Table 3 delineates the simulation durations for both computational methodologies. It is demonstrated that, for two transient simulation scenarios, reductions are achieved in computational expense by 44% and 56% achieves, respectively, with global mesh refinement.

Table 3 Comparison of the calculation times of the two methods.

Method name	Without global mesh refinement	With global mesh refinement
Scenario	1 h	24 h	1 h	24 h	1 h	24 h
Method setting	⃞Δx = 1 km	Δx = 1 km	⃞ΔxCoarse = 10 km	ΔxFine = 1 km	ΔxCoarse = 10 km	ΔxFine = 1 km
Time (s)	201	345	36	76	53	98
Total time (s)	201	345	112	151

Case 3: A Complex Pipeline Network

Basic Data

The applicability of the new method is validated for a single pipe and for those incorporating branch lines in Cases 1 and 2, respectively. Acknowledging the inherent complexity of natural gas pipeline networks, characterized by intricate interfaces and operational modes, Case 3 extends the analysis. The XY natural gas pipeline network consists of the XY-first Line, the XY-second Line, the XY-third Line and the LT Line, amassing a total length of 6100 km. This extensive network includes four principal gas sources, 15 distribution nodes, 206 valve chambers, 45 compressor stations, and 292 individual pipelines. Specifically, the second and third western lines are configured in parallel, with compressor stations interconnected via a network of pipelines and valves, thereby constituting a sophisticated pipeline system. The pipe is made of X70 steel and has been in service for 19 years. This configuration is schematically represented in Figure 7.

Discussion of the Results

[IMAGE OMITTED. SEE PDF]

Segmental roughness correction values for different pipelines in XY line are shown in Figure 8. It can be seen that for the same pipeline, the segmental roughness values are similar, while between different pipes the roughness values have obvious differences.

[IMAGE OMITTED. SEE PDF]

Figure 8 presents correction values of the segmental roughness for various pipelines within the XY Line. It is observed that, within the same pipeline, the segmental roughness values exhibit a degree of similarity, while significant variations appear in roughness values between different pipelines.

The application of the corrected values of pipe roughness, as depicted in Figure 8, is implemented for the current pipeline network. Figure 9 illustrates a comparative analysis of the inlet and outlet pressures before and after the roughness correction, with the measured pressures, as well as the corresponding errors for each compressor station along the XY Line. Without correction, the average relative deviations between the simulated and measured values for the inlet and outlet pressures are recorded at 7.8% and 6.0%, respectively. Implementing roughness optimization strategy, these deviations are significantly reduced to 1.8% for the inlet pressure and 0.8% for the outlet pressure, respectively.

[IMAGE OMITTED. SEE PDF]

Two Scenarios of simulation durations of 1 h (Scenario 1) and 24 h (Scenario 2) are conducted. Table 4 shows the simulation time used with and without global mesh refinement. The results show that the global grid optimization method saves 52% and 65% of the computational time consumption, respectively.

Table 4 Comparison of the calculation times of the two methods.

Method name	Without global mesh refinement	With global mesh refinement
Scenario	1 h	24 h	1 h	24 h	1 h	24 h
Method setting	⃞Δx = 1 km	Δx = 1 km	⃞ΔxCoarse = 10 km	ΔxFine = 1 km	ΔxCoarse = 10 km	ΔxFine = 1 km
Time (s)	295	483	43	98	47	123
Total time (s)	295	483	141	170

Engineers need to make time-consuming manual adjustments to the data set to arrive at a reasonable estimate of pipe roughness. The proposed roughness correction method reduces simulation errors by correcting pipe roughness through input data collected by the SCADA system. The current method is validated through practical cases.

However, there are some limitations in the current method. the SCADA acquisition data is used as the boundary condition of the simulation, and the quality of the data will directly affect the accuracy and convergence time of the simulation. Poor data will lead to nonconvergence of the simulation. Therefore, the quality of the measured input data has a significant impact on the correction algorithm to solve the true roughness parameters. Further research on data governance will be of great help in this area. Leung et al. [15] applied the EVM to determine the internal pipe wall roughness values in a pipeline network. This method allows for measurement error correction of the input data collected by the SCADA system. Combining genetic algorithm and EVM may have better processing results.

Pipe roughness parameters are usually determined in physical systems. The genetic algorithm treats these parameters as random variables with a given range of values. The final set of roughness correction values converges to the nearest minimum of the objective function. However, this may not represent the true actual roughness value of the pipe. Therefore, if there is some a priori knowledge of the system and its model parameters so that reasonable initial estimates can be obtained, the number of variables can be effectively reduced and the error of parameter correction can be reduced.

In addition, the number of parameter variables to be corrected increases with the number of pipes. The correction algorithm tends to fall into local convergence. Worse still, roughness correction is not necessarily required for some pipelines. Therefore, the pipe network can be divided into several small pipe networks or pipes according to the data collection points. The roughness correction algorithm is then executed.

Conclusion

Based on the implicit finite difference method, a new method is proposed for natural gas pipeline network simulation, with the roughness correction coupled with the measured data of SCADA system, and the global mesh refinement to improve the speed of the simulation. The accuracy of the proposed method is verified by using the actual pipeline examples, and the following conclusions are obtained:

• 1.

  The pipe roughness deviates from the design value because the pipe is subjected to internal erosion and corrosion. This leads to a large deviation of the pipeline simulation data. The roughness correction algorithm proposed in this paper improves the pipeline simulation accuracy. Three cases were considered to verify the proposed method, and it is observed that the average relative errors of simulated values with roughness optimization in the three cases are reduced by 3.87%, 5.06%, and 6.0%, respectively, compared with those without roughness optimization by comparing with the measured value from SCADA.

• 2.

  The increase of the pipeline size leads to a surge in the number of nonlinear systems of equations formed by the simulation, which reduces the computational speed of the simulation. A solution method is proposed to obtain the initial solution by coarse grid computation and fine solution by fine grid. Numerical experiments show that the computation time of the three cases is saved by 39%, 56%, and 65% respectively by applying a global mesh refinement strategy.

• 3.

  The analysis indicates that the pipe roughness is significantly reduced at locations where gas flow changes such as lows, elbows and tees in contrast to other regions. Consequently, to further improve the accuracy of pipeline roughness optimization, the uncertainty caused by pipeline corrosion and solid particle erosion should be considered, which will be investigated in the next phase.

• 4.

  These findings provide the following insights for future research: The internal diameter of the pipe for each grid characterizes the degree of clogging in the pipe. Each grid section plus a download value characterizes the occurrence of leakage at that location. A global optimization algorithm is used to solve for pipe fault detection using measured and simulated data as objective functions, and the inner diameter of each gridded pipe and the size of the download value of each cross-section as variables.

References

B. Hong, X. Cui, B. Wang, D. Fan, X. Li, and J. Gong, “Long‐Term Dynamic Allocation and Maintenance Planning of Modular Equipment to Enhance Gas Field Production Flexibility,” Energy 252 (2022): [eLocator: 123920].

B. Hong, X. Li, S. Song, S. Chen, C. Zhao, and J. Gong, “Optimal Planning and Modular Infrastructure Dynamic Allocation for Shale Gas Production,” Applied Energy 261 (2020): [eLocator: 114439].

M. Antonio and Z. Enrico, “A Framework for the Resilience Analysis of Complex Natural Gas Pipeline Networks From a Cyber‐Physical System Perspective,” Computers & Industrial Engineering 162, no. 1 (2021): [eLocator: 107727].

Z. H. Li and L. H. Yang, Pipeline Real‐time Data Integration and Pipeline Network Virtual Reality System: Digital Oil & Gas Pipeline: Research and Practice (Cham: Springer, 2021).

J. F. Helgaker and T. Ytrehus, “Coupling Between Continuity Momentum and Energy Equation in 1D Gas Flow,” Energy Procedia 26 (2012): 82–89.

E. Tentis, D. Margaris, and D. Papanikas, “Transient Gas Flow Simulation Using an Adaptive Method of Lines,” Comptes Rendus Mécanique 331, no. 7 (2003): 481–487.

P. Wang, B. Yu, Y. Deng, and Y. Zhao, “Comparison Study on the Accuracy and Efficiency of the Four Forms of Hydraulic Equation of a Natural Gas Pipeline Based on Linearized Solution,” Journal of Natural Gas Science and Engineering 22, no. 1 (2015): 235–244.

X. Wu, C. Li, Y. He, and W. Jia, “Operation Optimization of Natural Gas Transmission Pipelines Based on Stochastic Optimization Algorithms: a Review,” Mathematical Problems in Engineering 2018 (2018): 1–18.

G. Vitaliy and Z. Anatoly, “An Explicit Staggered‐Grid Method for Numerical Simulation of Large‐Scale Natural Gas Pipeline Networks,” Applied Mathematical Modelling 65, no. 1 (2018): 34–51.

P. Wang, B. Yu, D. X. Han, D. L. Sun, and Y. Xiang, “Fast Method for the Hydraulic Simulation of Natural Gas Pipeline Networks Based on the Divide‐And‐Conquer Approach,” Journal of Natural Gas Science & Engineering 50, no. 1 (2018): 55–63.

Graham Fleming and Qinghua Wang, “A Parallel Solution For Large Surface Networks in a Fully Integrated Reservoir Simulator,” in SPE Reservoir Simulation Conference (Galveston, Texas: 2017) SPE‐182634‐MS.

C. J. Li, W. L. Jia, and X. Wu, “A Steady State Simulation Method for Natural Gas Pressure‐Relieving Systems,” Journal of Natural Gas Science & Engineering 19 (2014): 1–12.

F. Abdolahi, A. Mesbah, R. B. Boozarjomehry, and W. Y. Svrcek, “The Effect of Major Parameters on Simulation Results of Gas Pipelines,” International Journal of Mechanical Sciences 49 (2007): 989–1000.

X. Yin, K. Wen, W. Huang, et al., “A High‐Accuracy Online Transient Simulation Framework of Natural Gas Pipeline Network by Integrating Physics‐Based and Data‐Driven Methods,” Applied Energy 333 (2023): [eLocator: 120615].

T. Leung, C. Hartloper, J. Smith, K. K. Botros, H. Golshan, and T. Glen, “Application of EVM to Determine Pipe Internal Surface Conditions in Complex Natural Gas Lateral Networks,” Journal of Natural Gas Science and Engineering 24 (2015): 217–227.

F. E. Uilhoorn, “State‐Space Estimation With a Bayesian Filter in a Coupled PDE System for Transient Gas Flows,” Applied Mathematical Modelling 39, no. 2 (2015): 682–692.

R. L. Meira, G. M. N. Costa, R. A. Kalid, and M. A. F. Martins, “Improving the Centrifugal Compressor Map Through Rigorous Thermodynamic Modeling: an Analysis on a Natural Gas Compression Station Pipeline,” Journal of Natural Gas Science and Engineering 92, no. 1 (2021): [eLocator: 104006].

I. Durgut and M. K. Leblebicioğlu, “State Estimation of Transient Flow in Gas Pipelines by a Kalman Filter‐Based Estimator,” Journal of Natural Gas Science and Engineering 35 (2016): 189–196.

H. Y. Ren, C. J. Li, and G. Xie, “Evaluation of State Equations of Natural Gas in Pipeline Transportation,” Advanced Materials Research 463–464 (2012): 936–939.

J. R. Wang, T. Wang, and J. Z. Wang, “Application of π Equivalent Circuit in Mathematic Modeling and Simulation of Gas Pipeline,” Applied Mechanics and Materials 496 (2014): 943–946.

A. J. Wiid, J. D. le Roux, and I. K. Craig, “Modelling of Methane‐Rich Gas Pipeline Networks for Simulation and Control,” Journal of Process Control 92, no. 6 (2020): 234–245.