An inexact Petrov-Galerkin approximation for gas transport in pipeline networks

This paper studies the discretization of gas transport in pipeline networks by an inexact Petrov-Galerkin method. A full convergence analysis is presented for single pipes under the assumption of a linear friction law and the possible extension to pipe networks is discussed. The generalization to nonlinear gas transport models and the efficient implementation by hybridization is investigated numerically.


Introduction
The flow of gas in a horizontal pipeline of constant cross section is described by [2] A∂ t ρ + ∂ x m = 0 (1) Here A and D are the cross section and diameter of the pipe, and λ is a dimensionless friction parameter. The functions ρ, p, and m describe the density, pressure, and mass flow rate of the gas. Under isothermal flow conditions, one has with constant c denoting the speed of sound. In practically relevant scaling regimes, the nonlinear term on the left hand side of (2)  justified by an asymptotic analysis [2,7]. Using this simplification and equation (3) to eliminate the density, one arrives at evolution problems of the general form where a and b are positive constants and d = d(p, m) denotes a state dependent friction coefficient. For our analyis, we will consider d = d(x) as a function depending only on space which can be justified, e.g., by linearization around a steady state.
Corresponding models for the gas flow on pipe networks are obtained by coupling the flow equations for single pipes via algebraic conditions [9,10]; see below.
The discretization of (4)- (5) and its extension to pipeline networks has been discussed intensively in the literature. In [9], a Galerkin approximation for (1)-(2) with cubic Hermite polynomials is investigated numerically. The discretization of transient gas flow models is also studied [2,8]. An entropy stable finite volume method is proposed in [10], and an energy stable mixed finite element approximation is investigated in [3]. Apart from [9], all methods discussed above are of lowest order and no rigorous convergence analysis is given.
In this paper, we study the discretization of (4)-(5) by a Petrov-Galerkin approach of potentially high order. The resulting scheme is shown to be stable which allows us to prove order optimal convergence rates. By using an appropriate functional analytic setting, the convergence results can be generalized almost verbatim to pipeline networks. A hybridization strategy will be discussed that facilitates the implementation and that allows to incorporate non-standard coupling conditions. The proposed method formally also allow to treat nonlinear models of gas transport and, in principle, high order convergence can be obtained in practically relevant regimes.

Notation and preliminaries
Let x L < x R and denote by L p (x L , x R ) and W k, p (x L , x R ), k ≥ 0 the standard Lebesgue and Sobolev spaces. The scalar product and norm of L 2 (x L , x R ) are written as (v, w) and v = v L 2 . Other norms will be designated by subscripts. We write H k (x L , x R ) = W k,2 (x L , x R ) for the Hilbert spaces and define for convenience. By L p (0, T; X) and W k, p (0, T; X) we denote the Bochner spaces of The value of f (t) may then itself be a function. In the following, we consider the linear system of differential equations for x L < x < x R and t > 0 with homogeneous boundary conditions Inhomogeneous and more general boundary conditions can be considered as well and our analysis applies with minor modifications. We will assume that (A1) a, b are positive constants, and For given f , g ∈ L 2 (0, T; L 2 (x L , x R )) and initial values p(0) ∈ H 1 0 , m(0) ∈ H(div), existence of a unique solution follows from semigroup theory. Any smooth solution of problem (6)-(8) also satisfies p(t) ∈ H 1 0 , m(t) ∈ H(div), and for all v, q ∈ L 2 (x L , x R ) and all 0 < t < T. This variational characterization will be the starting point for our discretization approach introduced in the next section.

Petrov-Galerkin approximation
We call T h := {T n : 1 ≤ n ≤ N } the mesh and denote by h n = |x n − x n−1 | and h = max n h n the local and global mesh size, respectively. Let be the space of piecewise polynomials on the mesh T h . We fix k ≥ 1 and search for approximations for the solutions p(t), m(t) of problem (6)- (8) in the spaces of continuous piecewise polynomials with appropriate boundary conditions. As finite dimensional test spaces for the variational problem (9)-(10), we choose consisting of discontinuous piecewise polynomials of lower order k − 1. We denote by I k h : and and let π k−1 h : Note that both projection operators I k h and π k−1 h can be defined locally on every element. Moreover, they are mutually related other by the commuting diagram property For the approximation of problem (6)-(8), we then use the following approximation.

