BOUNDARY FEEDBACK STABILIZATION OF A SEMILINEAR MODEL FOR THE FLOW IN STAR-SHAPED GAS NETWORKS∗

The flow of gas through a pipeline network can be modelled by a coupled system of 1-d quasilinear hyperbolic equations. In this system, the influence of certain source terms that model friction effects is essential. Often for the solution of control problems it is convenient to replace the quasilinear model by a simpler semilinear model. In this paper, we analyze the behavior of such a semilinear model on a star-shaped network. The model is derived from the diagonal form of the quasilinear model by replacing the eigenvalues by the sound speed multiplied by 1 or −1 respectively. Thus in the corresponding eigenvalues the influence of the gas velocity is neglected, which is justified in the applications since it is much smaller than the sound speed in the gas. For a star-shaped network of horizontal pipes for suitable coupling conditions we present boundary feedback laws that stabilize the system state exponentially fast to a position of rest for sufficiently small initial data. We show the exponential decay of the L-norm for arbitrarily long pipes. This is remarkable since in general even for linear systems, for certain source terms the system can become exponentially unstable if the space interval is too long. Our proofs are based upon an observability inequality and suitably chosen Lyapunov functions. At the end of the paper, numerical examples are presented that include a comparison of the semilinear model and the quasilinear system. Mathematics Subject Classification. 93D15, 35L60. Received May 27, 2020. Accepted June 1, 2021.


Introduction
The flow of gas through pipelines is governed by a quasilinear system of balance laws (see for example [2]). We consider a model for a gas pipeline network where at the nodes the solutions for the adjacent pipes are coupled by algebraic node conditions that require the conservation of mass and the continuity of the pressure. The eigenvalues of the system have the form c + v and −c + v, where v denotes the velocity of the gas and c denotes the sound speed. In the operation of gas transportation systems the velocity of the gas flow is much smaller than the sound speed. Therefore, in order to obtain a semilinear model the eigenvalues are replaced by this sound speed, that is by −c and c. It is important to stress that the source term plays an essential role in the model of gas network flow, since if the gas is not at rest, the friction effects lead to a decrease of the pressure along the pipe in the direction of the flow.
In this paper, consider the system with absorbing Riemann boundary conditions of Dirichlet type. We show that on a given finite time horizon [0, T ], the gas flow in a star-shaped network can be steered to a position of rest exponentially fast in the sense of the L 2 -norm if the L ∞ -norm of the initial state is sufficiently small. Moreover, we also show that on an infinite time interval the H 1 -norm decays exponentially fast if the H 1 -norm of the initial state is sufficiently small. From the point of view of applications it is desirable to know that is is possible to steer the flow close to a standstill since sometimes flow reversal is required, for example in the network in Belgium that is used for gas transit in different directions. For a smooth transition this requires to bring the gas almost a standstill before driving the flow in reverse direction. From the point of view of mathematical control theory this is important since while many stabilization results for systems of conservation laws have been studied before, for systems of balance laws with linear source terms exponential stabilization is in general not possible (see [9,14], Thm. 2).
The exponential stabilization of the gas flow governed by the isothermal Euler equations in fan-shaped networks in the L 2 -sense has been studied in [12]. For a single pipe, a strict H 1 -Lyapunov function and feedback stabilization for the quasilinear isothermal Euler equations with friction have been studied in [8]. These results about the quasilinear system have required restrictive upper bounds on the lengths of the pipes. Similar results have been obtained for H 2 -Lyapunov functions, that allow an extension to an infinite time horizon (see [13]).
The finite time stabilization of a network of strings is studied in [1,15]. Finite-time control for linear evolution equation in Hilbert space has been studied in [21]. The finite-time stabilization of hyperbolic systems with zero source term over a bounded interval and the exponential decay for small source terms has been analyzed in [20]. The stabilization of the wave equation on 1-d networks has been considered in [23]. The exponential stability of a linear model for the propagation of pressure waves in a network of pipes has been studied in [7]. The limits of stabilizability in networks of strings are studied in [9]. This paper has the following structure. In Section 2 we introduce the quasilinear isothermal Euler equations and node conditions that govern the flow through a gas pipeline network. In Section 3, we present the corresponding Riemann invariants and transform the system in diagonal form. Then we derive the semilinear model that provides an approximation for small gas velocities.
In Section 4, a well-posedness result is presented. We are working in the framework of solutions that are defined through a fixed point iteration along the characteristics. The definition of the fixed point mapping is derived from the integral equations along the characteristic curves that are known a priori for the semilinear model.
We show that with an absorbing Riemann Dirichlet boundary feedback, the flow through a star-shaped network of pipelines is driven to zero in H 1 exponentially fast. For the proof we first introduce a quadratic L 2 -Lyapunov function and show that it decays exponentially fast on finite time intervals without additional constraints on the lengths of the pipes. In the proof we use an observability inequality for the L 2 -norm. Then we define a quadratic Lyapunov function with exponential weights to show that the time derivatives also decay exponentially fast. This yields the exponential decay of the H 1 -norm of the state for initial states with sufficiently small H 1 -norm.
In Section 7, numerical experiments are presented that illustrate the theoretical findings. We also show simulations for the original quasilinear model with the suggested boundary feedback that indicate that also in this case the system decays exponentially.

