1 Introduction

Gas transport through long pipes is usually dominated by friction at the pipe walls. For large friction \(O(1/\varepsilon ^2)\) and on the practically relevant time scales \(t=O(1/\varepsilon )\), the governing balance equations for mass and momentum can be phrased in as

$$\begin{aligned} a \partial _{\tau }\rho + \partial _xm&= 0 \end{aligned}$$
(1)
$$\begin{aligned} \varepsilon ^2 \partial _{\tau }w + \partial _xh + \gamma |w| w&= 0. \end{aligned}$$
(2)

Here a is the cross section of the pipe, \(\rho \) is the density of the gas, and w, \(\tau \), \(\gamma >0\) denote the rescaled velocity, time, and friction coefficient, respectively. Furthermore

$$\begin{aligned} m = a \rho w \qquad \text {and} \qquad h = \frac{\varepsilon ^2 w^2}{2} + P'(\rho ) + g z \end{aligned}$$
(3)

are the rescaled mass flow rate and the total specific enthalpy with \(P(\rho )\) denoting the pressure potential, g the gravity constant, and z the elevation of the pipe. The system (1)–(3) is equivalent to the barotropic Euler equations with friction and hence of quasi-linear hyperbolic type. For convenience of the reader, a detailed derivation of the above equations, which are assumed to hold for \(0< x < \ell \) and for time \(\tau \ge 0\), is given in the appendix.

In rescaled variables, the physical energy of the system is given by

$$\begin{aligned} \mathcal {H}(\rho ,w)&= \int _0^\ell a \left( \varepsilon ^2 \frac{\rho w^2}{2} + P(\rho ) + g z \rho \right) \, dx. \end{aligned}$$
(4)

One can observe that the state variables \((\rho ,w)\) and the co-state variables (h, m), defined in Eq. (3), are directly linked via the derivative of this energy functional; details will be given below. As a direct consequence of the particular problem structure, regular solutions of (1)–(3) can be shown to satisfy the following energy-dissipation identity

$$\begin{aligned} \frac{d}{d\tau }\mathcal {H}(\rho ,w) + \mathcal {D}(\rho ,w)&= -m h \big |_0^\ell \end{aligned}$$
(5)

with dissipation functional \(\mathcal {D}(\rho ,w) = \int _0^\ell \gamma \rho |w|^3 \, dx \ge 0\). The free energy of the system thus changes only due to friction at the pipe walls and energy transfer across the boundary.

To fully describe the evolution, the problem has to be complemented by appropriate initial and boundary conditions. For instance, one may prescribe

$$\begin{aligned} h = h_\partial \qquad \text {on } \{0,\ell \}; \end{aligned}$$
(6)

alternative conditions will be briefly discussed later.

1.1 Scope and summary of main results

In this paper, we investigate the stability of solutions to the nonlinear system (1)–(3) with respect to perturbations of the initial and boundary values \(h_\partial \) in (6) as well as on the model parameters \(\varepsilon \) and \(\gamma \). Our analysis is done under the following structural assumptions, which characterize standard operating conditions for gas transport in pipeline networks.

  1. (A1)

    The pressure potential \(P: \mathbb {R}_+ \rightarrow \mathbb {R}\) is smooth and strictly convex and there exist positive constants \(\underline{\rho },\bar{\rho }, {\bar{w}}\) and \({\bar{\varepsilon }}\) such that

    $$\begin{aligned} \rho P''(\rho ) \ge 4 {\bar{\varepsilon }}^2 |{\bar{w}}|^2 \qquad \forall \underline{\rho }\le \rho \le {\bar{\rho }}. \end{aligned}$$
    (7)

    Moreover, \(0<\underline{a} \le a \le {\bar{a}}\) and \(|gz| \le {\bar{g}}{\bar{z}}\) for some constants \(\underline{a},{\bar{a}}\), and \({\bar{g}} {\bar{z}}\).

  2. (A2)

    Sufficiently smooth solutions \((\rho ,w)\) and \((\hat{\rho },{\hat{w}})\) to the system (1)–(3) exist for parameters \(0 \le \varepsilon , {\hat{\varepsilon }} \le {\bar{\varepsilon }}\) and \(0 < \underline{\gamma }\le \gamma ,{\hat{\gamma }} \le {\bar{\gamma }}\) with \(\underline{\gamma },{\bar{\gamma }}\) constant, which additionally satisfy

    $$\begin{aligned} \underline{\rho } \le \rho ,{\hat{\rho }} \le {\bar{\rho }} \qquad \text {and} \qquad - {\bar{w}} \le w,{\hat{w}} \le {\bar{w}}. \end{aligned}$$
    (8)

Note that the pressure potential \(P=P(\rho )\) plays the role of an abstract state equation with \(c=\sqrt{p'(\rho )}=\sqrt{\rho P''(\rho )}\) denoting the density dependent speed of sound and \(p=p(\rho )\) the pressure law. Conditions (A1) and (A2) imply subsonic flow and uniform bounds for the states that are bounded away from vacuum; see [6] for further details.

Under the conditions stated above, we will prove the following general stability result: Let \((\rho ,w)\) and \(({\hat{\rho }},{\hat{w}})\) be sufficiently regular solutions of (1)–(3) with parameters \(\varepsilon ,\gamma \) and \({\hat{\varepsilon }},{\hat{\gamma }}\), and boundary values \({\hat{h}}_\partial \) and \(h_\partial \) as described in (6). Then

$$\begin{aligned}&\Vert \rho (\tau ) - {\hat{\rho }}(\tau )\Vert _{L^2}^2 + \varepsilon ^2 \Vert w(\tau ) - {\hat{w}}(\tau )\Vert ^2_{L^2} + \int _0^\tau \Vert w(s)-{\hat{w}}(s)\Vert ^3_{L^3} ds \\&\quad \le {\hat{C}} e^{{\hat{c}} \tau } \Big (\Vert \rho (0) - {\hat{\rho }}(0)\Vert _{L^2}^2 + \Vert \varepsilon w(0) - \varepsilon {\hat{w}}(0)\Vert _{L^2}^2 \\&\qquad + |\gamma - {\hat{\gamma }}|^{3/2} + |\varepsilon ^2 -{\hat{\varepsilon }}^2| + \int _0^\tau |h_\partial (s) - {\hat{h}}_\partial (s)| \, ds \Big ) \end{aligned}$$

with constants \({\hat{c}},{\hat{C}}\) depending only on the bounds in assumptions (A1)–(A2) and on bounds for time derivatives of \(\hat{\rho }\) and \({\hat{w}}\). A precise statement of our results and of the regularity assumption on the solutions is given in Sect. 4.

Let us emphasize that the energy functional (4), the dissipation functional (5), and also the co-state variables (3) depend explicitly on the parameters \(\varepsilon \) and \(\gamma \), and hence their definition changes for the perturbed problem. As two immediate consequences of the above stability estimate, we obtain

  • uniqueness of regular subsonic bounded state solutions for specified initial and boundary values, and their stability with respect to perturbations in these problem data as well as the model parameters;

  • convergence of solutions to those of the parabolic limit problem which results by formally setting \(\varepsilon =0\) in equations (1)–(6).

In contrast to the hyperbolic problem under consideration, the existence and uniqueness of solutions to the parabolic limit problem for \(\varepsilon =0\) can be proven rigorously by variational arguments; see [1, 25] or [27]. This parabolic problem also serves as the basis for simulation codes utilized in the gas network community [2, 24]; our stability and asymptotic analysis yields a quantitative justification for the use of this model. Finally, due to the energy-based modelling, our results can be generalized almost verbatim to gas networks by utilizing appropriate coupling conditions; details will be given in Sect. 5.

1.2 Main tools

The proof of our main results is based on the observation that (1)–(3) can be written as an abstract dissipative Hamiltonian system

$$\begin{aligned} \mathcal {C}\partial _{\tau }{{\varvec{u}}}+ (\mathcal {J}+\mathcal {R}({{\varvec{u}}})) {{\varvec{z}}}({{\varvec{u}}})&= \mathcal {B}_\partial {{\varvec{z}}}({{\varvec{u}}}), \end{aligned}$$
(9)

with \({{\varvec{u}}}=(\rho ,w)\) and \({{\varvec{z}}}({{\varvec{u}}})=(h,m)\) denoting the state and co-state variables, which are linked via the energy functional \(\mathcal {H}({{\varvec{u}}})=\mathcal {H}(\rho ,w)\); moreover, \(\mathcal {C}\), \(\mathcal {J}\), \(\mathcal {R}({{\varvec{u}}})\), and \(\mathcal {B}_\partial \), are appropriate operators, the last one incorporating the boundary conditions; see Sect. 2. We refer to [23, 28, 29] for similar structures arising in general port-Hamiltonian systems.

Any smooth function \({\widehat{{{\varvec{u}}}}}=({\hat{\rho }},{\hat{w}})\), e.g., a solution of (1)–(3) with perturbed model parameters, may be considered to be a solution of the perturbed system

$$\begin{aligned} \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}+ (\mathcal {J}+ \mathcal {R}({\widehat{{{\varvec{u}}}}})) {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})&= \mathcal {B}_\partial {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) + {\widehat{{{\varvec{r}}}}}, \end{aligned}$$
(10)

with the same operators \(\mathcal {C}\), \(\mathcal {J}\), \(\mathcal {R}(\cdot )\), \(\mathcal {B}_\partial \), and the same energy functional \(\mathcal {H}(\cdot )\) and state to co-state mapping \({{\varvec{z}}}(\cdot )\), up to some perturbation \({\widehat{{{\varvec{r}}}}}\). The coincidence of the underlying Hamiltonian structure will allow us to estimate the difference between \({{\varvec{u}}}\) and \({\widehat{{{\varvec{u}}}}}\) in terms of the perturbations \({\widehat{{{\varvec{r}}}}}\) by means of relative energy estimates. For the stability analysis on the abstract level, we require some general conditions that will be verified in detail for the gas transport problem under investigation using the assumptions (A1)–(A2).

The use of an abstract problem structure greatly simplifies the analysis and allows to generalize our results in various directions. We briefly discuss the incorporation of perturbations in other model parameters and the extension to gas networks.

1.3 Review of related work

Relative entropy or energy techniques have been used intensively for the existence, stability, and discretization error analysis of time dependent partial differential equations; we refer to [17] for a recent summary of corresponding results for parabolic evolution problems. In the present paper, we are interested in hyperbolic problems, where the use of relative entropy arguments goes back to the seminal works of DiPerna [7] and Dafermos [5]; also see [6] for an introduction to the field. Typical aspects addressed are: convergence to steady states, stable dependence of solutions on initial data and parameters, and asymptotic limits. Examples for the latter include the low Mach limit of the Euler and Navier–Stokes equations, which are investigated in [10] for instance. Long time convergence of solutions to damped Euler equations to Barenblatt solutions were investigated by Huang and coworkers in a series of papers [11].

