Classical solutions and feedback stabilization for the gas flow in a sequence of pipes

We consider the subcritical flow in gas networks consisting of a 
finite linear sequence of pipes coupled by compressor stations. Such 
networks are important for the transportation of natural gas over 
large distances to ensure sustained gas supply. We analyse the 
system dynamics in terms of Riemann invariants and study stationary 
solutions as well as classical non-stationary solutions for a given 
finite time interval. Furthermore, we construct feedback laws to 
stabilize the system locally around a given stationary state. To do 
so we use a Lyapunov function and prove exponential decay with 
respect to the $L^2$-norm.


Introduction.
Recently, there has been intense research on the system dynamics in gas networks (see [1,2,3,4,5,8,10,11,17,19] and publications of the Pipeline Simulation Interest Group [18]). Several models for the flow in gas networks have been proposed and discussed both from an analytical and a numerical point of view [10,17]. The most important physical phenomenon is the pressure loss along the pipe due to the friction at the inner surface of the pipe.
In this paper our starting point are the isothermal Euler equations which model the gas flow by a 2 × 2-PDE-system of conservation laws for each pipe. Motivated by the fact that shocks occurring for non-classical solutions can damage the pipes we study the existence and uniqueness of classical stationary and non-stationary solutions for a linear sequence of N gas pipes coupled by N − 1 compressor stations. Furthermore, we develop linear feedback laws to stabilize this system locally around a stationary state.
Our work is based upon [8] where the feedback stabilization is analysed for a system of two pipes coupled by one compressor. However, the approach of [8] does not work for a network of more than two pipes that are serially connected by compressors. In this paper we develop new feedback laws that allow feedback stabilization for a network consisting of N pipes coupled by N − 1 compressor stations (see also Remark 3.4). Related problems of optimal control of networked hyperbolic systems have been considered in [7,20].
2. Analysis of the gas flow in a single pipe. In this first section we analyse the gas flow in a single pipe. We study the isothermal Euler equations, their transformation to Riemann invariants, stationary and non-stationary classical solutions and a Lyapunov function for this system of one pipe. Having analysed these aspects for a single pipe, we can easily transfer them to the network case in section 3.
2.1. Isothermal Euler equations. As the lengths of the pipes in real gas networks are much larger than their diameters, we can use a one-dimensional model for a gas pipe and parametrise its length by the space interval [0, L] (x-coordinate). We see the pipe as a graph consisting of one edge and two boundary nodes. One boundary node is at the end x = 0 of the pipe and the other one at the end x = L. Furthermore, we consider the system on a finite time interval [0, T ] (t-coordinate) with T > 0.
The state of the gas in a pipe can be characterised by the density of the gas ρ = ρ(t, x) and the mass flux in the pipe q = q(t, x). The density ρ is always positive whereas the sign of q depends on the direction of the mass flux. A common model to describe the relationship between these two variables are the isothermal Euler equations [1,2,8,10,11,15]. For one pipe they have the form The first equation states the conservation of mass and the second one is the momentum equation. The right hand side of the momentum equation is a friction term with the friction factor f g > 0 and the diameter of the pipe δ > 0. Note that the friction force always acts in the opposite direction as the direction of the mass flux q. The constant a > 0 is the sonic velocity in the gas. It depends on the gas under consideration and the gas temperature according to the equation where Z is the natural gas compressibility factor, R the universal gas constant, T the absolute gas temperature and M g is the molar weight of the gas. The relationship between the density ρ and the gas pressure p is described by the pressure law p = a 2 ρ.
In this paper we consider subsonic or subcritical states, i.e. states with |q|/ρ < a. If |q|/ρ becomes equal to a, the state becomes critical. Equations (1) and (2) together form the p-system which is a hyperbolic balance law. In the following we will always assume that the mass flux is positive, i.e. q > 0. Using the chain rule we get that the equations (1) and (2) together can equivalently be written as The matrixÂ is not diagonal and, thus, the derivatives ∂ x ρ and ∂ x q are coupled in equation (4). By transforming the system to Riemann invariants in section 2.2 we get a diagonal system matrix.