The isothermal Euler equations
The isothermal Euler equations as a model for the flow through gas pipelines have already been stated in [2]. We use the model for real gas as described in [17]. Let a directed star-shaped graph G = (V, E) of a pipeline network be given. Here V denotes the set of vertices and E denotes the set of edges. Each edge e ∈ E corresponds to an interval [0, L e ] that represents a pipe of length L e > 0. Let D e > 0 denote the diameter of the pipe and λ e f ric > 0 the friction coefficient in pipe e. Define θ e = λ e f ric D e . Let ρ e denote the gas density, p e the pressure and q e the mass flow rate. Let α e ∈ (−0.25, 0] be given and define the compressibility factor as z e (p e ) = 1 + α e p e . (2.1) We assume that α e is independent of e. Equation (2.1) is also stated in [22] as the model of the American Gas Association (AGA). In [5] it is stated that it is sufficiently accurate within the network operating range. Let R e s denote the gas constant and T e the temperature that are independent of e. We assume that We study the isothermal Euler equations that govern the flow through a single pipe.

The Node conditions for the network flow
In this section we introduce the coupling conditions that model the flow through the nodes of the network. The node conditions that determine the flow dynamics are given in [2,10]. At the vertices v ∈ V , the flow is governed by the node conditions that require the conservation of mass and the continuity of the pressure. Let E 0 (v) denote the set of edges in the graph that are incident to v ∈ V and x e (v) ∈ {0, L e } denote the end of the interval [0, L e ] that corresponds to the edge e that is adjacent to v. We assume that at the central node v ∈ V of the star-shaped graph we have x e (v) = 0 for all e ∈ E 0 (v).
The continuity of the pressure at v means that for all e, f ∈ E 0 (v) we require the equation The conservation of mass is guaranteed by the Kirchhoff condition