One particular aspect that we want to address in the present study are parabolic limits of quasilinear hyperbolic equations; see [16, 18, 19, 22] for some exemplary results in this direction. The latter reference as well as [4, 12] strongly rely on the underlying dissipative Hamiltonian structure of many equations in fluid mechanics, which will also play a major role in our analysis below. The previous papers, however, use formulations in conservative variables, which allows to deal with one solution being a weak solution in the classical sense of hyperbolic conservation laws [6]. Parabolic limits of hyperbolic \(2\times 2\) systems have also been studied in [21, 31] using compensated compactness arguments [14, 15]. This has the advantage that no additional regularity of solutions to the parabolic limit problem is required, but no quantitative information about the speed of convergence is obtained. Spectral estimates for the linear part of the problem are used in [8] to derive convergence in the parabolic limit.

In contrast to the work mentioned above, our study is based on a formulation in primitive variables, which requires both solutions to have some minimal smoothness. On the other hand, this formulation allows us to incorporate boundary conditions more naturally and to extend our results to networks in a straight-forward manner using appropriate coupling conditions that guarantee conservation or dissipation of energy at network junctions [9, 26]. The proper handling of the network junctions plays an important role and therefore many results that are available for problems in \(\mathbb {R}^n\) or on periodic domains cannot be applied to the network setting without significant modifications.

Similar formulations for compressible flow were also considered in the context of port-Hamiltonian systems; see [13, 28, 29] for modeling and [3, 20, 23] for corresponding discretization strategies. In contrast to these works, which emphasize the geometric modelling of the underlying problems, we here focus on the derivation of analytical results taking advantage of the specific problem structure. Other systems that fit into the general framework developed in this paper include the Euler-Korteweg system, the system of quantum hydrodynamics and the Euler-Poisson equations; see [12] for an overview.

1.4 Outline of the manuscript

The remainder of the manuscript is organized as follows: Sect. 2 is concerned with the perturbation analysis for the abstract system (9). In Sect. 3, we then briefly review the derivation of the model equations (1)–(3) starting from the usual balance equations of gas dynamics, and we show that they fit into the abstract form (9). In Sect. 4, we verify the assumptions required for our abstract analysis for the gas transport problem under investigation, which allows us to state and prove our main results. Their generalization to gas networks is discussed in Sect. 5.

2 An abstract stability estimate

In this section, we consider abstract evolution problems of the form

$$\begin{aligned}&\mathcal {C}\partial _{\tau }{{\varvec{u}}}+ (\mathcal {J}+ \mathcal {R}({{\varvec{u}}})) \, {{\varvec{z}}}({{\varvec{u}}}) = \mathcal {B}_\partial {{\varvec{z}}}({{\varvec{u}}}), \end{aligned}$$
(11)
$$\begin{aligned}&{{\varvec{z}}}({{\varvec{u}}}) = \mathcal {C}^{-1} \mathcal {H}'({{\varvec{u}}}), \end{aligned}$$
(12)

with state and co-state variables \({{\varvec{u}}}\) and \({{\varvec{z}}}({{\varvec{u}}})\) that are directly connected via the derivative \(\mathcal {H}'({{\varvec{u}}})\) of an associated energy functional \(\mathcal {H}({{\varvec{u}}})\). Let us note that non-linearity enters in the friction operator \(\mathcal {R}({{\varvec{u}}})\) through the non-linear state to co-state map \({{\varvec{z}}}={{\varvec{z}}}({{\varvec{u}}}).\) Abstract problems of similar structure are considered in [23, 29] under the name port-Hamiltonian systems. After briefly introducing a reasonable functional analytic setting, we derive an abstract stability estimate which will be the basis for our further considerations.

2.1 Notation and basic assumptions

Let \(\mathbb {W}\subset \mathbb {V}\) be real Hilbert spaces, with \(\mathbb {V}' \subset \mathbb {W}'\) denoting the corresponding dual spaces. We use \(\langle \cdot ,\cdot \rangle \) to denote the duality product on \(\mathbb {V}' \times \mathbb {V}\) and \(\mathbb {W}' \times \mathbb {W}\). Furthermore, let \(\mathcal {H}: \mathbb {D}\subset \mathbb {V}\rightarrow \mathbb {R}\) be a smooth and strictly convex energy functional, with \(\mathbb {D}\subset \mathbb {V}\) denoting some appropriate, e.g., convex and closed, subset. We consider abstract evolution problems of the form (11)–(12), with operators satisfying the following conditions.

Assumption 1

\(\mathcal {C}:\mathbb {V}\rightarrow \mathbb {V}'\) is linear, bounded, self-adjoint, and elliptic on \(\mathbb {V}\), i.e.,

$$\begin{aligned} \langle \mathcal {C}{{\varvec{u}}}, {{\varvec{v}}}\rangle = \langle \mathcal {C}{{\varvec{v}}}, {{\varvec{u}}}\rangle \qquad&\forall {{\varvec{u}}},{{\varvec{v}}}\in \mathbb {V}, \end{aligned}$$
(13)
$$\begin{aligned} c \Vert {{\varvec{v}}}\Vert _\mathbb {V}^2 \le \langle \mathcal {C}{{\varvec{v}}}, {{\varvec{v}}}\rangle \le C \Vert {{\varvec{v}}}\Vert _{\mathbb {V}}^2 \qquad&\forall {{\varvec{v}}}\in \mathbb {V}, \end{aligned}$$
(14)

with positive constants c, C independent of \({{\varvec{v}}}\). For any \({{\varvec{u}}}\in \mathbb {D}\) the operator \(\mathcal {R}({{\varvec{u}}}): \mathbb {W}\rightarrow \mathbb {W}'\) is linear, bounded, self-adjoint, and non-negative, in particular

$$\begin{aligned} \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{w}}}, {{\varvec{z}}}\rangle&= \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}, {{\varvec{w}}}\rangle \qquad{} & {} \forall {{\varvec{w}}},{{\varvec{z}}}\in \mathbb {W}, \end{aligned}$$
(15)
$$\begin{aligned} \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{w}}}, {{\varvec{w}}}\rangle&\ge 0 \qquad{} & {} \forall {{\varvec{w}}}\in \mathbb {W}. \end{aligned}$$
(16)

Furthermore, \(\mathcal {J}: \mathbb {W}\rightarrow \mathbb {W}'\) is linear and anti-symmetric, i.e.,

$$\begin{aligned} \langle \mathcal {J}{{\varvec{w}}}, {{\varvec{z}}}\rangle = - \langle \mathcal {J}{{\varvec{z}}}, {{\varvec{w}}}\rangle \qquad \forall {{\varvec{w}}},{{\varvec{z}}}\in \mathbb {W}, \end{aligned}$$
(17)

and finally, the operator \(\mathcal {B}_\partial : \mathbb {W}\rightarrow \mathbb {W}'\) is linear and bounded.

From conditions (13)–(14), one can immediately see that

$$\begin{aligned} \langle {{\varvec{u}}}, {{\varvec{v}}}\rangle _\mathcal {C}:=\langle \mathcal {C}{{\varvec{v}}}, {{\varvec{u}}}\rangle \end{aligned}$$
(18)

defines a scalar product on \(\mathbb {V}\) and the associated norm \(\Vert {{\varvec{v}}}\Vert _\mathcal {C}^2 = \langle {{\varvec{v}}},{{\varvec{v}}}\rangle _\mathcal {C}= \langle \mathcal {C}{{\varvec{v}}}, {{\varvec{v}}}\rangle \) is equivalent to the standard norm \(\Vert \cdot \Vert _\mathbb {V}\). The expression \({{\varvec{z}}}({{\varvec{u}}}) = \mathcal {C}^{-1} \mathcal {H}'({{\varvec{u}}}) = {\text {grad}}_\mathcal {C}\mathcal {H}({{\varvec{u}}})\) then denotes the gradient of the functional \(\mathcal {H}\) at \({{\varvec{u}}}\) with respect to this scalar product. We further introduce the symbol

$$\begin{aligned} \mathcal {G}({{\varvec{u}}}) = \mathcal {C}^{-1} \mathcal {H}''({{\varvec{u}}}) \end{aligned}$$
(19)

for the Hessian operator \(\mathcal {G}({{\varvec{u}}}): \mathbb {V}\rightarrow \mathbb {V}'\), and note that

$$\begin{aligned} \langle \mathcal {G}({{\varvec{u}}}) {{\varvec{v}}}, {{\varvec{w}}}\rangle _\mathcal {C}= \langle \mathcal {H}''({{\varvec{u}}}) {{\varvec{v}}}, {{\varvec{w}}}\rangle = \langle \mathcal {H}''({{\varvec{u}}}) {{\varvec{w}}}, {{\varvec{v}}}\rangle = \langle \mathcal {G}({{\varvec{u}}}) {{\varvec{w}}}, {{\varvec{v}}}\rangle _\mathcal {C}, \end{aligned}$$
(20)

i.e., the Hessian is symmetric with respect to the scalar product induced by \(\mathcal {C}\).

Notation 2

By a classical solution of (11)–(12) on [0, T], we mean a function

$$\begin{aligned} {{\varvec{u}}}\in C^1([0,T];\mathbb {V}) \cap C^0([0,T];\mathbb {D}) \quad \text {with} \quad {{\varvec{z}}}( {{\varvec{u}}}) \in C^0([0,T];\mathbb {W}) \end{aligned}$$

that satisfies (11)–(12) for all \(0 \le \tau \le T\) in the sense of \(\mathbb {W}'\).

2.2 Power balance

As a direct consequence of the underlying port-Hamiltonian structure and our assumptions on the operators, we obtain the following power balance relation.

Lemma 3

Let \({{\varvec{u}}}\) be a classical solution of (11)–(12). Then

$$\begin{aligned} \frac{d}{d\tau }\mathcal {H}({{\varvec{u}}}) = -\mathcal {D}({{\varvec{u}}}) + \langle \mathcal {B}_\partial {{\varvec{z}}}({{\varvec{u}}}), {{\varvec{z}}}({{\varvec{u}}})\rangle . \end{aligned}$$
(21)

with \(\mathcal {D}({{\varvec{u}}}):= \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}}), {{\varvec{z}}}({{\varvec{u}}})\rangle \) denoting the dissipation functional, i.e., the total energy of the system can only change via dissipation or power flowing over the ports.

Proof

By formal computation and using equations (11)–(12), we obtain

$$\begin{aligned} \frac{d}{d\tau }\mathcal {H}({{\varvec{u}}})&= \langle \mathcal {H}'({{\varvec{u}}}), \partial _{\tau }{{\varvec{u}}}\rangle = \langle \mathcal {C}\partial _{\tau }{{\varvec{u}}}, {{\varvec{z}}}({{\varvec{u}}})\rangle \\&= -\langle \mathcal {J}{{\varvec{z}}}({{\varvec{u}}}), {{\varvec{z}}}({{\varvec{u}}})\rangle - \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}}), {{\varvec{z}}}({{\varvec{u}}})\rangle + \langle \mathcal {B}_\partial {{\varvec{z}}}({{\varvec{u}}}), {{\varvec{z}}}({{\varvec{u}}})\rangle . \end{aligned}$$

The result then follows immediately from the assumptions on the operators. \(\square \)

2.3 Evolution of the relative energy

We now study the stability of solutions to (11)–(12) with respect to perturbations. To this end, let \({\widehat{{{\varvec{u}}}}}\) denote a classical solution of

$$\begin{aligned} \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}+ [\mathcal {J}+\mathcal {R}({\widehat{{{\varvec{u}}}}})] {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})&= \mathcal {B}_\partial {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) + {\widehat{{{\varvec{r}}}}}, \end{aligned}$$
(22)
$$\begin{aligned} z({\widehat{{{\varvec{u}}}}})&= \mathcal {C}^{-1} \mathcal {H}'({\widehat{{{\varvec{u}}}}}), \end{aligned}$$
(23)

