Transient gas pipeline flow: analytical examples, numerical simulation and a comparison to the quasi-static approach

The operation of gas pipeline flow with high pressure and small Mach numbers allows to model the flow by a semilinear hyperbolic system of partial differential equations. In this paper we present a number of transient and stationary analytical solutions of this model. They are used to discuss and clarify why a PDE model is necessary to handle certain dynamic situations in the operation of gas transportation networks. We show that adequate numerical discretizations can capture the dynamical behavior sufficiently accurate. We also present examples that show that in certain cases an optimization approach that is based on multi-period optimization of steady states does not lead to approximations that converge to the optimal state.


Introduction
The isothermal Euler equations are a well-known model for gas flow in pipelines. They form a quasilinear system of hyperbolic partial differential equations. For the range of flows that occurs in the operation of gas networks, often a simpler model is sufficient. In this paper we study a semilinear model that is derived from the isothermal Euler equations and generates reasonable states for sufficiently small Mach numbers and sufficiently high gas pressures without abrupt changes. We present some analytical solutions of the semilinear model. The analytical example solutions help to improve the understanding of the model. For example, the semilinear model has some ghost-solutions that are outside of the range where the model is physically 1 3 valid. These are traveling waves that travel with the speed of sound. For all subsonic velocities, there exist exponential traveling waves solutions. Moreover, for certain parameters there also exist exponential profiles that travel with constant acceleration. These new example solutions complement the product solutions already presented in Gugat and Ulbrich (2017).
As a first step in the optimization of gas networks, it makes sense to consider stationary optimization, see Schmidt et al. (2016). However, since supply and demand are functions of time, the dynamic control for gas pipeline systems allows a more efficient operation. In Groß et al. (2018), instationary flows are considered as a sequence of stationary gas flows. In our paper, the stationary flows for the semilinear model are given in Sect. 4. A comparison between transient optimal control and multi-period optimization of steady states is presented in Sect. 5 using a simplified example. The discussion clarifies the differences between optimal states that are obtained as solutions of the dynamic optimal control problem and optimal states that are generated with the steady state approach. If the problem data have only small variability as a function of time, the quasi-static approach is very useful. However, for problem data with large variations the dynamic model has some advantages.
In Grimm et al. (2019) a multilevel model of the European entry-exit gas market is studied where the transmission system operator (TSO) sets maximal technical capacities at the first level. At the fourth level, the TSO minimizes transport cost. Both levels require to take into consideration the transient dynamics of the gas flow in the networks as soon as the consumer demand varies with time.
Questions about the robustness of the gas transmission system including effects on the power grid are studied in Ahumada-Paras et al. (2021). The simulation of gas networks coupled to power grids is studied in Fokken et al. (2019). State and parameter estimation for gas pipeline networks is considered in Sundar and Zlotnik (2019). This paper has the following structure. In Sect. 2 we present the isothermal Euler equations. In Sect. 3 we derive the semilinear model and present some analytical solutions. In Sect. 4 we consider the stationary solutions of the semilinear model. In Sect. 5, we show how in a dynamic optimal control problem non-monotone pressure profiles can occur, whereas for a quasi-static control, the pressure profiles along a horizontal pipe are always monotone.
In Sect. 6 we present an illustration of the errors caused by a finite differences discretization of the partial differential equation. Moreover, we show that optimal steady states generated with the quasi-static approach in general do not converge to the optimal state for the PDE-constrained optimal control problem. This is a contribution to the discussion of the limits of the multi-period steady state approach. The quasi-static model is not suitable if the problem data change rapidly over time.

The isothermal Euler equations
Let an interval [0, L] be given that represents a pipe e of length L. Let D > 0 denote the diameter of the pipe, fric (x) ≥ 0 the space-dependent Lipschitz continuous friction coefficient and (x) ∈ (−∞, ∞) the space-dependent Lipschitz continuous slope. Define s lope (x) = sin( (x) ) and (x) = fric (x) D . Let g denote the gravitational constant. Let > 0 denote the gas density, p > 0 the pressure and q the mass flow rate. We assume that we have an ideal gas that satisfies the state equation where R e s is the gas constant and T e is the temperature. We study the well-known isothermal Euler equations that govern the flow through a single pipe (see for example Banda et al. 2006). In our analysis also the velocity and the speed of sound c that is given by appear. Thus we have For the Mach number M this yields We consider the case of subsonic flow where the absolute value of the velocity of the gas is strictly less than the speed of sound in the gas, that is for the Mach number we have This is the case that is relevant for gas transportation networks where upper bounds for the velocity of the gas are prescribed in the operation of the gas pipelines. In order to state (2) in terms of the dimensionless Mach-number M and the pressure p, we observe that |M| < 1.

