{'fraction_non_alphanumeric': 0.04016637851111278, 'fraction_numerical': 0.01866746849081661, 'mean_word_length': 4.2721268163804496, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An explicit staggered-grid method for numerical simulation of large-scale natural gas pipeline networks

We present an explicit second order staggered finite difference (FD) discretization scheme for forward simulation of natural gas transport in pipeline networks. By construction, this discretization approach guarantees that the conservation of mass condition is satisfied exactly. The mathematical model is formulated in terms of density, pressure, and mass flux variables, and as a result permits the use of a general equation of state to define the relation between the gas density and pressure for a given temperature. In a single pipe, the model represents the dynamics of the density by propagation of a non-linear wave according to a variable wave speed. We derive compatibility conditions for linking domain boundary values to enable efficient, explicit simulation of gas flows propagating through a network with pressure changes created by gas compressors. We compare Kiuchi's implicit method and an explicit operator splitting method with our staggered grid method, and perform numerical experiments to validate the convergence order of the new method. In addition, we perform several computations to investigate the influence of non-ideal equation of state models and temperature effects into pipeline simulations with boundary conditions over various time and space scales.


Introduction
Recent years have witnessed a renewal of interest in efficient methods for accurate simulation of transient flows through pipeline systems [5,9,27,2,12]. The physical models and mathematical descriptions of fluid flow in pipelines under various regimes are well-understood [28], and a variety of available numerical methods have been developed over decades [26,15]. In general, the Euler equations in one dimension are used to describe compressible fluid flow in a pipe using two hydrodynamic partial differential equations (PDEs) and an equation of state that relates gas density to pressure. The Darcy-Wiesbach formula is used to approximate the nonlinear energy dissipation caused by three-dimensional turbulent flow within the pipe, while diffusive effects of turbulence are ignored.
A range of models and numerical methods is available for solution of the initial boundary value problem (IBVP) of gas flow through a pipeline system, depending on the time and space scales of interest. A variety of high-fidelity, explicit methods are used to accurately approximate the effects of rapid variations in pressure and flow on time-scales corresponding to the speed of wave propagation [30,4,14]. These methods solve the full isothermal Euler equations and account for nonlinear self-advection, so that the formation and propagation of shock waves can be captured. They are typically second order accurate in regions where pressure and flux are continuous. These methods are designed for accurate solution over short time scales and are impractical for simulation of large gas pipelines in the typical operating regime in which significant boundary changes occur on time scales of hours and longer. Other studies have evaluated non-isothermal models [22,7,1]. On the other hand, several lowerfidelity, implicit methods are better suited to characterizing dynamics of large pipelines in the operating regime where pressure and mass flow change on time scales much longer than the rate of wave advection [21,17,16,13]. The terms in the Euler equations that may be omitted in this regime of interest have been characterized in previous empirical studies [8,21], and the result is a nonlinear first order hyperbolic PDE system. Such implicit methods are efficient, scalable, and adaptable to larger systems with complex network structure and nodal components such as gas compressors and regulators. However, even slight changes on fast time-scales can present obstacles to stability of these methods [16].
A major aim of recent work has been to bridge the gap between high-fidelity explicit methods that are computationally costly but stable for very fast transients [1], and lower accuracy schemes that use simplified physics but can be applied to large-scale pipeline networks [13]. One recently proposed method applies the well-known operator-splitting scheme to solve the simplified equations that describe one-dimensional dynamics of gas flow in networks of long pipes [10]. This method accurately represents rapid, small amplitude changes while resolving slowly-varying dynamics on the time scale of normal gas pipeline operations, and is extensible to complex pipeline networks with nodally-located time-varying compressors that boost gas flow into adjoining pipes. Although it is unconditionally stable and second order accurate, the operator-splitting discretization is based on the propagation of characteristics through a uniformly-spaced grid, which requires the assumption of uniform and constant wave propagation speed, and thus an ideal gas model. Because high pressure transmission pipeline flows exhibit highly non-ideal gas states, an explicit method for simulation of pipeline flows using non-ideal models is of significant interest.
In this paper, our aim is to develop a computationally efficient and explicit numerical method that is capable of simulating gas pipeline dynamics with non-ideal equation of state models. We formulate a model for gas transport in a large-scale pipeline network. For a single pipe, solution of an IBVP for this model could be viewed as representing the propagation of a non-linear wave in the density of the gas. The non-linearity is present in the wave speed and in the non-linear dissipation term. We present a new second-order accurate staggered grid discretization, with finite differences in space and time, for this model on a network and compare the resulting IVP solutions with other numerical schemes including that of Kiuchi [16], explicit operator splitting [10], and implicit lumped elements [31].
We focus in particular on the connection between the staggered grid approach and the implicit method of Kiuchi, which uses a standard (non-staggered) discretization where the staggered quantities are approximated using trapezoidal rule averages. The implicit method of Kiuchi is unconditionally stable even when using an arbitrarily large time discretization step relative to the spatial discretization, although accuracy is reduced in this case. Kiuchi's implicit method requires the solution of a system non-linear equations for each time step, which arise from the dissipation term and the generally complicated form of the compressibility factor in the non-ideal case. Although our proposed explicit staggered finite difference discretization is restricted by a Courant-Levy stability condition, it employs simple time stepping that does not require solution of a large-scale implicit system. Moreover, the data dependence on each time step is highly localized, which facilitates parallelization of computation of IBVP solutions for very large-scale (e.g. continental) pipeline networks.
In Section 2 we present a one-dimensional mathematical formulation for gas transport in a long pipe and various possibilities for modeling the dependence between gas density, pressure and temperature using a compressibility factor function. We present the derivation of compressibility factor formulas in Appendix A. In Section 3, we derive an explicit staggered grid discretization for non-ideal gas flow on a single pipe, then validate its second order accuracy, and derive a stability criterion for the grid steps. In Section 4 we develop standardized compatibility rules for boundary conditions at junctions where two or more individual pipes are linked, so that the scheme can be applied to complex pipeline networks with gas compressors. In Section 5 we present numerical experiments on a single pipe to verify that the method is better than second order accurate, as well as to evaluate the significance of non-uniform temperature and non-ideal gas modeling for simulation of flows on various time and space scales. In Section 6 we examine a number of numerical experiments on a network of pipes using the proposed method as well as the operator splitting method and a lumped-element implicit method. Finally, in Section 7 we conclude with a discussion and open questions for future work.

Mathematical Modeling of Pipeline Dynamics
The Euler equations in one dimension describe compressible fluid flow in a long pipe using two hydrodynamic PDEs and an equation of state that relates gas density to pressure [15]. The Darcy-Wiesbach formula is used to approximate the nonlinear energy dissipation caused by three-dimensional turbulent flow within the pipe. Because the effects of turbulent diffusion are ignored in that approximation, these equations can be considered as a phenomenological average of the three-dimensional Navier-Stokes hydrodynamic equations that captures effects over spatial scales much longer than the pipe diameter and over time scales longer than those required for dissipation of acoustic waves. The rapid dissipation of turbulent effects in this regime of interest precludes the formation of shockwaves, which are strongly nonlinear acoustic events that develop and dispel over much shorter time scales as a result of significant abrupt changes [28,24].
The Euler equations can be stated as two PDEs that describe the conservation of mass and momentum, respectively, as well as an equation of state. The PDEs are The variables are gas velocity v, pressure p, and density ρ, defined on a domain x ∈ [0, L] at time t. The parameters are the friction factor λ, pipe diameter D, gravitational acceleration g, and pipe angle θ. The terms on the right hand side of (2) aggregate friction and gravity effects. We assume that gas pressure p and density ρ satisfy the equation of state relation where R and T are the ideal gas constant and temperature, respectively. Here Z(p, T ) is a compressibility factor that defines the deviation of the equation of state from the ideal gas law. We will assume that the equation (3) can be uniquely solved for pressure as a function of density. Equation (2) is valid in the regime when changes in gas consumption and injections are sufficiently slow to not excite propagation of sound waves. Formally, the term ∂ t (ρv) on the left hand side is much smaller than ∂ x (ρv 2 + p). The ratio of the pressure gradient term ∂ x p to the term ∂ t (ρv) is typically on the order of 1 : 0.01, [21]. In addition, the flow velocities are much smaller than the nominal speed of sound, c = √ ZRT , where Z is a nominal compressibility factor. In this regime, the gas advection term ∂ x (ρv 2 ) can be omitted, even in comparison with the already small term ∂ t (ρv). For simplicity, we assume that the pipeline is level, thus ignoring the gravity term ρg sin θ on the right hand side of (2). We will usually assume that the gas temperature is uniform along the pipeline, so that T is constant. One exception to this will be the Section 5.4, where we consider the effects of of heat generation by a compressor station.
With the above assumptions, (1) and (2) can be written in terms of (per area) mass flux φ = ρv (which has units kg/m 2 /s) as The gas dynamics on a long pipe segment of a large gas transmission pipeline are commonly represented using (4)-(5), subject to the initial and boundary conditions [16,30]. Simple boundary conditions on each end include one of ρ(t, 0) = ρ 0 (t) or φ(t, 0) = φ 0 (t), and one of ρ(t, L) = ρ L (t) or φ(t, L) = φ L (t). More complicated boundary conditions can arise when we consider junctions of several pipes, as well as devices such as compressors and regulators that are used to actuate flow.
The compressibility factor Z(p, T ) is typically given as a formula with parameters that have been fitted to measured data obtained during early engineering studies [6,11]. A simple and widely-used formula is the CNGA method [19,18], given by where a 1 = 344400, a 2 = 1.785, and a 3 = 3.825 are nominal coefficients and G is the gas gravity (with air = 1 as reference). Here p is pressure in units of psig (gauge pressure), and T is in units of • R (degrees Rankine). We consider a nominal gas gravity of G = 0.650784 for 80% methane and 20% ethane. The compressibility factor relation formula (6) can be used with the equation of state (3) to create a bijection between density and pressure. We let P : R ≥0 → R ≥0 denote a mapping from density to pressure acting p = P (ρ), and let P −1 be the inverse, acting ρ = P −1 (p). This mapping is derived for several cases of physical assumptions in Appendix A.

Numerical discretization scheme for a single pipe
We now describe a numerical scheme for finite difference (FD) discretization of the dynamic equations (4-5) on a single pipe. The proposed scheme is motivated by the fact that these PDEs reduce to the linear wave equation under the ideal gas assumptions for the constitutive law and zero friction, i.e. Z(p) ≡ 1 and f ≡ 0. Therefore, we will adopt the classical second order accurate staggered discretization for a wave equation written in a mixed form. This discretization could be considered a prototype of the classical Yee scheme for time domain computational electrodynamics [29]. Consider the staggered FD discretization shown in Figure 1, where space grid points are indexed by two indices i and j = i + 1 2 and time layers are indexed by two additional indices n and m = n + 1 2 . Variables ρ and p will be discretized on the staggered grid in "blue" points (circles), while the flux variable φ will be discretized on the staggered grid in "green" points (squares). The equation (4) will be discretized, at a point (m, i), as The equation (5) will be discretized at a point (n + 1, j), as where β = λ/(2D). Observe that in the fraction on the RHS of (8), ρ and φ are evaluated between staggered grid points where their values are not available. Therefore, we approximate the RHS of (8) by Time, t Figure 1: Illustration of the staggered grid discretization. The pressure, p, and density, ρ, are approximated at points x i and t n (marked by circles). The flux, φ, is approximated at points x j and t m (marked by squares).
By substituting (9) in the equation (8) we get a nonlinear implicit equation Although the LHS is implicit in φ m+1 j , (10) needs to be solved only pointwise, i.e. there is no coupling between φ m+1 j and φ m+1 j+1 . The function that needs to be inverted to solve for φ m+1 From the form of the RHS in (11) one can see that the function F (x) is a bijection on the real line. Indeed, both factors in the RHS of (11) are increasing for x > 0 and change sign at x = 0. Hence, we make the following trivial, but useful observation Using this observation we can write the inverse of F (x) = y explicitly using the quadratic formula where a = β t ρ i+1 n + ρ i+1 n+1 > 0 and therefore sign(y)ay = a|y| ≥ 0.
From (12) and (14) we conclude that only the positive solution should be used in (13). Thus, the inverse of F (x) has the closed form expression Remark 1. Note that in the limit of vanishingly small friction coefficient λ equation (10) transforms into

Discretization order analysis
It is critical to observe that the discretization (7-8) conserves mass in the system, as evident from (7). Indeed, summing (7) for all space points i, one obtains where φ m 0 and φ m N are the left-most and right-most values of the mass flux, respectively. The discretization (7-8) is second order accurate in space and time, which can be validated by verifying that each of the terms in the discrete equations (7-8) constitutes a second order accurate approximation to the corresponding terms in the PDEs (4-5). We provide this validation below by showing that the nonlinear term (9) is approximated with at least second order of accuracy. Specifically, we aim to show that We refer to Figure 1 for the illustration of the staggered grid. First, let us consider the numerator and the denominator of (17) independently, and show that they constitute a second order approximation. Consider a Taylor expansion of ρ(t, x) at points (t n+1 , x i ) and (t n+1 , x i+1 ), which yields where where ξ i,j and ξ i+1,j are some points in the intervals (x i , x j ) and (x j , x i+1 ), respectively.
and taking an average of (18) and (19) yields For the numerator, we assume that φ does not change sign and treat φ |φ| as one function. As above, we then perform a Taylor expansion in time, similar to (18)(19), and obtain Substituting (20) and (21) into the LHS of (17) and performing a Taylor expansion in the correction terms O( x 2 ) and O( t 2 ) we obtain the RHS of (17).

Stability analysis
Although the discretization of the friction term (9) causes equation (8) to appear implicit, there is actually no spatial coupling between the φ m+1 terms as seen in (15), and therefore one should expect only conditional stability. Furthermore, the nonlinearity of the PDE (5) precludes the performance of a Von-Neumann type analysis. Instead, we will analyze what happens with a perturbation to a numerical solution on one time step and derive some error estimates that will include t, x and solution properties. The stability criterion that we derive below provides the condition under which we will be able to guarantee that the perturbation is not amplified so greatly as to dominate the solution.
To this end, we assume that (ρ n i , p n i , φ m j ) and (ρ n i +ρ n i , p n i +p n i , φ m j +φ m j ) are two solutions of the numerical scheme, and (ρ n i ,p n i ,φ m j ) is the perturbation. In order to maintain simplicity of notation, we introduce the discrete operators Because the underlying PDE (4) and the corresponding discrete equation (7) are linear, the perturbation satisfies the similar discrete equation On the other hand, because the underlying PDE (5) and its discrete form (8) are nonlinear, the perturbation (in general) is not itself a solution of the numerical scheme. Instead, we have where we assumed that φ > 0 and φ +φ in the domain of interest. This is done mostly to simplify the notations, but also because this is where the physical assumptions of the model reductions are valid in order to obtain PDEs (4)(5). With regard to (23), we note thatp does not satisfy the same equation as p. Instead, we havẽ where P (ρ) is the point-wise dependence of p on the value of ρ. Assuming the perturbations are small, we neglect higher order terms in (24). The pressure term in (23) becomes By omitting higher order terms, the RHS of (23) simplifies to Again, neglecting higher order terms inφ andρ we obtain Summing up, the perturbations (ρ n i ,p n i ,φ m j ) satisfy equations (22) and (27), which correspond to discretization of wave equation with a variable advection speed P (ρ) n+1 j , and with the forcing term in the RHS of (27). In fact, P (ρ) n+1 j is the local speed of the non-linear equation as well. Thus, the stability condition on the time step is formulated in terms of this wave speed: Because the density changes as a function of space and time, we require this condition to be satisfied for the range of expected densities ρ, see e.g. method of frozen coefficients [25].

Treatment of boundary conditions on a pipeline network
In this section we describe a standardized approach for construction of pipeline network models with gas compressors for the purpose of simulating transient flows over large spatial scales. This reduction extends the previous modeling work [31,32,20]. We represent a pipeline network a metric graph with a system of PDEs on each edge, subject to timedependent boundary conditions and nodal compressor actions, in a manner similar to what has been done for quantum graphs [3]. The graph is defined by a set of nodes (or vertices) V and a set of directed edges E ⊂ V × V with elements (k, l) ∈ E that represent pipes between the nodes k, l ∈ V. Every edge (k, l) ∈ E is associated with a diameter D kl , length L kl , cross-section area S kl , and friction coefficient λ kl , which constitute a metric on E. Gas flow through each pipe (k, l) ∈ E is characterized by the hydrodynamic equations (4)-(6) with solution denoted by ρ kl , p kl and φ kl , and where p kl = P (ρ kl ). Here P (·) is the map defined in one of (55) or (60). For the pipe (k, l), the sign of gas flow from node k to node l is denoted as positive by convention. In our notation, we define the collection of pipes incoming to or outgoing from a node l ∈ V by The balance of incoming and outgoing physical mass flows and injection of gas into the network must be preserved for each network node l ∈ V. Such mass flow balance is maintained by applying the Kirchhoff-Neumann boundary conditions at network nodes, at which gas may also be injected or withdrawn from the system. We denote the mass flow rate of gas withdrawn from the network at a node l ∈ V by the time-dependent function q l . The balance condition can then be written as where sgn kl = 1 if flow φ kl through (k, l) ∈ E is denoted as incoming into l, and sgn kl = −1 if flow φ lk through (l, k) ∈ E is denoted outgoing from 1. Note that because mass flux φ kl is expressed in per area units, while gas withdrawal q l is in mass flow units, the cross-section areas S kl are required to weight mass balance in the above equation. That is, φ kl (t, L kl ) is the (per area) mass flux at the end of pipe (k, l), and S kl φ kl (t, L kl ) is the total mass flow into node l from pipe (k, l) if sgn kl = 1. Similarly, φ lk (t, 0) is the (per area) mass flux at the start of pipe (l, k) and S lk φ lk (t, 0) is the total mass flow out of node l into pipe (l, k) if sgn kl = −1. Gas flow through the network is controlled by compressors that boost pressure while preserving mass flow. In our model, such compressors are located at the interfaces between nodes and adjoining pipes, and are able to boost flow going from a node into an adjoining pipe, see Fig. 4. This can be generalized by where α ≥ 1 is a compression ratio that relates the suction (inlet) and discharge (outlet) pressures p in and p out , respectively. As with the other boundary conditions, we reformulate this condition in terms of density according to To model the action of such a compressor at the boundary of a pipe, we define auxiliary nodal pressure variables p k and p l so that the function of compressors at the endpoints of a pipe (k, l) can be described by where α kl α kl are compression ratios that relate nodal pressures to boundary pressure at the inlet and outlet of the pipe, respectively. In the method presented here, we assume that no gas is removed from the system by compression processes. Complete boundary conditions can be provided for the pipeline network by specifying flow withdrawal or pressure at each node. If flow q l removed from the network is specified, the pressure at that node is a dependent variable, and if pressure p l at a node is given, flow into or out from the network at that node is a free variable. The pressure at at least one so-called "slack" node must be specified. This endpoint boundary information, together with initial conditions, defines a well-posed IBVP for a pipeline network. The equations (30) and (33) together constitute nodal compatibility laws that provide consistent boundary conditions for all the pipes in the network. Next, we describe how the staggered grid discretization scheme can be applied to link pipes together with these compatibility conditions.

Boundary conditions for a single pipe
The most natural boundary conditions are on the values of the flux φ, which are implemented for a single pipe simply by assigning the value to the flux at the ends. Alternatively, the boundary conditions on pressure p and density ρ are reformulated at each time step by computing the appropriate value of the mass flux at that step that would result in those values of pressure and density, respectively. Thus, for each pipe, each time step consists of the following sub-steps: 1. Compute the new value of the flux inside the domain.

Compute the appropriate values for the flux at the ends of the domain (left-most and right-most values)
. This is the only step that incorporates the boundary conditions.
3. Compute the values of the density (and hence pressure) at the next time step.
In case that the boundary conditions specify values of mass flux at an endpoint of a pipe, step 2 above is straightforward. However, if the boundary conditions specify the value of density at a pipe endpoint, the appropriate value of the flux is found from the conservation of mass equation (7). For example, if φ j m = φ kl (t m , L kl ), i.e., φ j m gives the right boundary value for the mass flux on a pipe (k, l), then we have at time t m where φ m j−1 = φ kl (t m , L kl − x) is the internal value computed in step 1, ρ n i = ρ kl (t n , L kl ) is the density at a known time t n , and ρ n+1 i = ρ kl (t n , L kl ) is specified by the boundary conditions. In case the boundary conditions specify the value of pressure, this condition is reformulated in terms of the value of density using the map P (·), and we proceed as above.

Boundary conditions at a node connecting two pipes
In this section we describe incorporation of the multiplicative pressure gain condition (31), which is used to model the function of a gas compressor.
We first reformulate the condition on density as a condition on the mass flux. Here we note that in order to satisfy mass conservation properties we must have a balance of mass flow through the joint. This is expressed in terms of per area mass flux variables by where S in , D in and S out , D out are the cross sectional areas and diameters of the pipes entering and leaving the joint (node), respectively. However, we still can choose the value of the flux. If φ j m is the value to be found, we then have the conditionss where φ m j−1 and φ m j+1 are found during the step 1, as described in   Thus, the pressure gain boundary condition (31) is implemented by solving (32-36) for the value of the flux φ m j , for which we have the two expressions From (37-38) we can obtain a solution for p in using the relation Note that the solution of (39) is unique when the function P (·) is one-to-one. In particular, this is the case when P (·) (and hence P −1 (·)) is a monotone function. For example, in the case (59) the equation (39) is a quadratic equation for p in that has a unique positive solution.

Boundary conditions for a joint with multiple pipes
In this section, we consider a generalization of the pressure gain boundary condition discussed in the previous section to the case of three and more pipes. For generality, we also assume that compressors are present at the endpoints of every pipe (k, l) ∈ E. Suppose that multiple pipes are adjoined at a node k ∈ V, which is associated with a unique pressure p k . Thus, the pressure gain condition will be written in the form p kl (t, 0) = α kl (t)p k (t) and p kl (t, L kl ) = α kl (t)p l (t), where the boost ratios α kl and α kl represent the actions of compressors propelling gas into the pipe from the endpoints. Here, p kl (t, 0) and p kl (t, L kl ) represent the pressures at the endpoints of the pipe (k, l), and p k and p l represent the pressures at the nodes on the other side of the compressors located at the endpoints. In this general setting, at least one of α kl and α kl are set to unity at any given time, because both compressors would not be used simultaneously.
The conservation of mass condition for the joint is given by the flow balance equation (30). As done in equation (37), applying the discretization to the conservation of mass equation, we obtain where φ m kl− is the flux at the node adjacent to the end node φ m lk (equal to φ kl (t m , L kl − x kl ) if sgn kl = 1 and φ kl (t m , x kl ) if sgn kl = −1); ρ n kl is the known value of density at the k node end of the pipe (k, l) and x kl is the spatial discretization step for the pipe (k, l). Multiplying (41) by S kl and summing up over all pipes going into the node and using the balance of flows equation (30), we get the new version of (39): This is an equation on the node pressure p n+1 l . All other quantities in (42) are known. Once p n+1 l is solved for, from (41) one finds the values of φ m kl , which are used as the boundary conditions on each of the pipes separately.
Remark 2. In this paper we consider only the case where the time discretization step t is the same for all pipes (k, l), while the spatial discretization step x kl may vary. Also note that on each pipe we have to satisfy an appropriate CFL condition relating x, t and the wave speed.

Numerical experiments on a single pipe
In this section we perform a number of numerical experiments designed to compare and contrast various models. One of the comparisons we perform is based on the choice of the constitutive law for Z(p). The two models to be considered are where b 1 = 1.00300865, and b 2 = 2.96848838 × 10 −8 .

Verification of the order of the method
In order to verify the order of convergence of the method we considered a single initial value problem on one pipe (with length L = 10 4 meters) solved numerically with different resolutions ( t = 1, 3 −1 , 3 −2 , . . . , 3 −5 ) but fixed time and space resolution x/ t = 454.55 ratio. The reason for this choice of the resolution ratio is to keep it larger than the maximum possible speed of sound in the pipe. As a reference we also use a numerical solution with a much higher resolution ( t = 3 −6 ). We chose to do a 3× refinement instead of the more standard 2× so that the discretization points for the pressure p and the density ρ for a coarser descritization are a subset of fine discretization. For 2× refinement, this is not the case.
We also note that in order to set up the discrete initial value problem we need to specify the values of ρ and φ which are defined at different times (t = 0 and t = t/2). We set up all problems with the same initial values of density which has an arctan() step profile centered at the center of the pipe and shifted upward by a mean value ρ. We then set up the fine-grid reference problem with the initial value for the flux The initial conditions for the flux φ on coarse grids are obtained from the solution for φ on the fine grid (by restriction) at the time t coarse /2.
These initial conditions were selected because in the case of a simple ideal gas model Z(p) ≡ 1 and zero friction f = 0 we get a traveling wave solution moving in positive direction with the speed c ref . The errors in ρ, p and φ are all computed after a single time step so as to approximate the L 2 -error and plotted on figure 4, e.g. the L 2 error for p is where in practice p exact is the high and p approx is a low resolution solution.
The convergence rates computed based on the last two data points as well as the first and last data points are shown in the table below. They indicate that all three variables converge at least with second order.  Table 1: Numerical estimates of order of accuracy shown in Figure 4.

Fast transients
Consider a 20km long pipe with diameter 0.9144m and friction coefficient 0.01 and gas at temperature of 288.706 K. The initial conditions are uniform throughout the pipe with zero flux φ(0, x) = 0, constant pressure p(0, x) = 6.5MPa, and corresponding density 56.816kg/m 3 . The boundary conditions are: • On the left boundary constant pressure is maintained, p(t, x L ) = 6.5MPa.
• On the right boundary the value of the flux is specified: for t < 10 in min, φ 0 , for 10 < t < 30 in min, φ 0 /10, for 30 < t in min, This initial value problem is descretized using the non-ideal gas model (44) and the ideal gas model (43). For the latter, the speed of sound is taken to be 338.25 m/s, so as to match the initial conditions for both pressure and density. This corresponds to a constant value of Z = 0.83616. The solutions for the IBVP on the left and right boundaries are shown on Figure 5. As could be seen from Figure 5 the ideal gas model (43) results in significantly smaller values of the maximum velocity on the right boundary compared to the non-ideal gas model (44). The smaller values of the velocity are due to significantly larger values of the density around t = 30min.
On the left boundary the ideal gas model produces somewhat larger values of the flux in the range 12 − 30min and delayed velocity profile in the times 32 − 40min.

Slow transients
Consider a 50km long pipe with diameter 0.9144m and friction coefficient λ = 0.01 at temperature 288.706 K. The initial conditions are taken to be uniform throughout the pipe with flux φ(0, x) = φ 0 = 240 kg m 2 s , constant pressure p(0, x) = 6.5MPa and density of 56.816kg/m 3 . The boundary conditions are as follows.
• On the left boundary, pressure is a slowly oscillating harmonic shifted by a constant p(t, x L ) = 6.5 (1 + .25 sin(πt/t scl ))mPa, where the time scale t scl = 6 hours.
• On the right boundary, the flux is constant φ(t, x R ) = φ 0 .
The initial value problem is solved using the non-ideal gas model (44) as well as the ideal gas model (43). For the latter, the speed of sound is taken to be 338.25 m/s, in order to match the initial conditions for both pressure and density. The solutions are shown after 50 periods on Figure 6.
During the first t 1 = 4hrs the boundary conditions are During this time the solution approaches a steady state. After t = t 1 we apply oscillatory boundary conditions as follows: where t scale = 12hrs. We consider two steady models of temperature distribution modeling heating of gas by the compressor. Both models assume exponential temperature decay to a steady value, but with different rates Here we took T ambient = 288.706K and T jump = 40K. In the first model it is assumed that the temperature spike decays to the ambient level on the scale of one kilometer, while in the second model the decay is assumed to happen on the scale of ten kilometers. The results of the simulation are shown in Figure 7. From the Figure 7 on the left boundary (where compressor is located and the pressure is given) one can see a significant difference in the flux between the two models. On the other hand, on the right boundary the differences between all variables are rather small.

Numerical simulations on a network of pipes
We simulate simple network illustrated on Figure 8. The network has five nodes (enumerated from 1 to 5) and the topology and parameters for pipes and compressors given in Tables 2  and 3.
The initial conditions for the pipe are steady-state. These are not given explicitly, but only specified through the conditions at the nodes. That is, initial nodal conditions are given by pressure at the pressure (slack) node, flow withdrawal at the flow (demand) node,     Table 3: Information on the location of compressors in 5-node pipe network. and compressor ratios. These values are given in Table 4, and the corresponding boundary pressures and flows on each pipe are given in 5.   Table 5: 5-node network initial data by pipe.
The compression ratios, as a functions of time, are given by where T = 86400 sec, and c 1 (0), c 2 (0), and c 3 (0) are given in Table 4. Note that c 2 (t) defines linear interpolation of the points (0, c 2 (0)), (21600, c 2 (0)), (25200, 1.4 · c 2 (0)), (64800, 1.4 · c 2 (0)), (68400, c 2 (0)), and (86400, c 2 (0)). The non-zero flows out of the network at the demand nodes, as functions of time, are given by where T = 86400 sec, and d 3 (0) and d 5 (0) are given in Table 4. Note that c 2 (t) defines linear interpolation of the points (0, d 5 (0)), (12000, d 5 (0)), (15600, 1.2 · d 5 (0)), (48000, 1.2 · d 5 (0)), (51600, d 5 (0)), and (86400, d 5 (0)). We simulate the gas dynamics in the network using three discretizations: (i) staggered discretization (presented in this paper), (ii) implicit lumped elements [31]. For these discretizations the conservation of mass is illustrated on fig. 10 while the descrepancy for the pressure and the flow for some of the nodes is shown on fig. 11. Relative mass error time in hours Staggered discretization Figure 10: Illustration of the discrepancy between the mass of the gas in the pipe and the mass of the gas that left the system through the nodes in the split-step and lumped element (left) and staggered discretization (right). The plot for the staggered discretization shows one data point per minute.
The split-step discretization assumes a constant compressibility factor and works in the pressure-flux formulation. It splits the each step into a linear and non-linear step. The linear step is integrated using characteristics. Compared to the split-step discretization, the explicit staggered discretization has the same efficiency for rather general forms of the compressibility factor and does not require it to be constant.
For the staggered discretization we use a non-ideal gas model and time step t = 1/8 sec. For the lumped element and the split step we use the ideal gas law model with speed of sound given by a = 377.9683 m/s. The time step for the split step discretization is the same t = 1/8 sec, while for the lumped element it is two points per kilometer. Pressure at node 5 (MPa) time in hours lumped element ∆t=4 split-step ∆t=.5 staggered ∆t=.5 staggered ∆t=.05 Figure 11: Illustration of the inflow through the node 1 (top) and pressure at the node 5 (bottom) using ideal gas model and split-step method with 16 points per km, lumped element method with 2 points per km, and staggered discretization with 16 points per km.

Conclusions
We have presented a computationally efficient, second-order accurate, explicit numerical method for simulating gas pipeline dynamics with non-ideal equation of state models. The method is applied to a model for gas transport in a transmission pipeline network. A new staggered grid discretization, with finite differences in space and time, was applied to this model to simulate initial boundary value problems (IBVPs) and make comparisons with other numerical schemes including that of Kiuchi [16], explicit operator splitting [10], and implicit lumped elements [31]. We demonstrated unconditional stability and proved second order accuracy of the approach. Although the proposed explicit staggered finite difference discretization is restricted by a Courant-Levy stability condition, it employs simple time stepping that does not require solution of a large-scale implicit system. Moreover, the data dependence on each time step is highly localized, which facilitates parallelization of computation of IVP solutions for very large-scale (e.g. continental) pipeline networks.
The numerical experiments on a single pipe verify that the method is better than second order accurate, and also showed significant effects arising from incorporating non-uniform temperature and non-ideal gas modeling for simulation of flows on various time and space scales. The comparison of the numerical scheme to an explicit operator splitting method and an implicit lumped-element method for a test network of pipes shows very close results. Our presented staggered-grid method accurately resolves fast transient dynamics on networks of arbitrary structure, and can be extended in the future to more comprehensive physical models of thermal effects and gas composition.
where b 1 = 1 + a 1 · p atm 10 a 2 G (1.8 · T ) a 3 , and b 2 = a 1 10 a 2 G 6894.75729 · (1.8 · T ) a 3 . (58) We may then write and solving for p yields Here the gas constant R g (in J/kg·K) is obtained from the gas gravity G by where M air = 28.9626 kg/k-mol, and R u = 8314.46 J/k-mol·K.

A.2 Non-isothermal approximation
Suppose that temperature is changing in the pipe.
where p is in Pascal, T is in Kelvin, and c 1 = 527588.415238, c 2 = 14.7, and c 3 = 6894.75729. We may then apply equation (55) where Relevant temperatures are 300 • K at inlet and 275 • K at outlet.