with appropriate perturbation described by the residual functional \({\widehat{{{\varvec{r}}}}}\in C^0([0,T];\mathbb {W}')\). As a measure for the difference of \({{\varvec{u}}}\) and \({\widehat{{{\varvec{u}}}}}\), we use the relative energy [6], defined by

$$\begin{aligned} \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) := \mathcal {H}({{\varvec{u}}}) - \mathcal {H}({\widehat{{{\varvec{u}}}}}) - \langle \mathcal {H}'({\widehat{{{\varvec{u}}}}}), {{\varvec{u}}}-{\widehat{{{\varvec{u}}}}}\rangle . \end{aligned}$$
(24)

Using the particular problem structure and some elementary computations, we can prove the following basic identity for the temporal change of the relative energy.

Lemma 4

Let \({{\varvec{u}}}\), \({\widehat{{{\varvec{u}}}}}\) be classical solutions of (11)–(12) and (22)–(23). Then

$$\begin{aligned} \frac{d}{d\tau }\mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}})&= -\langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}}) - \mathcal {R}({\widehat{{{\varvec{u}}}}}) {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle + \langle \mathcal {B}_\partial ({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle \\&\quad +\langle \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) - \mathcal {G}({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}})\rangle + \langle {\widehat{{{\varvec{r}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle . \end{aligned}$$

Proof

By formal differentiation of the relative energy, we obtain

$$\begin{aligned}&\frac{d}{d\tau }\mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}})\\&\quad = \langle \mathcal {H}'({{\varvec{u}}}), \partial _{\tau }{{\varvec{u}}}\rangle - \langle \mathcal {H}'({\widehat{{{\varvec{u}}}}}), \partial _{\tau }{\widehat{{{\varvec{u}}}}}\rangle - \langle \mathcal {H}'({\widehat{{{\varvec{u}}}}}),\partial _{\tau }{{\varvec{u}}}- \partial _{\tau }{\widehat{{{\varvec{u}}}}}\rangle - \langle \mathcal {H}''({\widehat{{{\varvec{u}}}}}) \partial _{\tau }\widehat{{\varvec{u}}}, {{\varvec{u}}}- {\widehat{{{\varvec{u}}}}}\rangle \\&\quad = \langle \mathcal {H}'({{\varvec{u}}}) - \mathcal {H}'({\widehat{{{\varvec{u}}}}}), \partial _{\tau }{{\varvec{u}}}- \partial _{\tau }{\widehat{{{\varvec{u}}}}}\rangle + \langle \mathcal {H}'({{\varvec{u}}}) - \mathcal {H}'({\widehat{{{\varvec{u}}}}}) - \mathcal {H}''({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}}), \partial _{\tau }{\widehat{{{\varvec{u}}}}}\rangle \\&\quad = \langle \mathcal {C}\partial _{\tau }{{\varvec{u}}}- \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle + \langle \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) - \mathcal {G}({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}})\rangle . \end{aligned}$$

Here we used the symmetry of \(\mathcal {H}''({\widehat{{{\varvec{u}}}}})\) in the second, and the definitions of the gradient and Hessian operators in the last step. We now use (11) and (22) to replace the time derivatives in the first term, and arrive at

$$\begin{aligned} \frac{d}{d\tau }\mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}})&= -\langle \mathcal {J}({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle -\langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}}) - \mathcal {R}({\widehat{{{\varvec{u}}}}}) {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle \\&\quad + \langle \mathcal {B}_\partial ({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle + \langle {\widehat{{{\varvec{r}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle \\&\quad +\langle \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) - \mathcal {G}({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}})\rangle . \end{aligned}$$

Due to the anti-symmetry property (17) of the operator \(\mathcal {J}\), the first term vanishes, and we already obtain the assertion of the lemma. \(\square \)

2.4 An abstract stability result

We now derive quantitative estimates for the difference of the solutions \({{\varvec{u}}}\) and \({\widehat{{{\varvec{u}}}}}\) to (11)–(12) and (22)–(23) with respect to perturbations in the right hand side and the initial and boundary values. To do so, we make some abstract assumptions that will later be verified for the gas transport problem under consideration.

Assumption 5

There exist constants \({\hat{c}}_0={\hat{c}}_0(\mathbb {D})>0\), \( {\hat{C}}_0={\hat{C}}_0(\mathbb {D})>0\) such that

$$\begin{aligned} {\hat{c}}_0 \Vert {{\varvec{u}}}- {\widehat{{{\varvec{u}}}}}\Vert _{\mathcal {C}}^2 \le \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) \le {\hat{C}}_0 \Vert {{\varvec{u}}}-{\widehat{{{\varvec{u}}}}}\Vert _{\mathcal {C}}^2 \qquad \text {for all } {{\varvec{u}}}, {\widehat{{{\varvec{u}}}}}\in \mathbb {D}\subset \mathbb {V}. \end{aligned}$$
(C0)

Moreover, there exists a relative dissipation functional \(\mathcal {D}(\cdot |\cdot ): \mathbb {D}\times \mathbb {D}\rightarrow [0,\infty )\) and perturbation functionals \(\mathcal {P}(\cdot ): \mathbb {W}' \rightarrow \mathbb {R}\) and \(\mathcal {P}_\partial (\cdot ): \mathbb {W}\rightarrow \mathbb {R}\) such that

$$\begin{aligned} -\langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}}) - \mathcal {R}({\widehat{{{\varvec{u}}}}}) {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}), {{\varvec{z}}}({{\varvec{u}}})-{{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle&\le {\hat{C}}_1 \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) - 2 \mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}), \end{aligned}$$
(C1)
$$\begin{aligned} \langle \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) - \mathcal {G}({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}})\rangle&\le {\hat{C}}_2 \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}), \end{aligned}$$
(C2)
$$\begin{aligned} \langle {\widehat{{{\varvec{r}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle&\le {\hat{C}}_3 \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) + \mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) + \mathcal {P}({\widehat{{{\varvec{r}}}}}), \end{aligned}$$
(C3)
$$\begin{aligned} \langle \mathcal {B}_\partial ({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle&\le \mathcal {P}_\partial ({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), \end{aligned}$$
(C4)

for classical solutions \({{\varvec{u}}},{\widehat{{{\varvec{u}}}}}\) of (11)–(12) and (22)–(23) with positive constants \({\hat{C}}_i\), which may depend on the set \(\mathbb {D}\), the constant \({\hat{c}}_0\), \({\hat{C}}_0\), and \({\widehat{{{\varvec{u}}}}}\) and its derivatives, but not on \({{\varvec{u}}}\).

Together with the relative energy identity stated in the previous lemma, these abstract conditions immediately lead to the following stability estimate.

Lemma 6

Let Assumptions 1 and 5 hold. Then any pair of classical solutions \({{\varvec{u}}}\) and \({\widehat{{{\varvec{u}}}}}\) to the evolution Eqs. (11)–(12) and (22)–(23) satisfies

$$\begin{aligned} {\hat{c}}_0 \Vert {{\varvec{u}}}(\tau )&- {\widehat{{{\varvec{u}}}}}(\tau )\Vert _{\mathcal {C}}^2 + \int _0^\tau e^{{\hat{c}} (\tau -\sigma )}\mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) d\sigma \\&\le {\hat{C}}_0 e^{{\hat{c}} \tau } \Vert {{\varvec{u}}}(0) - {\widehat{{{\varvec{u}}}}}(0)\Vert ^2_{\mathcal {C}} \\&\qquad + \int _0^\tau e^{{\hat{c}} (\tau - \sigma )} \left[ \mathcal {P}({\widehat{{{\varvec{r}}}}}(\sigma )) + \mathcal {P}_\partial ({{\varvec{z}}}({{\varvec{u}}}(\sigma )) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}(\sigma ))) \right] \, d\sigma , \end{aligned}$$

with constant \({\hat{c}} = {\hat{C}}_1+{\hat{C}}_2 + {\hat{C}}_3\) and \({\hat{c}}_0,{\hat{C}}_0\) obtained from Assumption 5.

Proof

From Lemma 4 and Assumption 5, we immediately obtain

$$\begin{aligned} \frac{d}{d\tau }\mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) \le - \mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) + ({\hat{C}}_1 + {\hat{C}}_2 + {\hat{C}}_3) \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) + \mathcal {P}({\widehat{{{\varvec{r}}}}}) + \mathcal {P}_\partial ({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})). \end{aligned}$$

The assertion then follows by application of the Gronwall lemma [30, Ch. 29], the definition of \({\hat{c}}\), and the estimates for the relative energy in condition (C0). \(\square \)

3 Application to gas networks

We now return to the gas transport problem stated in the introduction and show that it fits into the abstract framework discussed in the previous section. In addition, we collect some auxiliary results that will be useful for the stability analysis of the next section.

3.1 Variational formulation and canonical form

We may multiply (1)–(2) by appropriate test functions q, r, integrate over the spatial domain, and use integration-by-parts in the second equation, to obtain

$$\begin{aligned} (a \partial _{\tau }\rho , q) + (\partial _xm, q)&= 0, \end{aligned}$$
(25)
$$\begin{aligned} (\varepsilon ^2 \partial _{\tau }w, r) - (h, \partial _xr) + (\gamma \frac{|w|}{a\rho } m, r)&= - h r\big |_{0}^\ell , \end{aligned}$$
(26)