The system in terms of Riemann invariants
As pointed out in [17], for e ∈ E the Riemann invariants R e − , R e + of the system are given by R e − (p e , q e ) = ln(p e ) − R e s T e q e p e (1 + α e p e ) , R e + (p e , q e ) = ln(p e ) + R e s T e q e p e (1 + α e p e ) .
For the central node v of the star-shaped graph (that corresponds to x = 0 for all e ∈ E) the node conditions (2.4), (2.5), can be written in the form of the linear equation This can be seen as follows. Equation (3.1) implies that for all e ∈ E, the value of R e + (t, 0) + R e − (t, 0) is the same, which implies that the value of ln(p e ) is independent of e, hence (2.4) holds. Moreover, (3.1) implies Due to (2.4) this implies that equation (2.5) holds. In this case also for the difference of the squares of the Riemann invariants we have For a boundary node v of our star-shaped graph and e ∈ E 0 (v) we have x e (v) = L e . We state the Dirichlet boundary conditions in terms of Riemann invariants in the form Define the number For e ∈ E, let ν e = 1 8 R e s T e θ e be given and define Define∆ e as the diagonal 2 × 2 matrix that contains the eigenvalues (3.7) In terms of the Riemann invariants, the quasilinear system (2.3) has the following diagonal form: In order to simplify the model, we replace the eigenvalues by This definition implies that λ e − = −λ e + . Moreover, for all e, f ∈ E we have λ e + = λ f + Define ∆ e as the diagonal 2 × 2 matrix that contains the eigenvalues λ e + and λ e − . The approximation ofλ e + by λ e + andλ e − by λ e − is justified by the fact that in the practical applications, the velocity of the gas flow is much smaller than the sound speed. Indeed, in typical applications, the fluid velocity is several meters per second while the speed of sound is several hundred meters per second. Whereas v can be neglected relative to c, the friction term cannot be neglected as this would cause a large relative error. In this way, we obtain a semilinear model. We do not claim that solutions to the isothermal Euler equations and the semilinear system are close to each other for all times, but we do expect that solutions to both systems share important qualitative features such as decay to steady states. This expectation is in agreement with several numerical simulations that we have carried out, some of which are presented in Section 7. Let us note that the difference between the models becomes smaller the closer solutions get to the equilibrium.
With the diagonal matrix ∆ e , the semilinear model that we will consider in the remainder of this work has the following form: Note that for given (R e + , R e − ) we have In particular this implies p e > 0. On account of the physical interpretation of the pressure it is very desirable that for the solutions we have p e > 0. This is an advantage of the model that is given by system (S). A similar semilinear model for gas transport has been studied in [18] in the context of identification problems. The model in [18] has the disadvantage that the matrix of the linearization of the source term is indefinite. However, the results from [18] can be adapted to the model that we consider in this paper.

A well-posedness result
In the semilinear model that we consider, the constant eigenvalues in the diagonal system matrix define two families of characteristics with constant slopes c and −c. For e ∈ E, define the sets For t ≥ 0, e ∈ E and the space variable x ∈ [0, L e ] we define the R 2 -valued function ξ e ± (s, x, t) as the solution of the initial value problem ξ e ± (t, x, t) = (t, x), ∂ s ξ e ± (t, x, t) = (1, ±c).

Define the points
For the t-component of P e± 0 (t, x) we use the notation t e± 0 (t, x) ≥ 0. The solution of (S) can be defined by rewriting the partial differential equation in the system in the form of integral equations along these characteristic curves, that is Note that almost everywhere the values of R e ± (P e± 0 (t, x)) are given on Γ e ± either by the initial data, that is is zero, this case only occurs for R e − ) or else by the node conditions (3.1) at x = 0. For a finite time interval [0, T ], the characteristic curves that start at t = 0 with the information from the initial data reach a point at the terminal time after a finite number of reflections at the boundaries x = L e (e ∈ E) or the central node x = 0.
The definition of the solutions of semilinear hyperbolic boundary value problems based upon (4.1) is described for example in [3]. For L ∞ -solutions, we have the following theorem.
there exists a unique solution of (S) that satisfies the integral equations (4.1) for all e ∈ E along the characteristic curves with R e + , R e − ∈ L ∞ ((0, L e ) × (0, T )) (e ∈ E) and the boundary conditions at x = L e and the node conditions at x = 0 almost everywhere in [0, T ] such that for all e ∈ E we have This solution depends in a stable way on the initial and boundary data in the sense that for initial data for the corresponding solution S e ± we have the inequality The proof is based upon Banach's fixed point theorem with the canonical fixed point iteration. It has to be shown that this map is a contraction in the Banach space 2 on a set of the form In order to show this, it suffices to derive an upper bound for the source term in (4.1) that is given by the continuously differentiable function σ e (R e + , R e − ). Moreover, it has to be shown that the iteration map maps from B(M ) into B(M ). This is true if M and ε are chosen sufficiently small.
In this analysis, it has to be taken into account that the characteristic curves can cross the central node at x = 0 a finite number of times. Due to the linear node condition (3.1) in each such crossing the absolute value of the outgoing Riemann invariants can be at most three times as large as the largest absolute values of the ingoing Riemann invariants. For t ∈ [0, T ] almost everywhere we have the inequality The last assertion follows from the integral equations satisfied by R e ± − S e ± and an application of Gronwall's Lemma.