3
Transient gas pipeline flow: analytical examples, numerical… Note that the eigenvalues of the system matrix A(M) are −c and c. In particular, they are constant.

Remark 1
The eigenvalues of the quasilinear system (2) are v + c and v − c . This implies that the error in the eigenvalues can only be small as long as |v| (the absolute value of the Mach number M) is sufficiently small.
Since the first equation in (2) has not been modified, also (8) guarantees the conservation of mass. In terms of the physical variables p and q, the model (8) can be stated as the well-known semilinear model (see for example Herty et al. 2009;Gugat et al. 2018) Note that (9) is semilinear since the system matrix is constant whereas (8) is stated in quasilinear form.

Analytical solutions of the semilinear model
Now we present several transient analytical solutions of (9).
Example 1 Assume that c > 0 and For a function e ∈ L ∞ (−∞, ∞) with e < 0 consider the Cauchy problem with (9) and the initial conditions Then for t ≥ 0 the solution of the Cauchy problem is given by the traveling waves solution This can be seen as follows. We have 1 c 2 p t = − � e (x + c t) = −q x (t, x) , hence the first equation in (9) holds. Moreover, we have Since q |q| p = 1 c e (x + c t) , due to assumption (10) we have (10) c 2 − 2 g s lope = 0.
Hence also the second equation of (9) holds. In Example 1, the solution is traveling with the velocity c and an arbitrary wave profile. In the next Example 2, we present a solution that is traveling with a velocity v with |v| < c.
Example 2 Now we present a parametric family of traveling waves solution where the velocity v ∈ (−c, c) is a real parameter. Note however, in contrast to Example 1 and Example 6 below, traveling waves with these subsonic velocities only exist with a special shape given by the exponential function. Let a number 0 > 0 be given and define Consider the Cauchy problem with (9) and the initial conditions Then for t ≥ 0 the solution of the Cauchy problem is given by the traveling waves solution This can be seen as follows. We have hence the first equation in (9) holds. Moreover, we have and using the definition of in (12) we obtain Hence also the second equation of (9) holds.
Example 3 For certain parameters, system (9) also has exponential wave solutions that travel with increasing speed with constant acceleration. This can happen for sufficiently small values of if the gas flows downwards with constant slope. Let = − 2 . Assume that Let real numbers C 2 ≤ 0 and C 3 be given. Define and the time-dependent velocity of the traveling pressure-wave profile Let v(t) = V � (t) . Then the Cauchy problem for (9) with the initial conditions has the solution This can be seen as follows. It is easy to check that the initial conditions for t = 0 hold. For t ≥ 0 we have hence the first equation in (9) holds. Moreover, On the other hand, since v(t) ≤ 0 , we have .