2.2.
Transformation to Riemann invariants. In this section we transform the system described in section 2.1 to Riemann invariants. For subsonic states, the system matrixÂ(ρ, q) has the strictly positive eigenvalue λ + = λ + (ρ, q) = q ρ +a > 0 and the strictly negative eigenvalue −λ − = −λ − (ρ, q) = q ρ − a < 0. Note that λ − denotes the absolute value of the strictly negative eigenvalue. The corresponding left eigenvectors are Let L = L(ρ, q) := l + l − denote the matrix of the left eigenvectors. Then det(L) = 2 a ρ 2 = 0 and We obtain the Riemann invariants as Their partial derivatives with respect to ρ and q satisfy We can restate ρ, q and λ ± in terms of R ± : As we assume that the mass flux q is positive, we have R + + R − < 0 and a state is subsonic (i.e. q/ρ < a) if and only if We define the set S as Multiplying equation (4) by the matrix L and using equation (8) combined with the chain rule we can transform the system equation (4) to the following system equation in terms of the Riemann invariants R + and R − : with the diagonal system matrix and the source term Equation (7) expresses R + and R − in terms of ρ and q. Equations (10) and (11) give ρ and q in terms of R + and R − . Now we consider the question how R + can be expressed by R − and q. To do so, we observe the function If the state (R + , R − ) is subsonic (see (12)), a simple calculation shows that Hence, due to the monotonicity property of the function Q, there exists a C 1function F + such that for all fixed real values R − and q with the equation holds and the state (F + (q, R − ), R − ) is subsonic. Analogously, also R − can be expressed by R + and q.