Remark 4.2. An analogous existence result holds for solutions in
, (e ∈ E) for initial and boundary data that are compatible with each other and with the node conditions such that y e ± − J and u e − J have sufficiently small H 1 -norms.

An L 2 -observability inequality for a star-shaped network
In this section we derive an observability inequality for a star-shaped network. The aim is to get an upper bound for the L 2 -norm of the system state at the time t in terms of the norms of complete observations at the boundary nodes on the time interval [t − T, t + T ] with T > 0 sufficiently large. For such an inequality, the observation have to be integrated over a sufficiently long time-interval. In [6], an observability inequality for a star-shaped network of strings (without source term) is derived.
Such an observability is of interest since in the gas networks, usually very few sensors are available. The observability inequality implies that if the state is measured at the boundary nodes for a sufficiently long time, from these observations the system state can be determined completely.
and that there exists a constantM such that for all e ∈ E and x almost everywhere in [0, L e ] for the solution of (S) we have the inequality Then there exists a constant C 0 (M ) such that for all t > T , we have the inequality Remark 5.2. The proof of Theorem 5.1 is based upon the particular structure of the source term σ e . In fact, we use the properties that σ e is an increasing function of R e + − R e − with value zero if R e + − R e − = 0 that is Lipschitz continuous on bounded intervals.
Since (4.2) with M = 1 2M implies (5.1) , Theorem 4.1 yields sufficient conditions for (5.1) ifM is sufficiently small. An a priori upper bound (5.1) can also be obtained in the framework of semi-global classical solutions (even in the sense of a maximum), see [19]. Also if there is a solutions with C([0, T ], H 1 (0, L e )) regularity for all e ∈ E on [t − T, t + T ] × [0, L e ], then this solution automatically satisfies (5.1).
Then we have H e (0) = 0. For the derivative of H e with respect to x we have almost everywhere Due to the partial differential This implies the inequality for all x ∈ [0, L e ]. Since for real numbers r 1 , r 2 we have ( Thus we have Hence using the definition (5.3) of H e (x) and inequality (5.4) we obtain Now we prove the analogous inequality for R e − . We have By integration on [0, L e ] this yields the inequality Similarly as above by transformation of the 2-dimensional integral over a triangle this yields

An H 1 -semi-norm-observability inequality
Now we prove an observability inequality for the time derivative. Together with the L 2 -observability inequality this will finally yield an observability inequality for the H 1 -norm. Theorem 5.3. Let a real number J be given. Assume that T > max e∈E L e c and t > T . Assume that system (S) with u e (t) = J has a solution on [0, t + T ] such that for all e ∈ E and x almost everywhere in [0, L e ] we have R e + (·, x), R e − (·, x), ∂ t R e + (·, x), ∂ t R e − (·, x) ∈ L 2 (0, T + t) and that there exists a constantM such that for all e ∈ E and x almost everywhere in [0, L e ] for the solution of (S) we have the inequality (5.1) for s almost everywhere in [0, t + T ]. Then there exists a constant C 1 (M ) such that for all t > T , we have the inequality Then we have K e (0) = 0. For the derivative of K e with respect to x we have almost everywhere Due to the partial differential equation in system (S) this yields d dx K e (x) This implies the inequality for x ∈ [0, L e ] almost everywhere. With (5.1), due to Gronwall's Lemma this implies the inequality for all x ∈ [0, L e ]. Since for real numbers r 1 , r 2 we have (r 1 + r 2 ) 2 ≤ 2 r 2 1 + 2 r 2 2 , for x ∈ [0, L e ] almost everywhere we obtain Thus we have Hence using the definition (5.7) of K e (x) and inequality (5.8) we obtain Now we prove the analogous inequality for R e − . We have By integration on [0, L e ] this yields the inequality Similarly as above by transformation of the 2-dimensional integral over a triangle this yields

Stability of the state on the network
In this section we present a Dirichlet boundary feedback that can be used to stabilize a constant state of system (S) exponentially fast. Note that in physical terms the constant state of system (S) are the states where the gas is at rest. In terms of Riemann invariants, this means that both Riemann invariants are equal. Hence in terms of Riemann invariants, the constant states are given by real numbers J that are equal to the logarithm of the (constant) pressure. Theorem 6.1 has two parts. In the first part a sufficient condition for the exponential decay of the L 2norm (6.2) on a finite time interval [0, T ] is provided under the assumption (5.1). This first part of Theorem 6.1 can be applied to L ∞ -solutions as discussed in Theorem 4.1. For the proof of the first part, the observability inequality from Theorem 5.1 is used.
In the second part of Theorem 6.1 more regular H 1 solutions are considered. Together with (6.2), (6.3) implies the exponential decay of the H 1 -norm and allows a globalization in time. is exponentially stable in the sense that there exist constants C 1 > 0 and µ 0 > 0 such that for all t ∈ [0, T ] Hence the L 2 -norm of the difference between the constant state J and the actual state decays exponentially fast. Assume in addition that y e ± − J has a sufficiently small H 1 -norm and is compatible with the node condition and the boundary conditions such that a solution of system (S) exists on [0, T ] in × e∈E C([0, T ], H 1 (0, L e )) and satisfies (5.1). Then in addition to (6.2) also the L 2 -norm of the time-derivatives decay exponentially fast in the sense that there exist constants C 2 > 0 and µ 2 > 0 Remark 6.2. In terms of the physical variables the assertion in Theorem 6.1 implies that the state decays exponentially fast to a state with constant pressure p = exp(J). For all states with constant pressure we have the flow rates q e = 0. This means that the gas is at rest, that is the velocity is zero. It is important to note that for the exponential decay of the L 2 -norm we do not need any restrictions on the lengths of L e . In the frictionless case (that is for ν e = 0) System (S) can be seen as a system of coupled vibrating strings with absorbing boundary conditions. Similarly as in [13], it can be shown that in this case the system is steered to a position of rest in finite time.
If T is chosen sufficiently large andM is sufficiently small, the exponential decay of the L 2 -norm of the state and of the derivatives in (6.3) implies that the H 1 -norm of the terminal state at time T is smaller than the H 1 -norm of the initial state. This implies that the solution can be extended to the time interval [0, 2T ] and so forth. In this way, for a sufficiently small norm of the initial state, we obtain a solution of the closed loop system that is global in time.
First we state a Lemma in order to clarify the relation between the regularity and the exponential decay of ∂ t R e and ∂ x R e : Lemma 6.3. Let (R e + , R e − ) denote a solution of (S) such that for all e ∈ E we have R e ± ∈ C([0, T ], H 1 (0, L e )) and (5.1) holds. Then for all e ∈ E and t ∈ [0, T ] we have ∂ t R e ± ∈ L 2 (0, L e ). Moreover, we have the inequality Hence if ∂ t R e ± (t, ·) L 2 (0, L e ) and R e ± (t, ·) − J L 2 (0, L e ) decay exponentially fast, also ∂ x R e ± decays exponentially fast.
Proof. We have ∂ t R e ± = −(±c)∂ x R e ± − (±σ e (R e + , R e − )). Since σ e (R e + , R e − )(t, ·) ∈ L ∞ (0, L e ), the first assertion follows. Moreover, we have Hence since R e ± ∈ L 2 (0, L e ), the triangle inequality, (5.1) and the definition of σ e imply inequality (6.4). Proof of Theorem 6.1. Lett ∈ (0, t) be given. For e ∈ E the partial differential equation implies that We multiply this equation by (R e + − J) and integrate over the interval (t −t, t +t) to obtain

This yields
Similarly, we obtain For e ∈ E and t ∈ [0, T ], define the function Then we have Hence due to the definition of σ e we obtain the inequality Due to the node conditions at x = 0 and (3.2) we have Hence due to the boundary condition at x = L e we have In particular, we have L 0 (t +t) ≤ L 0 (t −t). Since the above inequality can be derived for allt ∈ (0, t), this implies in particular that L 0 is decreasing. Now we chooset = T 0 > 0. For all e ∈ E the observability inequality (5.2) implies with C 0 (M ) as defined in (5.5). This yields the inequality Since L 0 is decreasing this yields Hence we have Similar as in Lemma 2 from [16], this implies that L 0 decays exponentially fast. Now we consider the evolution of the time-derivatives. For e ∈ E and t ∈ [0, T ], we consider Due to the partial differential equation in system (S), for solutions with H 2 -regularity we have Lett ∈ (0, t) be given. For e ∈ E the partial differential equation implies that We multiply this equation by ∂ t R e + and integrate over the interval (t −t, t +t) to obtain This yields Similarly, we obtain Hence we have Thus we obtain the inequality Due to the node conditions at x = 0 we have e∈E (D e ) 2 ∂ t R e + (τ, 0)) 2 − (∂ t R e − (τ, 0)) 2 = 0.
Hence due to the boundary condition at x = L e we have In particular, we have L 1 (t +t) ≤ L 1 (t −t). Since the above inequality can be derived for allt ∈ (0, t), this implies in particular that L 1 is decreasing. Now we chooset = T 0 > 0. For all e ∈ E the observability inequality (5.6) implies with C 1 (M ) as defined in (5.9). This yields the inequality Since L 1 is decreasing this yields Hence we have Similar as in Lemma 2 from [16], this implies that L 1 decays exponentially fast.
In this inequality, the second derivatives do not appear and the constants are also independent of the second derivatives. Therefore, by a density argument we conclude that it also holds for solutions with H 1 -regularity. Hence (6.3) follows.