Thus the second equation in (9) holds if
Since 1 2 = − , this equation holds if But this is equivalent to the definition of . Note that the velocity is decreasing and if < 0 or C 2 < 0 it is strictly decreasing and the solution becomes supersonic after finite time. On the other hand, if | | and |C 2 | are sufficiently small, the velocity remains arbitrarily small for an arbitrarily long given finite time interval.
Example 4 There are also solutions with constant pressure where the flow rate is increasing. Assume that Let a constant pressure value p > 0 and a real number ≥ 0 be given. Define the parameters Then A = p c 2 g |s lope | = c 2 2 p A 2 . Consider the Cauchy problem with (9) and the initial conditions Then the solution of the Cauchy problem is the constant pressure p and the increasing flow rate This can be seen as follows. We have hence the first equation in (9) holds. Moreover, c 2 + = −g s lope .
(15) s lope < 0.  (9) holds. For the Mach number we obtain In particular, the Mach number is uniformly bounded and since it is increasing as an upper bound we have the limit Hence if this limit is less than 1, the flow remains subsonic all the time. In particular, for suitable data it remains within the range where the model is valid. This situation occurs if for a moderately sloped pipe the gas flows downward along the pipe. Then without additional compressor action, the flow is driven by gravity and in the limit approaches a constant stationary state. Note that in this solution the flow rate q changes permanently with time. So the behavior is completely different from stationary states, where the value of q is fixed all the time.
Example 5 Now we present a solution that is initially given by the trigonometric tangent function. Let (15) hold, A be as in (16) and as in (17). Let a real number < 0 be given. As long as t satisfies + t ∈ (− ∕2, 0) , that is the constant pressure p with the increasing flow rate solve (9). This can be seen as follows. It is easy to check that the first equation in (9) holds. Moreover we have Thus the second equation in (9) holds. Note that at the time where + t = 0 , the sign of q(t) changes. This means that the direction of the gas flow changes. For + t ≥ 0 , the solution continues with the flow rate given by the hyperbolic tangent function as in Example 4.
Product solutions of (9) are given in Lemma 3 in Gugat and Ulbrich (2017).
Example 6 For certain parameter values, the semilinear model (9) has traveling waves solutions that travel with the speed of sound c. Note that this analytical solution is not within the range of parameters where the model is physically valid. We discuss this example since it shows that nonphysical solutions exist. Unwanted solutions of this type can be avoided by additional constraints on the range of the pressure values and the gas velocities/Mach numbers in the optimal control problems. With appropriate bounds these additional constraints reduce the size of the feasible set in such a way that only states within the physical range of the model are admitted.
Similarly as in (10) assume that For a function e ∈ L ∞ (−∞, ∞) with e > 0 consider the Cauchy problem with (9) and the initial conditions Then for t ≥ 0 the solution of the Cauchy problem is given by the traveling waves solution This is similar as in Example 1, but with a different sign. This can be seen as follows. We have 1 , hence the first equation in (9) holds. Moreover, we have Since q |q| p = 1 c e (x − c t) , due to assumption (18) we have (18) c 2 + 2 g s lope = 0.

Transient gas pipeline flow: analytical examples, numerical…
Hence also the second equation of (9) holds.

The stationary states
In this section, the stationary states for (9) are studied. The stationary solutions of the quasilinear system (2) have been studied in Gugat et al. (2015). We present the stationary states for the semilinear model not only for horizontal pipes, but also for pipes with non-zero slopes. The details are presented in Gugat et al. (2021). For the stationary states, there exists a constant C s such that c |q| = C s . If s lope = 0 (that is for a horizontal pipe), we have (p) 2 2 x = − 1 2 sign(M) C 2 s . Note that if M > 0 , p is decreasing along the pipe and if M < 0 , the gas flows in the opposite direction and p is increasing along the pipe. We obtain the well-known Weymouth equation (see for example Herty et al. 2009, equation (6)). If M > 0 , the equation has real solutions only as long as x is sufficiently small. This implies that the model is well-defined only as long as the length L of the pipe is sufficiently small. Also for the quasilinear system (2) a blow up of the derivative for the stationary states occurs after a finite length, see (Gugat et al. 2015). Note that for the stationary states with q > 0 , the gas velocity is increasing along the pipe and the blow up occurs at the point where the flow becomes sonic, that is when the gas velocity reaches the speed of sound c. Hence this point is outside of the range of validity of the semilinear model. Now we consider the case that s lope < 0 . Similarly as for horizontal pipes, for a stationary state the flow rate q is constant along the pipe. We only obtain a stationary pressure that is constant along the pipe if the source term on the righthand side of (9) vanishes. This is the case for a constant flow rate q > 0 and the constant pressure p > 0 with We also have the non-constant stationary solutions with q > 0 and with a real number Ĉ such that p(0) =p √ 1 +Ĉ . This can be verified by checking that (q, p(x)) satisfies (9).