Stationary states.
In this section we analyse the existence of stationary subsonic continuously differentiable states. We denote the Riemann invariants for the stationary case which do not depend on the time variable t asR + =R + (x) and R − =R − (x). Equation (14) yields the following ordinary differential equation for R + andR − : As we want to consider subsonic stationary states with a positive mass flux, we have (R + (x),R − (x)) ∈ S (see (12), (13)). Hence, we can analyse the autonomous system d dx (20) for x ≥ 0 with corresponding boundary conditions (R + (0),R − (0)) = (R +,0 ,R −,0 ) ∈ S, which we prescribe at x = 0. For the existence and uniqueness of a C 1 -solution of (20) we obtain the following theorem: Theorem 2.1 (Existence and Uniqueness of Stationary C 1 -Solutions). For (R +,0 ,R −,0 ) ∈ S define the positive real number Then there exists a unique continuously differentiable solution Proof. For (R + (x),R − (x)) ∈ S the right hand side of (20) is continuous and continuously partially differentiable with respect toR + andR − and, thus, is locally Lipschitz-continuous. Hence, due to the Picard-Lindelöf theorem there exists a unique continuously differentiable solution (R + (x),R − (x)) of (20) that satisfies (R + (0),R − (0)) = (R +,0 ,R −,0 ). Now we prove that lim and lim which implies that the maximal interval of existence of (R + (x),R − (x)) as a C 1 - and, hence,R + (x) +R − (x) is strictly monotonically decreasing with respect to x. For −2a < z < 0 consider the function Φ(z) := 4a 2 z 2 + 2 ln(−z). The derivative Φ (z) can be calculated as Hence, Φ is strictly monotonically increasing on (−2a, 0). Using the chain rule equations (24) and (25) yield that Due to the continuity of Φ and due to the definition of x 0 we have Hence, the monotonicity property stated by equation (25) implies (22) and, thus, equation (20) yields (23).
Finally, the critical length x 0 can be written as and, hence, the monotonicity property of Φ (see (25)) and the inequalityR +,0 + R −,0 > −2a imply that x 0 is actually positive.
Remark 2.2. Theorem 2.1 says that the solution (R + (x),R − (x)) cannot be extended beyond x 0 , neither as a subsonic solution nor as a continuously differentiable solution. Hence, Theorem 2.1 implies that a stationary subsonic C 1 -solution exists along the whole pipe if and only if the length L of the pipe is less than the critical length x 0 , i.e. L < x 0 .
Having analysed stationary subsonic continuously differentiable states in terms of Riemann invariants we will now transfer our results to the densityρ =ρ(x) and the mass fluxq =q(x) for the stationary subsonic gas flow.
Equation (1) directly implies that in this case the mass flux is constant with (see (11) For the derivative d dxρ (x) equations (10) and (20) and the chain rule imply and, thus, in the stationary subsonic case the density is strictly monotonically decreasing on the interval [0, x 0 ), and for x → x 0 − a blow up in the derivative d dxρ (x) occurs. Hence, due to the pressure law (3) also the pressure p decreases along the pipe.
Furthermore, equations (10) and (11) yield that Using equation (7) one can easily show that the critical length x 0 can be expressed in terms ofρ(0) andq as For f g → 0 we have d dxρ (x) → 0 and x 0 → ∞, that is, in the frictionless case, the subsonic stationary solutions are exactly all constantsρ > 0,q > 0 withq/ρ < a and there is no limit for the length of the pipe.

Remark 2.3.
To analyse the stationary subsonic case we have considered the system equation (14) in terms of the Riemann invariantsR + andR − . However, one can also consider the original isothermal Euler equations (1) and (2) for the stationary case in terms ofρ andq. This is elaborately done in [8] and leads to the same results.

2.4.
Local existence of non-stationary semi-global solutions. In this section we analyse the existence of non-stationary solutions R ± (t, x) =R ± (x) + r ± (t, x) locally around a stationary stateR ± (x) on a given finite time interval [0, T ]. For these solutions we prescribe initial conditions R ± (0, x) =R ± (x) + r ± (0, x) in a sufficiently small C 1 -neighbourhood of the given stationary stateR ± (x) (i.e. ||r ± (0, x)|| C 1 ([0, L]) is sufficiently small). The solutions are called semi-global as they exist on the finite time interval [0, T ].
Thereby, we refer to a result of T. Li, B. Rao and Z. Wang [13] about the existence of semi-global C 1 -solutions for initial-boundary value problems with nonlocal boundary conditions for quasilinear hyperbolic systems and apply this result to our system in terms of Riemann invariants (see also [21]).
Let a stationary stateR = (R + ,R − ) be given. Consider a non-stationary state Then equation (14) and ∂ tR = 0 yield: Using equation (14) for the stationary case, we obtain the following equation for r + and r − : with the diagonal system matrix and the source term (29) where the notationR = (R + ,R − ), r = (r + , r − ) and ±λ ± = −R + +R− 2 ± a (see (9)) has been used. Furthermore, the equations (10) and (11) imply for the density ρ and the mass flux q in the non-stationary case with the densityρ and the mass fluxq belonging to the stationary state.
The following lemma states the existence of a C 1 -solution of (27) on a finite time interval for appropriate initial and boundary conditions. Since we will analyse networks of several pipes in section 3, in Lemma 2.4 we already consider a system of N (N ≥ 1) equations of the form (27). Thereby, we denote variables referring to equation i with a superscript (i) (i ∈ {1, ..., N }). The lengths L (i) (i ∈ {1, ..., N }) of the pipes may be different, but all pipes have the same diameter δ and the same friction factor f g . The description of a network of N pipes that are serially connected by N − 1 compressors can be found in section 3.1.
.., N }) be given. Then there exists ε > 0 (depending on T ) such that the following statement holds: Consider the system of 2N partial differential equations for r (i) Then, for all i ∈ {1, ..., N }, there exists a unique solution (r and the C 1 -compatibility conditions at the points (t, x) = (0, 0) and (t, Proof. Without restriction, we can assume that all pipes have the same length 8,12]. As a consequence of this change, the equations (32) become The equations (37) have the same form as equation (1.4) from [13] where in our case, using the notation of [13], the vectors l i are the canonical unit vectors. For r (i) .., N }) the right hand sides of (37) are equal to zero and, − < 0) and, thus, the inequality (1.6) from [13] holds (i ∈ {1, ..., N }).
Hence, all requirements of Theorem 2.1 from [13] are satisfied and Theorem 2.1 yields our assertion.
Remark 2.5. In Lemma 2.4 the usual C 1 -compatibility conditions ensure that the initial and the boundary data and their derivatives fit together at the points (t, x) = (0, 0) and (t, x) = (0, L (i) ). The C 1 -compatibility conditions follow from a straight forward calculation for ϕ (i) , g In [13], quasilinear hyperbolic systems with nonlocal boundary conditions are considered, where the eigenvalues and the source term of the system under consideration depend on the unknown variables but not explicitly on x and t. In [21], for the case of local boundary conditions, quasilinear hyperbolic systems with also space-and time-dependent eigenvalues and source terms are studied. The dependence of x and t can also be transferred to the case of nonlocal boundary conditions as we have used it in the proof of Lemma 2.4.