Problem 1. (Inexact Petrov-Galerkin method) Find functions
, and such that The well-posedness of this problem follows from the results of the next section.

Discrete stability estimates
We now derive some discrete stability estimates that yield well-posedness of the semidiscrete method and that allow us to establish error estimates of optimal order.
with constant C(T) ≤ CT and C independent of T and the solution.
Proof. Let us start by noting that By identity (17), integration-by-parts, and the boundary conditions (8), one can verify that and Young inequalities, and using positivity of d, we then obtain the estimate d dt The Gronwall lemma and the choice α = 1/T finally yields the assertion.
Note that the above estimate does not yet give full control over the solution. A repeated application, however, allows us to prove the following stability estimate.

Lemma 2.
Let p h , m h denote a solution of Problem 1. Then for all 0 ≤ t ≤ T with C (T) = C T and C independent of T and of the solution.
Proof. As a direct consequence of the Poincaré inequality, one has The first terms in these estimates are already covered by Lemma 1. From the two equations (18)-(19) with q h = ∂ x m h (t) and v h = ∂ x p h (t), we further deduce that Bounds for π k−1 h ∂ t p h (t) and π k−1 h ∂ t m h (t) can be obtained by formally differentiating (18)-(19) with respect to time and applying Lemma 1 for the resulting system. A combination of the above estimates then yields the assertion of the lemma. Remark 1. Let us note that Problem 1 formally amounts to a finite dimensional system of differential algebraic equations. From the stability estimates of Lemma 2 and [6, Theorem 4.12], one can deduce that this system is solvable for any choice of admissible initial values. The semidiscretization is thus well-defined.

Error estimates
We start by decomposing the error via p − p h ≤ p − I k h p + I k h p − p h and m − m h ≤ m − I k h m + I k h m − m h into approximation and discrete error components. For the first part, we can utilize the following well-known estimates. For Here By the a-priori estimates of Lemma 2, one then obtains the following result.
with a constant C (T) = C T and C independent of T and of the solution.
Proof. We apply Lemma 2 for p h (t) = I k h p(t) − p h (t) and m h (t) = I k h m(t) − m h (t) and then estimate the terms on the right hand side of the result step by step. By definition of the initial values, we have p h (0) = m h (0) = 0. Moreover, where we used the definition of the initial value m h (0) in the second and (17) in the third step. Thus π k−1 h ∂ t p h (0) ≤ I k h ∂ t p(0) − ∂ t p(0) , and in a similar manner, . This explains the first two terms in the estimate in the lemma. The terms under the integral are derived by estimating π k−1 h f (t) , π k−1 h g(t) and the derivatives π k−1 h ∂ t f (t) , π k−1 h ∂ t g(t) via the triangle inequality, and noting that where we used that d is piecewise constant.

Remark 2.
A similar result can be proven for piecewise smooth d ∈ W 1,∞ (T h ) and additional terms of the form , the product of the two terms again has optimal approximation order.
By combination of the above estimates, we finally obtain the following result. Furthermore, let (p, m) be a sufficiently smooth solution of (6)- (8). Then for all 0 ≤ t ≤ T, one has For sufficiently smooth solutions, the proposed method thus yields convergence with the optimal order that can be expected.

Extension to networks
We now illustrate that our method and the convergence results of the previous section can be generalized easily to pipe networks. Let (V, E) denote a directed graph with vertices v ∈ V and edges e ∈ E; see Figure 1 for illustration. For any To every edge e, we associate a positive length e , and we identify e with [0, e ] in the sequel. This allows us to define spaces L p (E) = {v : v| e ∈ L p (e)} and H 1 (E) = {v ∈ L p (E) : v| e ∈ H 1 (e)} of, respectively, integrable and piecewise smooth functions on the graph. The flow of gas in a pipe network is then described as follows: On every edge e representing a pipe, we require that where f e = f | e denotes the restriction of a function f ∈ L p (E) to one edge. The equations for the individual pipes are coupled by algebraic conditions at the pipe junctions, and at the boundary vertices, we assume that Inhomogeneous and other types of boundary conditions can again be incorporated with minor modifications. For the analysis of the problem, we now utilize the spaces for all q ∈ L 2 (E), v ∈ L 2 (E), and all 0 < t < T. Here (v, w) = e (v e , w e ) e with (v e , w e ) e = ∫ e v e w e dx denotes the scalar product on L 2 (E).