Numerical experiments
In this section we will conduct numerical experiments in order to illustrate the theoretical findings. For the semi-linear model, we have used a finite difference upwind discretisation for the convective terms. The temporal discretisation of the convection is explicit Euler while the friction terms were discretized using implicit Euler. This allows us to use the maximal time step allowed for by the CFL condition so that discontinuities in the solution are not smoothed out by artificial diffusion. For the quasi-linear model (2.3), we have used a Lax-Friedrichs flux in space and the same discretisation in time as for the semi-linear model. In each of the experiments we have plotted L 0 (see Eq. (6.6)) over time. In order to test how sharp the decay estimates (6.2) are, we have also plotted In all computations the initial velocity field is constant zero. Our results confirm the solutions decay to equilibrium at least exponentially. Indeed, the experimentally observed decay rates are (as expected) larger than the lower bound from our analysis.

Discontinuous initial data with friction
In this experiment the initial pressure on pipes 1 and 3 is constant 60bar whereas in pipe 2 we have 60bar on the half of the pipe adjacent to the node and 80bar on the half of the pipe adjacent to the boundary. The results are displayed in Figure 1. The picture shows that, as predicted, L 0 decays exponentially. However, in our analysis we only obtain a lower bound for the decay rate. The picture shows that the observed rate is larger than our lower bound that holds for all initial data.
We also plot snapshots of the numerical solutions at t = 0s, t = 28s and t = 56s in Figure 2. Note that we have R e + (0) = R e − (0) on all pipes. It can be seen that the discontinuities in the solution remain, since there is no diffusion, but the plateau values change due to the friction terms.