Lyapunov function.
In the previous section we have considered the existence of a non-stationary solution r + (t, x) and r − (t, x) locally around a stationary state. In this section we introduce a Lyapunov function that can be seen as a weighted L 2norm of r + and r − with exponential weights that are defined using the eigenvalues λ + and −λ − for the stationary state (see (41)). We adopt this Lyapunov function and the estimate for its time derivative from [8]. In section 3.3 we will use this Lyapunov function to develop feedback laws for the stabilization of our system locally around a stationary state.
Assume that we have a given stationary stateR ± (x) with corresponding eigenvalues ±λ ± = −R + +R− 2 ± a (see (9)). In the previous section we obtained the partial differential equation (27) for a non-stationary solutionR ± (x) + r ± (t, x) locally around the stationary stateR ± . Equation (27) can also be written in the following form: where the notatioñ has been used. Note thatK + andK − depend on t and x. Consider the Lyapunov function with A > 0, B > 0 and for µ > 0 with exponential weights From [8] we obtain the following lemma for the estimate of the time derivative E (t):  (39) and (40). Define the positive real number Furthermore, define Assume that we have Then we can choose A > 0 and B > 0 such that both hold or such that both and hold.

GAS FLOW IN A SEQUENCE OF PIPES 701
The proof of this lemma uses integration by parts and the monotonicity property of h + (x)/h − (x). In this paper we do not want to repeat the complete proof of this lemma (It can be found in [8].), but we just would like to give the following remark on Lemma 2.6: In section 3.3 we will develop feedback laws which guarantee that the boundary terms in equation (54) are non-positive and, hence, 3. Gas flow through a sequence of pipes coupled by compressors. In this section we consider the gas flow through a linear sequence of N pipes coupled by N − 1 compressor stations. In section 2.3 we have seen that in the stationary case the densityρ decreases along the pipe and for x → x 0 − (with the critical length x 0 ) we haveρ →q a , i.e. the state becomes critical, and a blow up in the derivative ofρ occurs. Therefore, to increase the maximum distance of gas transportation, one can use a linear sequence of several pipes with sufficiently short single lengths and couple subsequent pipes by a compressor station. The compressor stations increase the density again which has decreased along the preceding pipe.
3.1. System equations and coupling conditions at the compressors. Assume that we have a linear sequence of N (N > 1) pipes coupled by N − 1 compressors. We number the pipes consecutively from pipe 1 to pipe N as they are arranged in the sequence. In the following we will always denote variables referring to pipe i with a superscript (i) (i ∈ {1, ..., N }).
The compressor coupling pipe i and i + 1 is denoted as compressor i, i + 1 and the values referring to the compressor i, i + 1 are denoted with a superscript (i, i + 1) The pipes are parametrised by space intervals [0, L (i) ] such that at the compressor i, i + 1 the ends x = L (i) of pipe i and x = 0 of pipe i + 1 meet. At the end x = 0 of pipe 1 there is the gas producer and at the end x = L (N ) of pipe N the consumer. Hence, we have a positive mass flux along the whole sequence from the boundary point x = 0 of pipe 1 to the boundary point x = L (N ) of pipe N . Furthermore, all pipes have the same diameter δ and the same friction factor f g (see Figure 1 for N = 3). According to our analysis for a single pipe in section 2.1 now for each pipe i we have the system equation in terms of ρ (i) and q (i) (see (4)): with the system matrixÂ(ρ (i) , q (i) ) and the source termĜ(ρ (i) , q (i) ) as in (5) (14)): (15) and the source term At the compressor i, i + 1 (i ∈ {1, ..., N − 1}) we have the following coupling conditions [1,8,11]: for t ∈ [0, T ] with the parameter κ ∈ [ 2 7 , 2 5 ] which depends on the gas under consideration. The first equation says that the mass flux does not change while passing the compressor. In the second equation u (i,i+1) (t) is the compressor power which is always non-negative. For u (i,i+1) (t) > 0 we have ρ (i+1) (t, 0) > ρ (i) (t, L (i) ) and, thus, the compressor increases the outlet density as well as the outlet pressure (see (3)). The compressor i, i + 1 is switched off if u (i,i+1) (t) = 0.
For the time derivative E S (t) we obtain the following theorem:    (46), (47) with L, µ,λ ± ,K ± replaced by L (i) , µ (i) , λ Then, for c := 1 2 min{µ (1) , ..., µ (N ) }, the following inequality holds (t ∈ [0, T ]):  − (·, L (i) ) ≡ 0, satisfy this feedback law. To keep this feedback law, the compressor i, i + 1 has to maintain the mass flux q (i) (t, L (i) ) as in (77), which is close to the mass fluxq for the stationary case if |r Remark 3.4. In [8] the feedback stabilization of a network consisting of two pipes coupled by one compressor is considered. Thereby, linear feedback laws are active at the producer's and the consumer's end of the network and the compressor only has to maintain the constant mass fluxq. Unfortunately, this approach does not work for a network of more than two pipes coupled by compressors as then there exist pipes with compressors at both ends. For such a more complex network also the compressors have to maintain more complex feedback laws. The technical reason is that, due to the coupling condition (58), the relationship between r Together with (85) this yields for i ∈ {1, ..., N }: and .
Hence, Lemma 2.6 yields the estimate In the following we show that the choice of A (i) , B (i) and k (i) in (75), (76), (79), (80), (83) (i ∈ {1, ..., N }) guarantees that the boundary terms in (87) are non-positive, which implies that the inequality (86) holds. First, for the boundary terms at the end x = 0 of pipe 1, we prove the inequality (t ∈ [0, T ]) A (1) (r For ε 1 sufficiently small we have (see (17)) q ∈ (0, a exp(1 +R  and, thus, we can use the function F + from (18) and get R − (0)) and r − (0) + r 4. Example. In section 2.3 we have seen that due to the friction force there is a loss of pressure along the pipe and that a subsonic stationary C 1 -solution exists along the pipe if its length is shorter than the critical length (see (26)) x 0 = δ f g a 2ρ (0) 2 q 2 + 2 ln q ρ(0) − 1 − 2 ln(a) .
The critical length can be increased by reducing the friction at the inner surface of the pipes or by increasing the inflow pressure. The offshore gas pipeline which is currently under construction in the Baltic Sea, for example, has a special interior antifriction coating and will operate with pressures around 200bar [14]. Thus, this pipeline will transport natural gas over a distance of 1220km from Portovaya Bay, Russia, to Greifswald, Germany, without intermediate compressor control. 5. Summary and outlook. In this paper we have analysed the gas flow through a network of serially connected pipes modelled by the isothermal Euler equations. We have considered the existence and uniqueness of subsonic stationary C 1 -solutions and the semi-global existence of non-stationary solutions locally around a stationary state. We have also developed boundary feedback laws to stabilize the system around a stationary state. To do so we have used a Lyapunov function which yields exponential decay of the L 2 -norm of the difference between the stationary state and a non-stationary solution locally around it.
We have presented our results for a linear sequence of pipes coupled by compressor stations. However, our results can be extended to the flow in tree-like networks by considering star-shaped junctions of pipes. Related problems of star-shaped string networks, that are governed by wave equations, have been studied in [9]. Furthermore, our Lyapunov function can be modified in order to obtain exponential decay of the H 2 -norm of the difference between the stationary state and the non-stationary solution (as in [6]) and to stabilize the system on an infinite time interval [0, ∞). Another open problem is the question whether the stabilization still works if one of the compressors breaks down. We will consider these aspects in a forthcoming paper.