A comparison between transient optimal control and multi-period optimization of steady states
While for demand functions that change very slowly in time, it makes sense to consider a quasi-static model, since at each moment the system state is very close to a steady state, for demand functions that move more rapidly with time, in general the state cannot be represented by a sequence of stationary states, An overview over the range of gas demand, for example the daily demand range, can be found in Gupta et al. (2019). Note that the domestic demand is strongly weather related, it depends for example on the temperature, see also (Hribar et al. 2019).
As in Gupta et al. (2019), as an approximation of the shape of the demand profiles (in m 3 ∕h ) for 22 hours (from 0:00 to 22:00) we use the sum of three Gaussian kernels Suitable values of (a k , b k , c k ) 3 k=1 are given in Table IV Mishra et al. (2016). Of course the total gas consumption on the time interval can easily be obtained by integrating q desi .
In order to compare the transient optimal control approach and the quasi-static optimization of steady states, let us first derive a PDE that is further simplified. From (9), under the assumptions that the second partial derivatives exist, s lope = 0 and = 0 we obtain Hence we arrive at the wave equation In order to determine a corresponding pressure, we use the first equation in (9), that is p t = −c 2 q x . For x ∈ [0, L] this yields 1 3

Transient gas pipeline flow: analytical examples, numerical…
In order to come to a clear comparison of transient versus multi-period steadystate modeling we consider the problem to satisfy the consumer demand at the end x = L of a single pipe with inflow at the end x = 0 . We consider the time interval [0, T] . Let a desired pressure value p desi > 0 be given. This is the desired pressure value at the point where the gas enters the system. Let the continuous initial flow q 0 , the corresponding time derivative q 1 and the continuous initial pressure p 0 with p 0 (0) = p desi be given. If the gas is initially at rest, we have q 0 = 0 and q 1 = 0 . We state this problem in the form of the optimal control problem The corresponding N-step steady-state model for a time grid For the solution of (SSM(N)) we obtain optimal the values u k = p desi = p(t (N) k , 0) . The corresponding steady pressures along the pipe are determined by the ordinary differential equation for the steady state that can be solved explicitly, see Sect. 4 and also (Herty et al. 2009). Note that for the steady states that are used in (SSM(N)) , the pressure is always decreasing along the pipe. For = 0 , we have p x (t (N) k , x) = 0 and we obtain constant pressure along the pipe for all time steps. Now we look at the optimal state for (OCP(T)) . The flow can be represented as the sum of traveling waves, . Then the initial conditions hold. For t > L c we define Then the boundary condition q(t, L) = q desi (t) is satisfied. For the pressure along the pipe, due to (22) we have For the spatial derivative this yields For the pressure at x = 0 , we obtain Let us assume that p � 0 (x) = 0 , that is initially the pressure is constant. This makes sense since, as stated in Gupta et al. (2019), from 12:00 h to 03:00 h, there is no significant flow inside the pipelines. Thus we can also assume q 0 = 0 and q 1 = 0 . Then the variation of the spatial derivative p x along the pipe depends on the variations of the derivative of q desi with time. Since these can be quite large, in particular close to the peaks, for the dynamic PDE model the pressure can be both in increasing and decreasing along the pipe. This is in contrast to the quasi-static model.
For the sake of illustration, let us look at an example.
So it turns out that for t > − L c , is 4 L c periodic with the value zero on half of the time intervals.

3
Transient gas pipeline flow: analytical examples, numerical… If p � 0 (L) = 0 for t ∈ (2 L c + 4 L c n, 3 L c + 4 L c n) this yields This illustrates the variations of p x for the solution of (OCP(T)).
In contrast to this, for the solution of (SSM(N)) with the quasi-static approach the pressure p is always decreasing along the pipe. Figure 1 shows the pressure waves for the optimal state for (OCP(T)) for p desi = 10 6 Pa , c = 340m∕s and L = 1000m , = 1000 . Note that the pressure profiles are not always monotone along the pipe.
It is important to emphasize that the PDE model is particularly relevant for dynamic situations, for example if a gas power plant is started. For long-term planning problems (say for 12 months, see for example Aras and Aras 2004) the quasistatic approach can be justified.

Comparison between discretized transient optimal control and multi-period optimization of steady states
In Sect. 6.1 we compare numerical results obtained by discretized models with the solutions obtained by the PDE model. In Sect. 6.2 we compare results that are obtained from a quasi-static approach with the solutions obtained by the PDE model and show that refinement of the time grid in the quasi-static approach does not lead to convergence to the optimal state.