Discontinuous initial data without friction
We use the same parameters as in the previous section but set friction to zero. The results are displayed in Figure 3. We can see that, as expected, L 0 goes to zero in finite time, see Remark 6.2.
We also plot snapshots of the numerical solutions at t = 0s, t = 28s and t = 56s in Figure 4. It can be seen that the discontinuities in the solution are not smoothed out and plateau values only change via interaction with nodes.

Continuous initial data with friction
In this experiment, the initial pressure on pipes 1 and 3 is p(x) = 60bar + 20bar sin(πx/km) while the initial pressure on pipe 2 is p(x) = 60bar + 20bar sin(πx/(2km)). All other parameters are as above. The results concerning decay of L 0 and L 1 are displayed in Figure 5. Note that in this case we have computed L 0 and L 1 for both the semilinear model (S) and the quasilinear isothermal Euler equations (2.3). Snapshots of solutions to both systems are shown in Figures 6 and 7, respectively. Figure 5 shows that L 0 and L 1 decay to zero exponentially. In both cases, the actually observed decay rates are larger than the lower bounds that we have proven. Our experiments, in particular, the snapshots at t = 28, show the difference between the solutions to the semi-and quasilinear systems due to the different wave speeds.   . Continuous initial data with friction: Snapshots of the semilinear solution at times t = 0s, t = 28s and t = 56s. Pipes 1 and 3 go along the positive and negative x-axis, pipe 2 goes along the positive y-axis. On all pipes R + is a yellow continuous line and R − is a blue dashed line.
Despite the quantitative differences, the solutions share important qualitative features such as exponential decay to the steady state. Indeed, Figure 5 shows that the isothermal Euler solution converges to the steady state nearly as quickly as the semilinear solution.