where \((a,b):= \int _0^\ell a(x) b(x) dx\) is used to abbreviate the \(L^2\)-scalar product. Note that boundary conditions (6) for h could be incorporated naturally in the last term of (26). The variational identities (25)–(26) hold for all time \(t>0\) of relevance and for all smooth test functions q, r, independent of time, and they are satisfied, in particular, by all smooth solutions of (1)–(3). In compact notation, the system (25)–(26) can be stated as

$$\begin{aligned} \langle \mathcal {C}\partial _{\tau }{{\varvec{u}}}, {{\varvec{w}}}\rangle + \langle \mathcal {J}{{\varvec{z}}}({{\varvec{u}}}), {{\varvec{w}}}\rangle + \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}}), {{\varvec{w}}}\rangle&= \langle \mathcal {B}_\partial {{\varvec{z}}}, {{\varvec{w}}}\rangle , \end{aligned}$$
(27)

with state variable \({{\varvec{u}}}=(\rho ,w)\), co-state variable \({{\varvec{z}}}({{\varvec{u}}}) = (h,m)\) defined by (3), time independent test function \({{\varvec{v}}}=(q,r)\), and operators \(\mathcal {C}\), \(\mathcal {J}\), \(\mathcal {R}({{\varvec{u}}})\), and \(\mathcal {B}_\partial \) given by

$$\begin{aligned} \langle \mathcal {C}\partial _{\tau }{{\varvec{u}}}, {{\varvec{v}}}\rangle&= (a \partial _{\tau }\rho , q) + (\varepsilon ^2 \partial _{\tau }w, r),&\langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}, {{\varvec{v}}}\rangle&= (\gamma \frac{|w|}{a\rho } m, r) \\ \langle \mathcal {J}{{\varvec{z}}}, {{\varvec{v}}}\rangle&= (\partial _xm,q) - (h, \partial _xr),&\langle \mathcal {B}_\partial {{\varvec{z}}}, {{\varvec{v}}}\rangle&= - h r\big |_0^\ell . \end{aligned}$$

These definitions hold for all appropriate arguments \({{\varvec{u}}}=(\rho ,w)\), \({{\varvec{z}}}=(h,m)\), and \({{\varvec{v}}}=(q,r)\). From these variational characterizations, it is not difficult to see that \(\mathcal {J}\) is anti-symmetric and that \(\mathcal {C}\) and \(\mathcal {R}({{\varvec{u}}})\) are symmetric and at least positive semi-definite, if the parameters a, \(\varepsilon ^2\), \(\gamma \) and the density \(\rho \) are positive. The operator \(\mathcal {B}_\partial \) associated with the boundary terms does not have a particular property, except being supported only at the boundary.

Remark 7

Equation (27) is the variational form of an abstract evolution problem (11) with state and co-state variables \({{\varvec{u}}}\) and \({{\varvec{z}}}({{\varvec{u}}})\), and energy functional \(\mathcal {H}({{\varvec{u}}})=\mathcal {H}(\rho ,w)\). Problem (1)–(3) thus is a dissipative Hamiltonian system of the form (11)–(12).

3.2 Auxiliary results

A quick inspection of the above derivations shows that the operators \(\mathcal {C}\), \(\mathcal {R}({{\varvec{u}}})\), and \(\mathcal {J}\), and \(\mathcal {B}_\partial \), can be formally identified with

$$\begin{aligned} \mathcal {C}= \begin{pmatrix} a &{} 0 \\ 0 &{} \varepsilon ^2 \end{pmatrix}, \qquad \mathcal {R}({{\varvec{u}}}) = \begin{pmatrix} 0 &{} 0 \\ 0 &{} \frac{\gamma |w|}{a \rho }\end{pmatrix}, \qquad \text {and} \qquad \mathcal {J}-\mathcal {B}_\partial = \begin{pmatrix} 0 &{} \partial _x\\ \partial _x&{} 0 \end{pmatrix}. \end{aligned}$$

The latter follows rigorously by reversing the order of arguments in the derivations of the weak formulation. The energy of the system is here given by

$$\begin{aligned} \mathcal {H}({{\varvec{u}}}) = \int _0^\ell a (\varepsilon ^2 \rho \frac{w^2}{2} + P(\rho ) + \rho g z ) dx, \end{aligned}$$

and by elementary computations, we obtain the formulas

$$\begin{aligned} {{\varvec{z}}}({{\varvec{u}}}) = \begin{pmatrix} \varepsilon ^2 \frac{w^2}{2} + P'(\rho ) \\ a \rho w \end{pmatrix} \qquad \text {and} \qquad \mathcal {G}({{\varvec{u}}}) = \begin{pmatrix} P''(\rho ) &{} \varepsilon ^2 w \\ a w &{} a \rho \end{pmatrix} \end{aligned}$$

for the gradient \({{\varvec{z}}}({{\varvec{u}}}) = \mathcal {C}^{-1} \mathcal {H}'({{\varvec{u}}})\) and the Hessian \(\mathcal {G}({{\varvec{u}}}) = \mathcal {C}^{-1} \mathcal {H}''({{\varvec{u}}})\) of the energy functional.

Remark 8

Let us emphasize that the operators \(\mathcal {C}\) and \(\mathcal {R}(\cdot )\), the energy functional \(\mathcal {H}(\cdot )\), and thus the functions \({{\varvec{z}}}(\cdot )\), \(\mathcal {G}(\cdot )\) explicitly depend on the model parameters \(\varepsilon \) and \(\gamma \).

3.3 Functional analytic setting

As a next step, we briefly discuss the choice of suitable function spaces for the gas transport problem under consideration. We define

$$\begin{aligned} \mathbb {V}:= L^2(0,\ell ) \times L^2(0,\ell ) \qquad \text {and} \qquad \mathbb {W}:=H^1(0,\ell ) \times H^1(0,\ell ), \end{aligned}$$
(28)

where \(H^1(0,\ell )\) denotes the standard Sobolev space on the interval \((0,\ell )\). We further introduce the set of admissible states

$$\begin{aligned} \mathbb {D}= \{(\rho ,w) \in \mathbb {V}: (A1)-(A2) \text { are satisfied}\}. \end{aligned}$$

Recall that classical solutions of (1)–(3) satisfy \({{\varvec{u}}}=(\rho ,w) \in C^1([0,T];\mathbb {V})\cap C^0([0,T];\mathbb {D})\) and \({{\varvec{z}}}({{\varvec{u}}}) = (h,m) \in C^0([0,T];\mathbb {W})\). Then \(\mathcal {C}\) and \(\mathcal {R}({{\varvec{u}}})\) with \({{\varvec{u}}}\in \mathbb {D}\) can be understood as self-adjoint and positive semi-definite bounded linear operators mapping from \(\mathbb {V}\) or \(\mathbb {W}\) to the dual spaces \(\mathbb {V}'\) or \(\mathbb {W}'\). For \(\varepsilon \), a uniformly positive and bounded, \(\mathcal {C}\) induces a norm

$$\begin{aligned} \Vert {{\varvec{u}}}\Vert _\mathcal {C}^2 = \Vert \sqrt{a} \rho \Vert ^2_{L^2(0,\ell )} + \Vert \varepsilon w\Vert ^2_{L^2(0,\ell )}, \end{aligned}$$
(29)

which is equivalent to the standard norm on \(\mathbb {V}= L^2(0,\ell ) \times L^2(0,\ell )\). Adopting the previous notation, we write \(\langle \cdot , \cdot \rangle \) for the duality products on \(\mathbb {V}' \times \mathbb {V}\) and \(\mathbb {W}' \times \mathbb {W}\), respectively, and use \((a,b)=\int _0^\ell a(x) b(x) \, dx\) to denote the scalar product of \(L^2(0,\ell )\). From the variational definition of the operator \(\mathcal {J}\), one can see that

$$\begin{aligned} \langle \mathcal {J}{{\varvec{w}}}, {\widetilde{{{\varvec{w}}}}} \rangle :=&(\partial _xm, {\tilde{h}}) - (h, \partial _x{\tilde{m}}) = -\langle \mathcal {J}{\widetilde{{{\varvec{w}}}}}, {{\varvec{w}}}\rangle , \end{aligned}$$
(30)

for all \({{\varvec{w}}}=(h,m)\), \({\widetilde{{{\varvec{w}}}}}=({\tilde{h}}, {\tilde{m}}) \in \mathbb {W}\), i.e., \(\mathcal {J}\) is skew-symmetric on \(\mathbb {W}\). The formula \(\langle \mathcal {B}_\partial {{\varvec{z}}}({{\varvec{u}}}), {{\varvec{v}}}\rangle = - h r\big |_0^\ell \) with \({{\varvec{v}}}=(q,r)\) finally shows that the boundary operator \(\mathcal {B}_\partial \) acts on the co-state variables. Note that the required boundary values are well-defined for functions \({{\varvec{z}}}({{\varvec{u}}}) = (h,m)\) and \({{\varvec{v}}}=(q,r) \in \mathbb {W}= H^1(0,\ell ) \times H^1(0,\ell )\).

In summary, we thus have shown that (1)–(3) can be interpreted as an abstract port-Hamiltonian system (11)–(12), and we verified that under conditions (A1)–(A2) also Assumption 1 is valid.

4 Stability analysis and parabolic limit

We will now verify the conditions of Assumption 5 for the gas transport problem (1)–(3) and then utilize the abstract stability results to prove stability of solutions with respect to perturbations in the model parameters as well as initial and boundary values. As a case of particular interest, we will study convergence in the parabolic limit \(\varepsilon \rightarrow 0\).

4.1 Perturbed problem

Let \(({\hat{\rho }},{\hat{w}})\) denote a solution to the perturbed equations

$$\begin{aligned} a \partial _{\tau }{\hat{\rho }} + \partial _x(a {\hat{\rho }} {\hat{w}})&= 0, \end{aligned}$$
(31)
$$\begin{aligned} {\hat{\varepsilon }}^2 \partial _{\tau }{\hat{w}} + \partial _x( P'({\hat{\rho }}) + {\hat{\varepsilon }}^2 \frac{{\hat{w}}^2}{2} + g z) + {\hat{\gamma }} \frac{|{\hat{w}}|}{ a\hat{\rho }} {\hat{m}}&= 0, \end{aligned}$$
(32)

which are again assumed to hold for all \(0< x < \ell \) and time \(t>0\). The terms carrying spatial derivatives are the corresponding co-state variables

$$\begin{aligned} {\hat{h}}({\hat{\rho }},{\hat{v}}) = {\hat{\varepsilon }}^2\frac{{\hat{w}}^2}{2} + P'(\hat{\rho }) + g z \qquad \text {and} \qquad {\hat{m}}({\hat{\rho }},{\hat{v}}) = a {\hat{\rho }} {\hat{w}}, \end{aligned}$$
(33)

which are again directly related to the derivatives of the corresponding energy functional \({\hat{\mathcal {H}}}({\hat{\rho }},{\hat{w}}) = \int _0^\ell a ({\hat{\varepsilon }}^2 \frac{{\hat{w}}^2}{2{\hat{\rho }}} + P(\hat{\rho }) + {\hat{\rho }} g z ) \, dx\). By the some elementary manipulations and the same arguments as employed above, this system can again be written in the abstract form

$$\begin{aligned} \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}}+ (\mathcal {J}+ \mathcal {R}({\widehat{{{\varvec{u}}}}})) {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})&= \mathcal {B}_\partial {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) + {\widehat{{{\varvec{r}}}}}\end{aligned}$$
(34)
$$\begin{aligned} {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})&= \mathcal {C}^{-1} \mathcal {H}'({\widehat{{{\varvec{u}}}}}), \end{aligned}$$
(35)