Comparison of the solutions of the PDE with the solutions for discretized models
As in Sect. 5 we consider a problem of optimal nodal profile control, where the profile desired by the consumer at the output node is given. We look for boundary controls such that for the generated flow the distance of the values at the output node to the desired profile is as small as possible.
Problems of optimal dynamic control for gas pipeline systems are studied in Osiadacz and Chaczykowski (2016) with the linear heat equation as a model for the system dynamics. In Krug et al. (2020), the time-domain decomposition for a general class of optimal control problems governed by semilinear hyperbolic systems is presented which could also be applied to gas network problems.
Let us start with the transient solution with constant pressure in Example 4 where the hyperbolic tangent appears. We consider a sloped pipe that goes slightly downwards, that is, (15) holds.
Let ≥ 0 be given. Define the desired pressure as p d and the desired outflow rate at the output node x = L as with A as in (16) and as in (17).
The control problem is to find an input pressure p(t, 0) that minimizes subject to for all t ∈ [0, T] and (9). So the problem is to control the input pressure in such a way that the generated output pressure is as close as possible to the given desired value. From Example 4 we know that the optimal value of this problem is zero, that is the desired pressure value can be reached exactly. The analytical solution for the flow q for Example 4 is shown in Fig. 2. Note that the grid in Fig. 2 is not a discretization grid. For our computations we discretize the time horizon [0, T] and the spatial domain [0, L] equidistantly by with K, N ∈ ℕ . We introduce the notation p k,n = p(t k , x n ). Now we fully discretize the equations in (9) by using the implicit Euler method for the time and the explicit and implicit finite difference method for the space: This discretization resembles of a first-order upwind scheme where we consider the location of the boundary data instead of the sign of the characteristic speeds.
To approximate the objective function we simply use the right Riemann sum since p k,0 = p(0, x) = p d is given data. The optimization problem now reads − g s lope p k+1,n+1 c 2 .
We used MATLAB and its function "fmincon" to solve (29) and then plot the results. For the flow Fig. 2 shows the analytical solution q(t, x) = A tanh( + t) and Fig. 3 shows its numerical counterpart. In the case of the pressure the numerical solution shows a small but noticeable difference to the constant analytical solution p(t, x) = p , as can be seen in Fig. 4. This is due to the errors that are made by discretizing the problem. The objective function value of the solution is 0.1165, so note that it is quite small compared to the value of p. 1 3

The quasi-static approach
For the corresponding N-periodic quasi-static we need a time grid with t (N) In the quasi-static approach, the flow rate is constant along the pipe. Hence the optimization problem is to minimize It is easy to see that for the solution we have the pressure values at x = 0 . In the quasi-static model the generated pressure along the pipe is given by the steady state that is determined by the ordinary differential equation ( √ 1 +Ĉ(t) exp 2 g |s lope | c 2 x .

3
Transient gas pipeline flow: analytical examples, numerical… Therefore, for N → ∞ , the error in the approximation of the constant optimal pressure p d will not converge to zero with refinement of the grid, that is for all x ∈ (0, L] we have Hence the quasi-static model does not lead to approximations that converge to the correct optimal state, for which the pressure is constant along the pipe.

Conclusions
We have presented examples that illustrate the dynamic behavior of the gas flow in pipes.
In the examples, instationary analytical solutions are presented. Numerical examples illustrate that a reasonable discretization can approximate the dynamical behavior with a small numerical error. We have shown that in general, the multiperiod steady state approach is unable to capture this dynamic behavior. We have presented an example where the optimal states generated with the multi-period steady state approach do not converge to the optimal transient state. Moreover, we have illustrated that in a dynamic situation with rapidly varying demand, pressure profiles can occur that are non-monotone along the pipe (see in particular Example 7).
Funding Open Access funding enabled and organized by Projekt DEAL. This work was supported by DFG in the Collaborative Research Centre CRC/Transregio 154, Mathematical Modelling, Simulation and Optimization Using the Example of Gas Networks, Project A05, C03, C05.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.