Remark 3.
Let us note that (29)-(30) has exactly the same form as the variational problem (18)-(19) on a single pipe. The inexact Petrov-Galerkin method and all results derived in the previous sections therefore translate almost verbatim to the network setting; let us refer to [4] for details and similar results for a different method, and to Section 9 for numerical illustration.

Remarks on the efficient implementation
In the discretization of (29)-(30), also compare with (18)-(19), the continuity and boundary conditions (24)-(26) are directly incorporated in the definition of the spaces Q h ⊂ H 1 0 and V h ⊂ H(div). For the implementation, it may be more convenient to use larger spaces Q h , V h ⊂ H 1 (E), and to enforce some of the boundary and coupling conditions (24)-(26) explicitly by additional equations. Using the wording of [1], this approach of relaxing continuity conditions might be called hybridization. Since the resulting method is algebraically equivalent to the original scheme based on function spaces with incorporated coupling and boundary conditions, all results of the previous sections apply verbatim also to the method obtained after hybridization.

Nonlinear problems
The formal extension of the Petrov-Galerkin method to nonlinear problems is straight-forward. The discrete variational formulation for (1)-(2), for instance, reads Numerical quadrature can be used in practice to facilitate the handling of the nonlinear terms. We do not give a complete convergence analysis here, but instead, we will demonstrate by numerical tests that for smooth solutions, the convergence results of Theorem 1 remain valid, at least in the practically relevant case of nonlinear friction.

Numerical results
We now illustrate the theoretical results of Section 5 by numerical tests. For our computations, we consider the pipe network depicted in Figure 1. As a first test case, we consider the linear problem (22)-(25) with inhomogeneous boundary conditions and we set p v 1 (t) = 1 and p v 6 (t) = 1 + 1 2 sin(πt) in the following. All pipes are chosen of unit length = 1 and the model parameters are set to a ≡ b ≡ d ≡ 1. The simulation is started from a stationary state for the boundary values at initial time. The results of the computations are summarized in the left column of Table 1. As predicted by our theoretical results, we observe second order convergence.
We now repeat our numerical tests for the same network but with a semilinear gas flow model resulting from (1)-(3) by dropping the nonlinear term ∂ x ( m 2 Aρ ) in equation (2). The model parameters are chosen as A = 1, c = 1, and λ/(2D) = 7/2; the latter was selected such that average of the resulting mass flow was similar to that of the linear model considered above. The computational results are depicted in the middle column of Table 1. Also for this nonlinear friction model, we observe second order convergence. These results can be explained theoretically in a similar way as those for the linear case by using a perturbation argument. In the right column of Table 1 Errors e h = (a p h (T ) − p h/2 (t) 2 + b m h (T ) − m h/2 (T ) 2 ) 1/2 at time T = 10 obtained with the Petrov-Galerkin approximation for the network problem with different gas flow models: linear model (left), semilinear model (middle), and quasilinear model (right  Table 1, we display the corresponding results for the quasilinear flow model (1)-(3) with the same parameters as used in the semilinear case. Note that a decrease in the convergence rates to first order is observed here. This is no surprise, since our analysis heavily relied on the anti-symmetry of the spatial derivative terms in (18)-(19), which is no longer valid for the quasilinear model (1)-(2). In Figure 2, we display the flow rates m| v at the boundary vertices v 1 and v 6 for the three different gass flow models discussed above as function function of time. The results are in reasonable agreement. In summary, the semilinear model seems to yield the best compromise between modelling errors and convergence order.