with functionals \(\mathcal {H}(\cdot )\) and \({{\varvec{z}}}(\cdot )\), and operators \(\mathcal {C}\), \(\mathcal {J}\), and \(\mathcal {R}(\cdot )\) denoting those for the unperturbed problem with parameters \(\varepsilon \) and \(\gamma \), and residual \({\widehat{{{\varvec{r}}}}}=({\widehat{{{\varvec{r}}}}}_1,{\widehat{{{\varvec{r}}}}}_2)\) given by

$$\begin{aligned} {\widehat{{{\varvec{r}}}}}_1 = 0 \qquad \text {and} \qquad {\widehat{{{\varvec{r}}}}}_2 = (\varepsilon ^2- {\hat{\varepsilon }}^2) (\partial _{\tau }{\hat{w}} + \tfrac{1}{2} \partial _x|{\hat{w}}|^2) + (\gamma - {\hat{\gamma }}) |{\hat{w}}| {\hat{w}}. \end{aligned}$$
(36)

In the following section, we verify the abstract assumptions required for our stability analysis, without explicitly taking into account the special form of \({\widehat{{{\varvec{r}}}}}_1\) and \({\widehat{{{\varvec{r}}}}}_2\).

4.2 Verification of conditions (C0)–(C4)

We always assume in the following that assumptions (A1)–(A2) are valid. Constants arising in the estimates may depend on the bounds of these conditions. Since we later consider the case \(\varepsilon \rightarrow 0\), we will make explicit the dependence of the constants on this scaling parameter.

4.3 Condition (C0)

From Taylor’s formula, we know that

$$\begin{aligned} f({{\varvec{u}}}|{\hat{{{\varvec{u}}}}})&= f({{\varvec{u}}}) - f({\widehat{{{\varvec{u}}}}}) - f'({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}-{\widehat{{{\varvec{u}}}}}) \\&= \frac{1}{2} \int _0^1 \langle (1-s) f''({\widehat{{{\varvec{u}}}}}+ s ({{\varvec{u}}}-{\widehat{{{\varvec{u}}}}})) ({{\varvec{u}}}-{\widehat{{{\varvec{u}}}}}), {{\varvec{u}}}-{\widehat{{{\varvec{u}}}}}) \rangle ds. \end{aligned}$$

Now let \(f({{\varvec{u}}}) = a(\varepsilon ^2 \rho \frac{|w|^2}{2} + P(\rho ))\) denote the integrand of the energy functional defined above, and note that its Hessian is given by

$$\begin{aligned} f''(\rho ,w) = a \begin{pmatrix} P''(\rho ) &{} \quad \varepsilon ^2 w \\ \varepsilon ^2 w &{} \quad \varepsilon ^2 \rho \end{pmatrix}. \end{aligned}$$

By multiplying with \({{\varvec{v}}}=(x,y)\) from the left and right, one can see that

$$\begin{aligned} \frac{1}{a} \langle f''(\rho ,w) {{\varvec{v}}}, {{\varvec{v}}}\rangle&= P''(\rho ) x^2 + 2 \varepsilon ^2 w x y + \varepsilon ^2 \rho y^2 \\&= P''(\rho ) x^2 + 2 \varepsilon w x (\varepsilon y) + \rho (\varepsilon y)^2. \end{aligned}$$

The second term in the second line can be estimated by

$$\begin{aligned} |2 \varepsilon w x (\varepsilon y) |&\le 2 \frac{\varepsilon ^2 w^2}{\rho } x^2 + \frac{1}{2} \rho (\varepsilon y)^2. \end{aligned}$$

From condition (7), one can see that \(2 \frac{\varepsilon ^2 w^2}{\rho } \le P''(\rho )/2\), which leads to the lower bound

$$\begin{aligned} \langle f''(\rho ,w) {{\varvec{v}}}, {{\varvec{v}}}\rangle \ge \frac{a}{2} \left( P''(\rho ) x^2 + \rho |\varepsilon y|^2 \right) , \end{aligned}$$

and the corresponding upper bound

$$\begin{aligned} \langle f''(\rho ,w) {{\varvec{v}}}, {{\varvec{v}}}\rangle \le \frac{3a}{2} \left( P''(\rho ) x^2 + \rho | \varepsilon y|^2 \right) . \end{aligned}$$

Using the uniform bounds for \(\rho \), \(P''(\rho )\), and a, we arrive at the following result.

Lemma 9

Let (A1)–(A2) hold. Then

$$\begin{aligned} {\hat{c}}_0 \Vert {{\varvec{u}}}- {\widehat{{{\varvec{u}}}}}\Vert _\mathcal {C}^2 \le \mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) \le {\hat{C}}_0 \Vert {{\varvec{u}}}- {\widehat{{{\varvec{u}}}}}\Vert _\mathcal {C}^2 \end{aligned}$$

with positive constants \({\hat{c}}_0,{\hat{C}}_0\) only depending on the bounds on \(\rho \) and \(P''(\rho )\) in the assumptions.

Remark 10

Lemma 9 is strongly based on the fact that \(\mathcal {H}\) in (4) is a strictly convex functional of \((\rho ,w)\) on the set of subsonic bounded states. Let us mention that \(\mathcal {H}\) can also be viewed as a strictly convex functional of the conservative variables \((\rho ,\rho w)\). This is the viewpoint taken in [12, 19].

4.4 Condition (C1)

Using \({{\varvec{u}}}=(\rho ,w)\) and \({\widehat{{{\varvec{u}}}}}=({\hat{\rho }},{\hat{w}})\), and the definition of \(\mathcal {R}({{\varvec{u}}})\) for our particular problem, the left hand side of (C1) can be expressed as

$$\begin{aligned} -\langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}({{\varvec{u}}})&- \mathcal {R}({\widehat{{{\varvec{u}}}}}) {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}), {{\varvec{z}}}({{\varvec{u}}}) -{{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle \\&= - \int _0^\ell \gamma \, (|w| w - |{\hat{w}}| {\hat{w}}) \, a \, (\rho w - {\hat{\rho }} {\hat{w}}) \, dx =: (*). \end{aligned}$$

The first term in the integrand can be written as \(|w| w - |{\hat{w}}| {\hat{w}} = 2 \int _0^1 |{\hat{w}}+s(w-{\hat{w}})| ds (w - {\hat{w}})\), and one can observe that \(\frac{|w|+|{\hat{w}}|}{4} \le \int _0^1 |{\hat{w}} + s (w-{\hat{w}})| ds \le \frac{|w|+|{\hat{w}}|}{2}\). The second term in the integrand can be expanded as \(\rho w - {\hat{\rho }} {\hat{w}} = {\hat{\rho }} (w - {\hat{w}}) + (\rho - {\hat{\rho }}) w \), which leads to

$$\begin{aligned} (|w| w - |{\hat{w}}| {\hat{w}}) \, (\rho w - {\hat{\rho }} {\hat{w}})&\ge {\hat{\rho }} \frac{|w|+|{\hat{w}}|}{2} |w-{\hat{w}}|^2 - |w-{\hat{w}}| |w+{\hat{w}}| |w| |\rho - {\hat{\rho }}| \\&\ge {\hat{\rho }} \frac{|w|+|{\hat{w}}|}{4} |w - {\hat{w}}|^2 - (|w|+|{\hat{w}}|) \frac{|w|^2}{{\hat{\rho }}} |\rho - {\hat{\rho }}|^2. \end{aligned}$$

In summary, we thus arrive at the estimate

$$\begin{aligned} (*) \le - \frac{1}{4} \int _0^\ell \gamma a {\hat{\rho }} (|w|+|{\hat{w}}|) |w-{\hat{w}}|^2 dx + \int _0^\ell \gamma \frac{|w|^2}{{\hat{\rho }}} (|w| + |{\hat{w}}|) a |\rho - {\hat{\rho }}|^2 dx. \end{aligned}$$

The uniform bounds for \({\hat{\rho }}\), w, \({\hat{w}}\), \(\gamma \), a and condition (C0) then lead to the following result.

Lemma 11

Let assumptions (A1)–(A2) be valid. Then condition (C1) holds with

$$\begin{aligned} \mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) = \frac{1}{8} \int _0^\ell \gamma a {\hat{\rho }} (|w| + |{\hat{w}}|) (w- {\hat{w}})^2 ds, \end{aligned}$$
(37)

where \({{\varvec{u}}}=(\rho ,w)\) and \({\widehat{{{\varvec{u}}}}}=({\hat{\rho }},{\hat{w}})\), and with constant \({\hat{C}}_1=2 {\bar{\gamma }} |{\bar{w}}|^2 \underline{\rho }^{-1} {\hat{c}}_0^{-1}\).

Remark 12

By the elementary fact that \(|w|+|{\hat{w}}| \ge |w-{\hat{w}}|\), one can see that

$$\begin{aligned} \mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}}) \ge c_D \Vert w-{\hat{w}}\Vert _{L^3}^3, \end{aligned}$$
(38)

with positive constant \(c_D=\frac{\underline{\gamma }\underline{a} \underline{\rho }}{8}\). Thus, the relative dissipation \(\mathcal {D}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}})\) provides control over the velocity perturbation even in the case \(\varepsilon \rightarrow 0\), where the velocity contribution to the relative energy \(\mathcal {H}({{\varvec{u}}}|{\widehat{{{\varvec{u}}}}})\) disappears. This will later be used in our stability analysis.

4.5 Condition (C2)

Let us start with noting that

$$\begin{aligned} {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) - \mathcal {G}({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}}) = \begin{pmatrix} P'(\rho |{\hat{\rho }}) + \varepsilon ^2 (w - {\hat{w}})^2/2\\ a (\rho - {\widehat{\rho }}) (w - {\widehat{w}}) \end{pmatrix}, \end{aligned}$$
(39)

which follows directly from the definitions of the gradient \({{\varvec{z}}}({{\varvec{u}}})\) and the Hessian \(\mathcal {G}({{\varvec{u}}})\) for the problem under investigation. By assumption (A1)–(A2), the pressure potential P is smooth and \(\rho ,{\hat{\rho }}\) are uniformly bounded, and consequently