Continuous initial data without friction
We study the case without friction. In this case, it can be shown that both L 0 (t) and L 1 (t) vanish (for the semilinear solution) for t > 2 · T 0 ≈ 59, see Remark 6.2. Numerical results showing the decay of L 0 and L 1 for both the semilinear and the quasilinear (isothermal Euler) system are displayed in Figure 8. Snapshots of the solutions to both systems are shown in Figures 9 and 10, respectively.
We observe, as expected, for the semilinear solution, both L 0 and L 1 have decayed to zero at t = 2T 0 ≈ 59. Our experiments, in particular, the snapshots at t = 28s and t = 56s, show the difference between both solutions due to the different wave speeds. It should be noted that, without friction, the semilinear model (S) is, in fact, linear and the evolution equations for both Riemann invariants are only coupled via the boundary conditions while in the isothermal Euler system there is coupling via the flux term. However, although there are quantitative differences, the solutions share important qualitative features such as exponential decay to the steady state. Indeed, Figure 5 shows that the solution to the isothermal Euler equations (2.3) converges to the steady state even quicker than the solution (S). It should be noted, however, that for the isothermal Euler solution L 1 is not monotone which makes sense since for quasilinear models (such as isothermal Euler) certain slopes grow steeper which, given sufficient time, might lead to the appearance of shocks. Figure 9. Continuous initial data without friction: Snapshots of the semilinear solution at times t = 0s, t = 28s and t = 56s. Pipes 1 and 3 go along the positive and negative x-axis, pipe 2 goes along the positive y-axis. On all pipes R + is a yellow continuous line and R − is a blue dashed line. Figure 10. Continuous initial data without friction: Snapshots of the quasilinear solution at times t = 0s, t = 28s and t = 56s. Pipes 1 and 3 go along the positive and negative x-axis, pipe 2 goes along the positive y-axis. On all pipes R + is a yellow continuous line and R − is a blue dashed line.

Conclusions
We have considered the flow of gas in a star-shaped network of pipelines that is governed by a hyperbolic semilinear model of partial differential equations that can be understood as a simplification of the isothermal Euler equations. The model under consideration is useful for simulations concerning the operation of gas networks. We have shown that locally around a state where the gas is at rest the system state can be steered towards this position of rest exponentially fast with a suitably chosen boundary feedback control. In fact the boundary control consists of absorbing boundary conditions. This result is remarkable, since we show the exponential decay of the H 1 -norm of the state for a system with a source term on a network. In the earlier work [12] only the decay of the L 2 -norm on a star-shaped network has been considered for the quasilinear model. In the present paper, we have shown the exponential decay of the L 2 -norm on the network without additional assumptions on the lengths of the pipes. Also the H 1 -norm decays exponentially on the network for arbitrarily long pipes if the H 1 -norm of the initial state is sufficiently small. Our theoretical findings are in agreement with numerical experiments that also indicate that the solutions to the semilinear model share certain qualitative features with solutions of the isothermal Euler equations. This analysis is motivated by the problem of flow reversal that occurs in certain scenarios in gas pipeline network operations. In order to avoid possible turbulence, in this situation it is a good strategy to control the flow to rest first and then start it again in the opposite direction.
We expect that the result can be extended to more general tree-shaped graphs. This will be the subject of future research. We will also investigate to which extent similar estimates can be proven for (entropy) solutions of the quasilinear isothermal Euler equations. The treatment of discontinuous weak solutions will require new tools that go beyond what we have developed here.