$$\begin{aligned} P'(\rho |{\hat{\rho }}) \le C a |\rho - {\hat{\rho }}|^2, \end{aligned}$$

with some constant C only depending on the bounds of the coefficients, the density, and the pressure potential. The second line in (39) can be estimated by

$$\begin{aligned} a (\rho - {\widehat{\rho }}) (w - {\widehat{w}}) \le \frac{1}{\varepsilon } ( a^2 (\rho - {\hat{\rho }})^2 + \varepsilon ^2 (w - {\hat{w}})^2 ) \end{aligned}$$

via Young’s inequality. The left hand side of (C2) can then be bounded by

$$\begin{aligned} \langle \mathcal {C}\partial _{\tau }{\widehat{{{\varvec{u}}}}},&{{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) - \mathcal {G}({\widehat{{{\varvec{u}}}}}) ({{\varvec{u}}}- {\widehat{{{\varvec{u}}}}})\rangle \\&= \int _0^\ell a \partial _{\tau }{\hat{\rho }} \big (P'(\rho |{\hat{\rho }}) + \varepsilon ^2 (w-{\hat{w}})^2/2 \big )dx + \int _0^\ell \varepsilon ^2 \partial _{\tau }{\hat{w}} a (\rho - {\hat{\rho }}) (w - {\hat{w}}) dx \\&\le (\Vert a \partial _{\tau }{\hat{\rho }}\Vert _{L^\infty } + \Vert \varepsilon \partial _{\tau }{\hat{w}}\Vert _{L^\infty }) ((C+a) |\rho - {\hat{\rho }}|^2 + \tfrac{3}{2} \varepsilon ^2 |w - {\hat{w}}|^2). \end{aligned}$$

Together with the bounds (C0), we thus obtain the following result.

Lemma 13

Let (A1)–(A2) hold. Then (C2) is valid with \({\hat{C}}_2 = C (\Vert \partial _{\tau }{\hat{\rho }}\Vert _{L^\infty } + \Vert \varepsilon \partial _{\tau }{\hat{w}}\Vert _{L^\infty })\) and constant \(C \ge 0\) depending only on the bounds in the assumptions.

4.6 Condition (C3)

We start by expanding

$$\begin{aligned} \langle {\widehat{{{\varvec{r}}}}}, {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle&= \int _0^\ell {\widehat{{{\varvec{r}}}}}_1 \left( \frac{\varepsilon ^2}{2} (|w|^2-|{\hat{w}}|^2) + P'(\rho ) - P'({\hat{\rho }}) \right) dx \\&\qquad + \int _0^\ell {\widehat{{{\varvec{r}}}}}_2 a (\rho w - {\hat{\rho }} {\hat{w}}) \, dx = (i)+(ii). \end{aligned}$$

Using the uniform bounds in assumption (A2) and smoothness of \(P(\cdot )\), we obtain

$$\begin{aligned} (i)&\le ({\bar{w}} \varepsilon ^2 \Vert w - {\hat{w}}\Vert _{L^2} + C_P'' \Vert \rho - \hat{\rho }\Vert _{L^2}) \Vert {\widehat{{{\varvec{r}}}}}_1\Vert _{L^2} \\&\le \varepsilon ^2 \Vert w-{\hat{w}}\Vert ^2_{L^2} + \frac{1}{2}\Vert \sqrt{a}(\rho - \hat{\rho })\Vert ^2_{L^2} + C_1 \Vert {\widehat{{{\varvec{r}}}}}_1\Vert ^2_{L^2}, \end{aligned}$$

where we applied Hölder’s inequality in the last step. Condition (C0) then allows to bound the first two terms by the relative energy. For the second term, we obtain

$$\begin{aligned} (ii)&= \int _0^\ell {\widehat{{{\varvec{r}}}}}_2 a ( (\rho - {\hat{\rho }}) w +{\hat{\rho }} (w-{\hat{w}})) dx \\&\le {\bar{w}} \sqrt{{\bar{a}}} \Vert {\widehat{{{\varvec{r}}}}}_2\Vert _{L^2} \Vert \sqrt{a}(\rho - \hat{\rho })\Vert _{L^2} + {\bar{a}} {\bar{\rho }} \Vert {\widehat{{{\varvec{r}}}}}_2\Vert _{L^{3/2}} \Vert (w - {\hat{w}})\Vert _{L^3} \\&\le \frac{1}{2}\Vert \sqrt{a}(\rho - {\hat{\rho }})\Vert ^2 + C_2 \Vert {\widehat{{{\varvec{r}}}}}_2\Vert _{L^2}^2 + \delta \Vert w-{\hat{w}}\Vert _{L^3}^3 + C(\delta ) \Vert {\widehat{{{\varvec{r}}}}}_2\Vert _{L^{3/2}}^{3/2}. \end{aligned}$$

Choosing \(\delta =c_D\) and using (38) allows to bound the third term in the last line by the relative dissipation. In summary, we then obtain the following result.

Lemma 14

Let (A1)–(A2) be valid. Then condition (C3) holds with \({\hat{C}}_3 = 1\) and

$$\begin{aligned} \mathcal {P}({\widehat{{{\varvec{r}}}}}) = C_1 \Vert {\widehat{{{\varvec{r}}}}}_1\Vert _{L^2}^2 + C_2 \Vert {\widehat{{{\varvec{r}}}}}_2\Vert _{L^2}^2 + C_3 \Vert {\widehat{{{\varvec{r}}}}}_2\Vert _{L^{3/2}}^{3/2}, \end{aligned}$$

with constants \(C_1\), \(C_2\), and \(C_3\) only depending on the bounds in (A1)–(A2).

Using the specific form of the residual given in (36), we may further estimate the perturbation functional \(\mathcal {P}({\widehat{{{\varvec{r}}}}})\) as follows.

Corollary 15

Let the conditions of Lemma 14 be valid and \({\widehat{{{\varvec{r}}}}}\) be defined as in (36). Then

$$\begin{aligned} \mathcal {P}({\widehat{{{\varvec{r}}}}}) \le {\hat{C}}_3' |\varepsilon ^2 - {\hat{\varepsilon }}^2|^{3/2} + {\hat{C}}_3'' |\gamma - {\hat{\gamma }}|^{3/2}, \end{aligned}$$

with constants \({\hat{C}}_3'\), \({\hat{C}}_3''\) only depending on the bounds in assumptions (A1)–(A2) as well as on \(\Vert \partial _{\tau }{\hat{w}}\Vert _{L^\infty (0,T;L^2)}\) and \(\Vert \partial _x{\hat{w}}\Vert _{L^\infty (0,T;L^2)}\).

4.7 Condition (C4)

Using the definition of the state and co-state variables, as well as the variational characterization of the boundary operator, we immediately obtain

$$\begin{aligned} \langle \mathcal {B}_\partial ({{\varvec{z}}}({{\varvec{u}}})&- {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}}) \rangle \\&= - (h(\rho ,w) - h({\hat{\rho }},{\hat{w}})) (m(\rho ,w) - m({\hat{\rho }},{\hat{w}}))\big |_0^\ell \\&\le |h(\rho ,w) - h({\hat{\rho }},{\hat{w}})|_\partial |m(\rho ,w) - m(\hat{\rho },{\hat{w}})|_\partial . \end{aligned}$$

with \(h(\tilde{\rho },{\tilde{w}}) = \varepsilon ^2 \frac{{\tilde{w}}^2}{2} + P'(\tilde{\rho })\) and \(m(\tilde{\rho },{\tilde{w}})=a\tilde{\rho }{\tilde{w}}\) denoting the co-state mappings of the unperturbed problem, and \(|a|_\partial = \sqrt{|a(0)|^2 + |a(\ell )|^2}\) denoting the \(\ell ^2\)-scalar product on the space of boundary values. We thus obtain

Lemma 16

Let assumptions (A1)–(A2) hold. Then condition (C4) is valid with

$$\begin{aligned} \mathcal {P}_\partial ({{\varvec{z}}}({{\varvec{u}}})-{{\varvec{z}}}({\widehat{{{\varvec{u}}}}}))&= |h(\rho ,w) - h({\hat{\rho }},{\hat{w}})|_\partial |m(\rho ,w) - m(\hat{\rho },{\hat{w}})|_\partial . \end{aligned}$$

While \(h(\rho ,w)\) and \(m(\rho ,w)\) amount to the natural boundary values of the unperturbed problem, the evaluation \(h({\hat{\rho }},{\hat{w}})\) and \(m({\hat{\rho }},{\hat{v}})\) of the unperturbed co-state mappings at the perturbed states does not have a physical meaning. We therefore decompose

$$\begin{aligned} h({\hat{\rho }},{\hat{v}})&= {\hat{h}}({\hat{\rho }},{\hat{v}}) + (h({\hat{\rho }},{\hat{v}}) - {\hat{h}}({\hat{\rho }},{\hat{v}})) \end{aligned}$$

into the natural boundary value \({\hat{h}}({\hat{\rho }},{\hat{w}}) = \hat{\varepsilon }^2 \frac{{\hat{w}}^2}{2} + P'({\hat{\rho }})\) and a corresponding perturbation \(h({\hat{\rho }},{\hat{v}}) - {\hat{h}}({\hat{\rho }},{\hat{v}})\) of the state to co-state mapping. The latter can be estimated by the bounds in assumptions (A1)–(A2), leading to the following result.

Corollary 17

Let the assumptions of Lemma 16 hold and \(h_\partial =h(\rho ,v)|_{\{0,\ell \}}\) and \({\hat{h}}_\partial ={\hat{h}}({\hat{\rho }},{\hat{v}})|_{\{0,\ell \}}\) denote the boundary values of the co-state variables of the unperturbed and perturbed system, respectively. Then condition (C4) holds with

$$\begin{aligned} \mathcal {P}_\partial ({{\varvec{z}}}({{\varvec{u}}})-{{\varvec{z}}}({\widehat{{{\varvec{u}}}}}))&= {\hat{C}}_\partial \left( |h_\partial - {\hat{h}}_{\partial }| + |\varepsilon ^2 - {\hat{\varepsilon }}^2| \right) \end{aligned}$$

with constant \({\hat{C}}_\partial \) depending only on the bounds in (A1)–(A2) and the Lipschitz bounds for \(({\hat{\rho }},{\hat{w}})\).

4.8 Stability estimate

Having verified the conditions of our abstract stability analysis, we can now apply Lemma 6 to obtain the following stability estimate.

Theorem 18

Let assumptions (A1)–(A2) hold and let \((\rho ,w)\) and \(({\hat{\rho }}, {\hat{w}})\) denote corresponding classical solutions of (1)–(3) for parameters \(\varepsilon ,\gamma \) and \({\hat{\varepsilon }},{\hat{\gamma }}\), respectively. Further assume that \((\hat{\rho },{\hat{w}})\) is Lipschitz continuous. Then

$$\begin{aligned}&\Vert \rho (\tau ) - {\hat{\rho }}(\tau )\Vert _{L^2(0,\ell )}^2 + \varepsilon ^2 \Vert w(\tau ) - {\hat{w}}(\tau )\Vert _{L^2(0,\ell )}^2 + \int _0^\tau \Vert w(s) - {\hat{w}}(s)\Vert _{L^3(0,\ell )}^3 ds\\&\quad \le {\hat{C}} e^{{\hat{c}} \tau } \big (\Vert \rho (0) - {\hat{\rho }}(0)\Vert _{L^2(0,\ell )}^2 + \varepsilon ^2 \Vert w(0) - {\hat{w}}(0)\Vert _{L^2(0,\ell )}^2 \\&\qquad + |\gamma - {\hat{\gamma }}|^{3/2} + |\varepsilon ^2 - {\hat{\varepsilon }}^2| + \int _0^\tau |h_\partial (s) - {\hat{h}}_\partial (s)|_\partial ds \big ), \end{aligned}$$

where \(h_\partial = h(\rho ,w)|_{\{0,\ell \}}\) and \({\hat{h}}_\partial = {\hat{h}}({\hat{\rho }},{\hat{w}})|_{\{0,\ell \}}\) are the boundary values of the corresponding co-state variables. Moreover, the constants \({\hat{c}}\), \({\hat{C}}\) in this estimate only depend on the bounds in assumptions (A1)–(A2) and the Lipschitz bounds for \(({\hat{\rho }},{\hat{w}})\).

Remark 19

Let us briefly discuss the conditions and conclusions of the theorem: The additional regularity for the reference solution \((\hat{\rho },{\hat{w}})\) was already required in Lemma 13 and Corollary 15. The reduction in the convergence rate with respect to \(\gamma \), on the other hand, comes from the fact, that the error in w is estimated by the friction term rather than the kinetic energy, in order to obtain estimates that are uniform in \(\varepsilon \). The further reduction of the convergence order in \(\varepsilon \) is due to effects from the boundary. With minor modifications of the proofs, one could handle further perturbations, e.g., in the pressure potential \(P(\cdot )\), the cross-section a, or the height function z, and also deal with other boundary conditions. A careful inspection of the proofs would also allow to relax the smoothness assumptions on \((\rho ,w)\) and \(({\hat{\rho }},{\hat{w}})\) to some extent.

4.9 The parabolic limit problem

We now study convergence of solutions to (1)–(3) in the limit \(\varepsilon \rightarrow 0\). For \({\hat{\varepsilon }}=0\) and \({\hat{\gamma }}=\gamma \), the resulting limit problem reads

$$\begin{aligned} a \partial _t{\hat{\rho }} + \partial _x(a {\hat{\rho }} {\hat{w}})&= 0, \end{aligned}$$
(40)
$$\begin{aligned} \partial _x(P'({\hat{\rho }})) + \gamma |{\hat{w}}| {\hat{w}}&= 0. \end{aligned}$$
(41)

This is a degenerate parabolic problem, whose solvability can be deduced from the results in [1, 25]; a detailed analysis can be found in [27]. A formal application of Theorem 18 to this limiting case directly leads to the following result.

Theorem 20

Let (A1)–(A2) hold and let \((\rho ,w)\) and \(({\hat{\rho }},{\hat{w}})\) denote classical solutions of (1)–(3) and (40)–(41), respectively. Further assume that \(\rho ={\hat{\rho }}\) at time \(t=0\) and \(\varepsilon ^2 \frac{w^2}{2} + P'(\rho ) = P'({\hat{\rho }})\) at the boundary \(x\in \{0,\ell \}\) for all \(0 \le t \le T\). In addition, assume that \(({\hat{\rho }},{\hat{w}})\) is Lipschitz continuous. Then

$$\begin{aligned} \Vert \rho (\tau ) - {\hat{\rho }}(\tau )\Vert _{L^2(0,\ell )}^2 + \int _0^\tau \Vert w(s) - {\hat{w}}(s)\Vert _{L^3(0,\ell )}^3 ds&\le {\hat{C}} e^{{\hat{c}} \tau } \varepsilon ^2 , \end{aligned}$$

with constants \({\hat{c}},{\hat{C}}\) having the same properties as in Theorem 18.

Remark 21

On bounded time intervals, the quadratic norm difference between solutions of the hyperbolic problem (1)–(3) and the parabolic limit problem (40)–(41) can thus bounded at least by \(O(\varepsilon ^2)\). In simulations for a simple test case with length \(\ell =5\), cross section \(a=1\), friction parameter \(\gamma =1\), pressure law \(p(\rho )=\rho ^2/2\), boundary data \(h_\partial (0,t)=1+t^2\), \(h_\partial (\ell ,t)=1\), and maximal time \(T=5\), we observed

\(\varepsilon \)

\(2^{-1}\)

\(2^{-2}\)

\(2^{-3}\)

\(2^{-4}\)

\(2^{-5}\)

\(2^{-6}\)

\(\Vert \rho - {\hat{\rho }}\Vert _{L^\infty (L^2)}^2\)

4.25e-0

6.85e-1

6.06e-2

4.20e-3

2.70e-4

1.70e-5

Rate

–

2.63

3.50

3.85

3.96

3.99

Note that a formal asymptotic analysis would predict a rate \(O(\varepsilon ^4)\) and such a rate has actually been proven in [18] for linear friction and unbounded domains. A careful revision of our arguments and additional assumptions on the solution may allow to rigorously establish these improved rates also for the more general problems considered here.

5 Extension to gas networks

As a final step of our analysis, we now show that, in the spirit of port-Hamiltonian modelling [28, 29], the results of the previous section can be extended almost verbatim to gas networks. We start with extending the rescaled model (1)–(6) to gas networks and then discuss the modifications needed in the stability and asymptotic analysis.

5.1 Network topology

Let \((\mathcal {V},\mathcal {E})\) denote a directed and connected finite graph with vertices \(v \in \mathcal {V}\) and edges \(e \in \mathcal {E}\), which are identified with intervals \((0,\ell ^e)\). We denote by \(\mathcal {E}(v)\) the set of edges incident to the vertex v, and decompose \(\mathcal {V}=\mathcal {V}_0 \cup \mathcal {V}_\partial \) into the sets of interior and boundary vertices, characterized by \(\mathcal {V}_0=\{v \in \mathcal {V}: |\mathcal {E}(v)|>1\}\) and \(\mathcal {V}_\partial =\{v \in \mathcal {V}: |\mathcal {E}(v)|=1\}\). Here \(|\mathcal {E}(v)|\) denotes the cardinality of the set \(\mathcal {E}(v)\). We further associate to any vertex \(v \in \mathcal {V}\) and edge \(e \in \mathcal {E}(v)\) a number

$$\begin{aligned} n^e(v) = {\left\{ \begin{array}{ll} 1 &{} \text {if } e=(\cdot ,v), \\ -1 &{} \text {if } e=(v,\cdot ). \end{array}\right. } \end{aligned}$$

The vertex v thus corresponds to the end point \(\ell ^e\) or the start point 0 of the interval \((0,\ell ^e)\) representing the edge e.

5.2 Gas transport on networks

After rescaling as outlined in the introduction, also see Appendix 1, the gas transport on every edge of the network is described by

$$\begin{aligned} a^e \partial _{\tau }\rho ^e + \partial _xm^e&= 0 , \qquad e \in \mathcal {E} \end{aligned}$$
(42)
$$\begin{aligned} \varepsilon ^2 \partial _{\tau }w^e + \partial _xh^e + \gamma ^e |w^e| w^e&= 0, \qquad e \in \mathcal {E}. \end{aligned}$$
(43)

By superscript \(^e\) we here denote functions or parameters restricted to the edge e. The corresponding co-state variables are defined accordingly by

$$\begin{aligned} h^e&= \varepsilon ^2 \frac{|w^e|^2}{2} + P'(\rho ^e) + g z^e, \qquad e \in \mathcal {E}, \end{aligned}$$
(44)
$$\begin{aligned} m^e&= a^e \rho ^e w^e, \,~~~~~~\qquad \qquad \qquad \qquad \quad \! e \in \mathcal {E}. \end{aligned}$$
(45)

These equations, which correspond to the conservation of mass and the balance of momentum, describe the flow of gas in the individual pipes. As outlined in [9, 26], the coupling across pipe junctions can be modeled by the following conditions

$$\begin{aligned} \sum _{e \in \mathcal {E}(v)} m^e(v) n^e(v)&= 0, ~~\qquad v \in \mathcal {V}_0, \end{aligned}$$
(46)
$$\begin{aligned} h^e(v)&= h^v, \qquad e \in \mathcal {E}(v), \ v \in \mathcal {V}_0, \end{aligned}$$
(47)

which correspond to conservation of mass and continuity of the total specific enthalpy h at pipe junctions. Note that \(h^v\) thus corresponds to the unique value of the enthalpy at the junction \(v \in \mathcal {V}_0\). A combination of the two conditions allows to show that no energy is produced via flow over junctions [9, 26]. Similar to the case of a single pipe, we may again prescribe the enthalpy at the boundary vertices by

$$\begin{aligned} h^e(v) = h_\partial ^v, \qquad v \in \mathcal {V}_\partial . \end{aligned}$$
(48)

5.3 Weak formulation and canonical form

In a similar manner to Sect. 3, we multiply (42)–(43) with suitable test functions, integrate over the edges, use integration-by-parts, and then sum over all edges, which immediately leads to the variational equations

$$\begin{aligned} \sum _{e \in \mathcal {E}} (a^e \partial _{\tau }\rho ^e, q^e)_e + (\partial _xm^e, q^e)_e&= 0, \\ \sum _{e \in \mathcal {E}} (\varepsilon ^2 \partial _{\tau }w^e, r^e)_e - (h^e, \partial _xr^e)_e + (\gamma ^e \frac{|w^e|}{a^e \rho ^e} m^e, r^e)_e&{=} -\sum _{v \in \mathcal {V}} \sum _{e \in \mathcal {E}(v)} h^e(v) r^e(v) n^e(v), \end{aligned}$$

which hold for all time independent piecewise regular test functions q, r, and all \(t>0\) of relevance. In the last term, we used the elementary identity

$$\begin{aligned} \sum _{e \in \mathcal {E}} \sum _{v \in e} h^e(v) r^e(v) n^e(v)&= \sum _{v \in \mathcal {V}} \sum _{e \in \mathcal {E}(v)} h^e(v) r^e(v) n^e(v) \end{aligned}$$

to change the order of summation. In summary, we see that the weak formulation of (42)–(45) again amounts to an abstract port-Hamiltonian system (11)–(12) with operators

$$\begin{aligned} \langle \mathcal {C}{{\varvec{u}}}, {{\varvec{v}}}\rangle&= \sum _e (a^e \rho ^r, q^e)_e + (\varepsilon ^2 w^e, r^e)_e \\ \langle \mathcal {R}({{\varvec{u}}}) {{\varvec{z}}}, {{\varvec{v}}}\rangle&= \sum _e (\gamma ^e \frac{|w^e|}{a^e \rho ^e} m^e, r^e)_e, \\ \langle \mathcal {J}{{\varvec{z}}}, {{\varvec{v}}}\rangle&= \sum _e (\partial _xm^e, q^e)_e - (h^e, \partial _xr^e)_e, \end{aligned}$$

and boundary operator defined by

$$\begin{aligned} \langle \mathcal {B}_\partial {{\varvec{z}}}, {{\varvec{v}}}\rangle&= -\sum _{v \in \mathcal {V}} \sum _{e \in \mathcal {E}(v)} h^e(v) r^e(v) n^e(v). \end{aligned}$$

Let us note that the coupling and boundary conditions (46)–(48) have not been incorporated up to this point. At this moment, the above variational equations thus describe the gas transport in a collection of separated pipes. The corresponding function spaces \(\mathbb {V}\) and \(\mathbb {W}\) for the network are thus simply given by

$$\begin{aligned} \mathbb {V}= \prod _{e \in \mathcal {E}} \mathbb {V}^e \qquad \text {and} \qquad \mathbb {W}= \prod _{e \in \mathcal {E}} \mathbb {W}^e \end{aligned}$$

with \(\mathbb {V}^e = L^2(0,\ell ^e) \times L^2(0,\ell ^e)\) and \(\mathbb {W}^e= H^1(0,\ell ^e) \times H^1(0,\ell ^e)\) denoting the corresponding spaces for the individual pipes.

5.4 Verification of conditions (C0)–(C4)

By the simple additive construction, conditions (C0)–(C3) can be obtained with the same arguments as for a single edge \(e \in \mathcal {E}\) and summation over all edges. It thus remains to consider condition (C4) in detail. We will assume that assumptions (A1)–(A2) hold uniformly for all edges \(e \in \mathcal {E}\) in the following.

We start with splitting the boundary term via

$$\begin{aligned}&\langle \mathcal {B}_\partial ({{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})), {{\varvec{z}}}({{\varvec{u}}}) - {{\varvec{z}}}({\widehat{{{\varvec{u}}}}})\rangle \\&\quad = -\sum _{v \in \mathcal {V}_0} \sum _{e \in \mathcal {E}(v)} (h^e(\rho ^e,w^e) - h^e(\hat{\rho }^e,{\hat{w}}^e)) (m^e(\rho ^e,w^e) - m^e({\hat{\rho }}^e,{\hat{w}}^e)) n^e \\&\qquad \quad -\sum _{v \in \mathcal {V}_\partial } \sum _{e \in \mathcal {E}(v)} (h^e(\rho ^e,w^e) - h^e({\hat{\rho }}^e,{\hat{w}}^e)) (m^e(\rho ^e,w^e) - m^e({\hat{\rho }}^e,{\hat{w}}^e)) n^e = (i)+(ii), \end{aligned}$$

into contributions coming from the junctions \(v \in \mathcal {V}_0\) and the boundary vertices \(v \in \mathcal {V}_\partial \). For the latter, we can use the results for a single pipe derived in Sect. 4, and hence

$$\begin{aligned} (ii) \le {\hat{C}}_\partial \, ( |h_\partial - {\hat{h}}_\partial |_\partial + |\varepsilon ^2 - {\hat{\varepsilon }}^2| ). \end{aligned}$$

Note that \(|h_\partial |_\partial ^2 = \sum _{v \in \mathcal {V}_\partial } |h_\partial (v)|^2\) now is the corresponding \(\ell ^2\)-norm on \(\mathbb {R}^{|\mathcal {V}_\partial |}\).

For the remaining junctions \(v \in \mathcal {V}_0\), we proceed as follows: We now make use of the coupling condition (47) and let \(h^v\) respectively \({\hat{h}}^v\) denote the uniquely determined values of the corresponding co-state variable at the junction \(v \in \mathcal {V}_0\). Then

$$\begin{aligned} \sum _{e \in \mathcal {E}(v)} (h^e(\rho ^e,w^e)&- h^e({\hat{\rho }}^e,{\hat{w}}^e)) (m^e(\rho ^e,w^e) - m^e({\hat{\rho }}^e,{\hat{w}}^e)) n^e \\&= \sum _{e \in \mathcal {E}(v)} (h^e(\rho ^e,w^e) - h^v) (m^e(\rho ^e,w^e) - m^e({\hat{\rho }}^e,{\hat{w}}^e)) n^e \\&\qquad + \sum _{e \in \mathcal {E}(v)} (h^v - {\hat{h}}^v ) (m^e(\rho ^e,w^e) - m^e({\hat{\rho }}^e,{\hat{w}}^e)) n^e \\&\qquad + \sum _{e \in \mathcal {E}(v)} ({\hat{h}}^v - h^e({\hat{\rho }}^e,{\hat{w}}^e)) (m^e(\rho ^e,w^e) - m^e({\hat{\rho }}^e,{\hat{w}}^e)) n^e \\&= (a) + (b) + (c). \end{aligned}$$

Due to the coupling condition (47), the first term (a) vanishes. Since \(h^v\) and \({\hat{h}}^v\) are both single valued on \(v \in \mathcal {V}_0\) and \(m^e({\hat{\rho }}^e,{\hat{w}}^e) = {\hat{m}}^e({\hat{\rho }}^e,{\hat{w}}^e)\), we further obtain

$$\begin{aligned} (b) = (h^v - {\hat{h}}^v) \sum _{e \in \mathcal {E}(v)} (m^e(\rho ^e,w^e) - {\hat{m}}^e({\hat{\rho }}^e,{\hat{w}}^e)) = 0, \end{aligned}$$

where we used the second coupling condition (46) for the perturbed and unperturbed problem, respectively. Using (47) for the perturbed problem, we finally get

$$\begin{aligned} (c)&= \sum _{e \in \mathcal {E}(v)} ({\hat{h}}^e({\hat{\rho }}^e,{\hat{w}}^e) - h^e({\hat{\rho }}^e,{\hat{w}}^e)) (m^e(\rho ^e,w^e) - m^e(\hat{\rho }^e,{\hat{w}}^e)) n^e \le C |\varepsilon ^2 - {\hat{\varepsilon }}^2|, \end{aligned}$$

with a constant C depending only on the bounds \({\bar{a}}, {\bar{\rho }}\) and \({\bar{w}}\) in assumption (A1)–(A2). By summation over all edges, we thus obtain the following result.

Lemma 22

Let (A1)–(A2) hold uniformly for all edges \(e \in \mathcal {E}\). Furthermore, let \({{\varvec{u}}}=(\rho ,w)\) and \({\hat{{{\varvec{u}}}}}=({\hat{\rho }},{\hat{w}})\) denote classical solutions of (42)–(48) with parameters and data \(\varepsilon ,\gamma ,h_\partial \) and \({\hat{\varepsilon }},\hat{\gamma }, {\hat{h}}_\partial \), respectively, and assume that \((\hat{\rho },{\hat{w}})\) is Lipschitz on every edge \(e \in \mathcal {E}\). Then condition (C4) holds with perturbation functional

$$\begin{aligned} P_\partial ({{\varvec{z}}}({{\varvec{u}}})-{{\varvec{z}}}({\widehat{{{\varvec{u}}}}})) = {\hat{C}}_\partial (| h_\partial - {\hat{h}}_\partial | + |\varepsilon ^2 - {\hat{\varepsilon }}^2|), \end{aligned}$$

and \({\hat{C}}_\partial \) depends only on the bounds in (A1)–(A2) and the Lipschitz bounds for \(({\hat{\rho }},{\hat{w}})\).

Remark 23

A key step in the derivation of Lemma 22 was to bound the term (ii), i.e., the contributions from interior nodes, in terms of perturbations \(|\varepsilon ^2 - {\hat{\varepsilon }}^2|\). The corresponding “relative energy flux” terms would look differently if we used the relative energy with respect to the conservative variables, as in [19]; it is not clear how to suitably bound those terms for the problem under consideration.

5.5 Stability and parabolic limit for gas networks

By the considerations of the previous sections, we immediately see that the assertions of Theorem 18 and 20 and also the remarks concerning possible generalizations remain valid verbatim also for gas networks. For the parabolic limit \(\varepsilon \rightarrow 0\), we state the corresponding result in detail.

Theorem 24

Let assumptions (A1)–(A2) hold and let \((\rho ,w)\), \(({\hat{\rho }},{\hat{w}})\) denote classical solutions of (42)–(48) with \(\varepsilon >0\), \({\hat{\varepsilon }}=0\), and \(\gamma ={\hat{\gamma }}\), \(h_\partial ={\hat{h}}_\partial \). Further assume that \(({\hat{\rho }},{\hat{w}})\) is Lipschitz on every edge \(e \in \mathcal {E}\). Then

$$\begin{aligned} \Vert \rho (\tau ) - {\hat{\rho }}(\tau )\Vert _{L^2(\mathcal {E})}^2 + \int _0^\tau \Vert w(s) - {\hat{w}}(s)\Vert _{L^3(\mathcal {E})}^3 ds&\le {\hat{C}} e^{{\hat{c}} \tau } \varepsilon ^2 , \end{aligned}$$

with constants \({\hat{c}},{\hat{C}}\) of the same form as in Theorems 18 and 20.

As similar generalization can also be made for the stability estimate of Theorem 18, and also the remarks concerning possible generalizations made for a single pipe generalize almost verbatim to gas networks.

6 Summary

In this paper, we studied the stability of classical solutions to the barotropic Euler equations with friction describing the gas transport in pipelines and pipeline networks. Our analysis was based on the formulation of a rescaled set of equations as an abstract port-Hamiltonian system involving state and co-state variables. Under some general smoothness assumptions on possible solutions, we established quantitative perturbation bounds via relative energy estimates. As a particular application, we proved convergence of sufficiently smooth solutions to solutions of the parabolic limit problem in the asymptotic high-friction, long-time, low-Mach limit \(\varepsilon \rightarrow 0\). A key ingredient of our analysis was the proper treatment of boundary conditions, which allowed us to generalize our results almost verbatim to gas-networks. Natural next steps for future investigation are the structure-preserving discretization and discretization error analysis, which can most probably be done with similar arguments. Another topic of interest might be the further investigation of the parabolic limit problem, which seems to be widely used in the gas network community.