Entropy-Preserving Coupling of Hierarchical Gas Models

TL;DR

This paper develops entropy-preserving coupling conditions for hierarchical gas flow models at pipeline junctions, ensuring well-posedness and solution uniqueness for complex transient flow scenarios.

Contribution

It introduces a novel entropy-preserving coupling framework for hierarchical gas models, including compressors, with proven existence and uniqueness of solutions.

Findings

Existence and uniqueness of solutions at junctions.

Well-posedness of Cauchy problems with small total variation.

Framework applicable to models with compressors and different fidelities.

Abstract

This paper is concerned with coupling conditions at junctions for transport models which differ in their fidelity to describe transient flow in gas pipelines. It also includes the integration of compressors between two pipes with possibly different models. A hierarchy of three one-dimensional gas transport models is built through the 3x3 polytropic Euler equations, the 2x2 isentropic Euler equations and a simplified version of it for small velocities. To ensure entropy preservation, we make use of the novel entropy-preserving coupling conditions recently proposed by Lang and Mindt [Netw. Heterog. Media, 13:177-190, 2018] and require the continuity of the total enthalpy at the junction and that the specific entropy for pipes with outgoing flow equals the convex combination of all entropies that belong to pipes with incoming flow. We prove the existence and uniqueness of solutions to generalised Riemann problems at a junction in the neighbourhood of constant coupling functions and stationary states which belong to the subsonic region. This provides the basis for the well-posedness of certain Cauchy problems for initial data with sufficiently small total variation.

Equations

$${\partial }_{t}U+{\partial }_x{F}_{{\mathcal{M}}_1}(U)$$

$$U(x,0)$$

$$U=\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right)\text{mbox}and{F}_{{\mathcal{M}}_1}(U)=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right).$$

$$p=(\gamma -1)\left(E-\frac{1}{2}\rho u^2\right).$$

$$s=c_v\ln \left(\frac{p}{{\rho }^{\gamma }}\right)+s_0,h=\frac{E+p}{\rho },\text{and}c=\sqrt{{\frac{\partial p}{\partial \rho }|}_s}=\sqrt{\frac{\gamma p}{\rho }},$$

$${\partial }_{t}U+{\partial }_x{F}_{{\mathcal{M}}_2}(U)$$

$$U(x,0)$$

$$U=\left(\begin{array}{c}\rho \\ \rho u\end{array}\right)\text{mbox}and{F}_{{\mathcal{M}}_2}(U)=\left(\begin{array}{c}\rho u\\ \rho u^2+p(\rho )\end{array}\right).$$

$$p(\rho )=\kappa {\rho }^{\gamma }\text{with}\kappa =e^{\frac{\bar{s}-s_0}{c_v}},$$

$$E=\frac{\kappa {\rho }^{\gamma }}{\gamma -1}+\frac{1}{2}\rho u^2,h=\frac{\kappa \gamma }{\gamma -1}{\rho }^{\gamma -1}+\frac{1}{2}{u}^2\text{and}c=\sqrt{\kappa \gamma {\rho }^{\gamma -1}}.$$

$$U=\left(\begin{array}{c}\rho \\ \rho u\end{array}\right)\text{mbox}and{F}_{{\mathcal{M}}_3}(U)=\left(\begin{array}{c}\rho u\\ p(\rho )\end{array}\right).$$

$$E=\frac{\kappa {\rho }^{\gamma }}{\gamma -1}\text{and}h=\frac{\kappa \gamma }{\gamma -1}{\rho }^{\gamma -1}.$$

$${\lambda }_1^{{\mathcal{M}}_1}(U)=u-c,{\lambda }_2^{{\mathcal{M}}_1}(U)=u,{\lambda }_3^{{\mathcal{M}}_1}(U)=u+c,$$

$${\lambda }_1^{{\mathcal{M}}_2}(U)=u-c,{\lambda }_2^{{\mathcal{M}}_2}(U)=u+c,$$

$${\lambda }_1^{{\mathcal{M}}_3}(U)=-c,{\lambda }_2^{{\mathcal{M}}_3}(U)=c.$$

$${\partial }_t{U}^{(j)}+{\partial }_x{F}_{m_j}(U^{(j)})$$

$$U^{(j)}(x,0)$$

$$D_{+}^{{\mathcal{M}}_1}$$

$$D_{-}^{{\mathcal{M}}_1}$$

$$D_{+}^{{\mathcal{M}}_i}$$

$$D_{-}^{{\mathcal{M}}_i}$$

$${\mathbb{I}}_i^{{\mathcal{M}}_k}:=\{j:U^{(j)}\in D_{-}^{{\mathcal{M}}_k}\}\text{and}{\mathbb{I}}_o^{{\mathcal{M}}_k}:=\{j:U^{(j)}\in D_{+}^{{\mathcal{M}}_k}\},k=1,2,3.$$

$${\mathbb{I}}_i\cup {\mathbb{I}}_o=\{1,\dots ,N\}.$$

$$\Phi \left(U^{(1)}(0^{+},t),\dots ,U^{(N)}(0^{+},t)\right)=\Pi (t),$$

$$\sum _{j=1}^N\Vert {\nu }_j\Vert q_j(0^{+},t)$$

$$h_j(0^{+},t)$$

$$s_j(0^{+},t)$$

$$s^{\ast }(t)=\frac{1}{{\sum }_{j\in {\mathbb{I}}_i}\Vert {\nu }_j\Vert q_j(0^{+},t)}\sum _{j\in {\mathbb{I}}_i}\Vert {\nu }_j\Vert (q_j{s}_j)(0^{+},t).$$

$$0=\Phi \left(U^{(1)}(0^{+},t),\dots ,U^{(N)}(0^{+},t)\right)=\left(\begin{array}{c}{\sum }_{j=1}^N\Vert {\nu }_j\Vert q_j(0^{+},t)\\ h_{j_0}(0^{+},t)-h_1(0^{+},t)\\ \cdots \\ h_{j_0}(0^{+},t)-h_{j_0-1}(0^{+},t)\\ h_{j_0}(0^{+},t)-h_{j_0+1}(0^{+},t)\\ \cdots \\ h_{j_0}(0^{+},t)-h_N(0^{+},t)\\ s_{j_1}(0^{+},t)-s^{\ast }(t)\\ \cdots \\ s_{j_{N_o}}(0^{+},t)-s^{\ast }(t)\end{array}\right)$$

$$q_1(0^{+},t)+q_2(0^{+},t)$$

$$\frac{\gamma }{\gamma -1}R{T}_1(0^{+},t)\left({\left(\frac{p_2(0^{+},t)}{p_1(0^{+},t)}\right)}^{\frac{\gamma -1}{\gamma }}-1\right)$$

$$s_1(0^{+},t)-s_2(0^{+},t)$$

$$\frac{T_2(0^{+},t)}{T_1(0^{+},t)}={\left(\frac{p_2(0^{+},t)}{p_1(0^{+},t)}\right)}^{\frac{\gamma -1}{\gamma }}.$$

$$\frac{\gamma }{\gamma -1}{C}_{p}R{q}_2(0^{+},t)T_1(0^{+},t)\left({\left(\frac{p_2(0^{+},t)}{p_1(0^{+},t)}\right)}^{\frac{\gamma -1}{\gamma }}-1\right)$$

$${\Phi }_1\left(U^{(1)}(0^{+},t),U^{(2)}(0^{+},t)\right)=\left(\begin{array}{c}{q}_1(0^{+},t)+q_2(0^{+},t)\\ \frac{\gamma R}{\gamma -1}{T}_1(0^{+},t)\left({\left(\frac{p_2(0^{+},t)}{p_1(0^{+},t)}\right)}^{\frac{\gamma -1}{\gamma }}-1\right)\\ s_1(0^{+},t)-s_2(0^{+},t)\end{array}\right),$$

$${\Phi }_2\left(U^{(1)}(0^{+},t),U^{(2)}(0^{+},t)\right)=\left(\begin{array}{c}{q}_1(0^{+},t)+q_2(0^{+},t)\\ \frac{\gamma C_{p}R}{\gamma -1}{q}_2(0^{+},t)T_1(0^{+},t)\left({\left(\frac{p_2(0^{+},t)}{p_1(0^{+},t)}\right)}^{\frac{\gamma -1}{\gamma }}-1\right)\\ s_1(0^{+},t)-s_2(0^{+},t)\end{array}\right)$$

Full text

Keywords: Conservation laws, networks, Euler equations at junctions, model hierarchy, coupling conditions of compressible fluids, compressor coupling

2010 Mathematics Subject Classification: 35L60, 35L65, 35Q31, 35R02, 76N10

1 Introduction

The transient flow of natural gas through pipeline networks in a dynamic supply-demand environment has been attracting increasing interest. Such distribution networks play an important role in future energy systems. They also allow the storage of renewable electric energy within a power-to-gas process chain. Simulation and optimisation of gas pipeline networks require the study of large scale models ranging from complex compressor stations to networks of a whole country. There exists a bunch of models based on the compressible Euler equations to predict the network behaviour with varying accuracy, see e.g. [2, 21, 22] and the nice overview in [4].

Since more accurate models are computationally more expensive, an appropriate use of a hierarchy of models is desirable. In a sequence of papers [13, 14, 15], we have developed adaptive strategies to automatically control the model selection, mainly depending on the dynamics of the gas flow. Generally, simplified models can be applied in regions with low activity, while sophisticated models have to be used in regions, where the dynamical behaviour has to be resolved in more detail.

A crucial point in the one-dimensional modelling process of gas networks is the determination of physically sound coupling conditions at junctions. Beside the natural mass and energy conservation, the equality of the dynamic pressure [12] or the pressure itself [1, 17] are widely used in the literature. The latter one is the usual choice in the engineering community. For isothermal and isentropic flows, investigations in [23, 24] showed that both pressure-based coupling conditions can deliver non-physical solutions characterized by the production of mechanical energy at a junction and should be replaced by the equality of enthalpy. Recently, we have extended this result to $3\times 3$ Euler systems with source terms at a junction of pipes with possibly different cross-sectional areas [18]. We additionally propose entropy-preserving coupling conditions, i.e., we require that the specific entropy for pipes with outgoing flow equals the convex combination of all entropies that belong to pipes with incoming flow.

In this paper, we generalize the design of entropy-preserving coupling conditions in order to account for varying models at a single junction. A hierarchy of three one-dimensional gas transport models which differ in their fidelity is built through the $3\!\times\!3$ polytropic Euler equations, the $2\!\times\!2$ isentropic Euler equations and a simplified version of it, where the kinetic energy is neglected. We also consider the practically important case of a compressor connected by two pipes with possibly different gas transport models. We first define solutions of generalised one-sided Riemann problems at a junction and show then by suitable application of the Implicit Function Theorem the existence and uniqueness of such solutions in the neighbourhood of constant coupling functions and stationary states which belong to the subsonic region. This provides the basis for the well-posedness of certain Cauchy problems for initial data with sufficiently small total variation.

The paper is organised as follows. In Sect. 2, we introduce a model hierarchy for polytropic Euler equations. Thermodynamically consistent coupling conditions are described in Sect. 3, including coupling at junctions and two models of compressor coupling. In Sect. 4, we study the solvability of generalized Riemann problems. The corresponding Cauchy problems and their solutions are studied in Sect. 5.

2 Model Hierarchy for Polytropic Euler Equations

We consider the one-dimensional polytropic Euler equations with source terms as our most accurate model $\mathcal{M}_{1}$ to describe the gas flow in a pipe of infinite length,

[TABLE]

with thermodynamic variables and flux functions

[TABLE]

Here, $\rho$ is the density, $u$ is the velocity, $p$ is the pressure, and $E$ is the total energy. Let $c_{v}>0$ and $c_{p}>0$ be the specific heat at constant volume and pressure, respectively. Then, $R=c_{p}-c_{v}$ is the gas constant and $\gamma=c_{p}/c_{v}>1$ is the adiabatic exponent. The relation between the specific internal energy $e$ and the temperature $T$ of a polytropic gas is described by $e=c_{v}T$. Together with the total energy $E=\rho e+\rho u^{2}/2$ and the ideal gas law $p=\rho RT$, the equation of state for an ideal polytropic gas in the common form reads

[TABLE]

For later use, we introduce the mass flux $q=\rho u$, the specific entropy $s$, the total enthalpy $h$, and the speed of sound $c$ defined by

[TABLE]

where $s_{0}\geq 0$ is a constant entropy value. More details about the underlying thermodynamic principles can be found, e.g., in [19, Sect. 14.4]. The right-hand side vector $G_{\mathcal{M}_{1}}(x,t,U)$ describes source terms, e.g., gravity and friction.

A first simplification for small disturbances around some background state is the use of an isentropic flow, where the entropy $s$ is taken as constant throughout the gas. In this case, we can drop the third equation in (1), i.e., the conservation of energy. The isentropic Euler equations taken as model $\mathcal{M}_{2}$ are

[TABLE]

with thermodynamic variables and flux functions

[TABLE]

Taking $s=\bar{s}$ in (5), we get an explicit relation between pressure and density,

[TABLE]

which serves now as equation of state for the isentropic Euler equations. Total energy, total entropy and speed of sound simplify to functions of $\rho$ and $u$,

[TABLE]

The isentropic equations are still nonlinear and shocks can appear, if we allow arbitrary data. Then entropy and energy will jump to a higher level across the shock, indicating the correct vanishing-viscosity solution. Although conservation of energy is no longer satisfied, such isentropic shocks may be a good approximation to reality, if they are weak enough. Further arguments are given in [19, Sect. 14.5].

In many practical situations, the spatial derivative of the kinetic energy $\partial_{x}(\rho u^{2})$ can be neglected [21]. This yields model $\mathcal{M}_{3}$ – a further simplification of (6) with thermodynamic variables and flux functions

[TABLE]

We formally set $G_{\mathcal{M}_{3}}=G_{\mathcal{M}_{2}}$ and use for the total energy and enthalpy the following approximations:

[TABLE]

The speed of sound in (10) remains unchanged.

The models $\mathcal{M}_{i}$, $i=1,2,3,$ define a hierarchy of models with decreasing fidelity. Their characteristic eigenvalues are given by

[TABLE]

In what follows, we will work within the subsonic region, i.e., $|u|<c$. In this case, only $\lambda_{2}^{\mathcal{M}_{1}}(U)$ can change its sign, depending on the velocity. It will be always clear to which model the state $U$ belongs.

3 Thermodynamically Consistent Model Coupling

3.1 Coupling at Junctions

In this section, we consider one-dimensional gas flow on a network consisting of a single junction connecting $N$ pipe sections of infinite length

[TABLE]

for $j=1,\ldots,N$. The possibly different models are identified by the parameters $m_{j}\in\mathcal{M}:=\{\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3}\}$. Each pipe is described by a vector, $\nu_{i}\in\mathbb{R}^{3}\setminus\{0\}$, originating from the common junction and parameterized by $x\in\mathbb{R}^{+}$, the real halfline $[0,\infty)$. The surface section of the pipe equals $\|\nu_{i}\|\!\neq\!0$. We assume $\nu_{i}\neq\nu_{j}$ for $i\neq j$.

For each model, we introduce two sets of subsonic data

[TABLE]

with $\accentset{\circ}{\mathbb{R}}^{+}=(0,\infty)$. Due to $\lambda_{2}^{\mathcal{M}_{1}}(U)=u$ and the orientation of the pipes, we can relate pipes modelled by $\mathcal{M}_{1}$ with a flow direction towards the junction with $D_{-}^{\mathcal{M}_{1}}$ (incoming flow), while $D_{+}^{\mathcal{M}_{1}}$ corresponds to pipes with flow direction away from the junction (outgoing flow). Since $\lambda_{1}^{\mathcal{M}_{i}}(U)<0$ and $\lambda_{2}^{\mathcal{M}_{i}}(U)>0$ for the isentropic models $\mathcal{M}_{i}$, $i=2,3$, a distinction between incoming and outgoing pipes is usually not necessary. However, this separation becomes crucial if these models are coupled with $\mathcal{M}_{1}$. This point is discussed in more detail in Section 4.

The corresponding index sets of the incoming and outgoing pipes are defined by

[TABLE]

For later use, we define the index sets $\mathbb{I}_{i}:=\mathbb{I}_{i}^{\mathcal{M}_{1}}\cup\mathbb{I}_{i}^{\mathcal{M}_{2}}\cup\mathbb{I}_{i}^{\mathcal{M}_{3}}$ and $\mathbb{I}_{o}:=\mathbb{I}_{o}^{\mathcal{M}_{1}}\cup\mathbb{I}_{o}^{\mathcal{M}_{2}}\cup\mathbb{I}_{o}^{\mathcal{M}_{3}}$ for incoming and outgoing pipes, and the special numbers $N_{o}^{\mathcal{M}_{k}}:=\dim\left(\mathbb{I}_{o}^{\mathcal{M}_{k}}\right)$ and $N^{\mathcal{M}_{k}}:=\dim\left(\mathbb{I}_{o}^{\mathcal{M}_{k}}\cup\mathbb{I}_{i}^{\mathcal{M}_{k}}\right)$. We will only consider cases with

[TABLE]

The coupling of the different model equations at the junction-pipe interface is prescribed by a set of coupling conditions of the form

[TABLE]

where $\Phi$ is a possibly nonlinear function of the traces $U^{(j)}(0^{+},t)=\lim_{x\rightarrow 0^{+}}U^{(j)}(x,t)$ of the unknown variables and $\Pi$ is a coupling constant, which depends only on time. We will use the entropy-preserving coupling conditions from [18] for $t>0$,

[TABLE]

with the entropy mix

[TABLE]

The function $h^{*}(t)$ in (H) is not prescribed and determined by the flow itself. It must be eliminated before the coupling function $\Phi$ is defined. We fix one of the enthalpy equations for a certain $j=j_{0}$, which will be specified later, and subtract from it all the other ones. This gives the coupling conditions

[TABLE]

with $N_{o}=\dim(\mathbb{I}_{o})$. Note that here $\Pi(t)\equiv 0$.

3.2 Compressor Coupling

Compressors in a network are typically placed between two pipes with equal surface section and have to be described by special coupling conditions. The task of a compressor is to increase the pressure which is permanently decreased through friction. We consider the resulting compression under adiabatic conditions, i.e., zero heat transfer between the gas and the surroundings, and as reversible process in which the entropy remains constant. This leads to the following coupling conditions (see also [20, Chapt. 4.4]):

[TABLE]

Here, $U^{(1)}\in D_{-}^{\mathcal{M}_{i}}$ and $U^{(2)}\in D_{+}^{\mathcal{M}_{j}}$, hence different models for the two pipes are allowed. The coupling constant $H^{*}(t)$ stands for the change in adiabatic enthalpy, necessary to raise the incoming pressure $p_{1}(0^{+},t)$ to the outgoing pressure $p_{2}(0^{+},t)$. The condition (CS) can be also expressed in the form

[TABLE]

In optimal control problems, $H^{*}(t)$ is often replaced by the theoretical compressor power, $P^{*}(t)=C_{p}\,qH^{*}(t),\;C_{p}=const.$ [16, 20], which can be also used as coupling constant. In this case, we have

[TABLE]

Note that $q_{2}(0^{+},t)>0$.

The coupling conditions defined above yield the two coupling functions

[TABLE]

[TABLE]

and the corresponding $\Pi$-functions in (24),

[TABLE]

4 Generalized Riemann Problems for Model Coupling

4.1 Coupling at Junctions

In this section we will show that the coupling conditions (M), (H) and (S) for the network system (16) – (17) with $G=0$ are well-defined. Following the theoretical framework applied in [5, 12, 18], we consider a generalised Riemann problem at a junction connecting pipes with different gas flow models, and show that there exists a unique self-similar solution in terms of the classical Lax solution to standard Riemann problems.

Let us denote by $\Omega_{j}=\{U^{(j)}\in D_{+}^{m_{j}}\cup D_{-}^{m_{j}}\}$ for $j=1,\ldots,N$ and $m_{j}\in\mathcal{M}$ nonempty sets and define the overall state space $\Omega=\Omega_{1}\times\Omega_{2}\times\ldots\times\Omega_{N}$.

Definition 4.1**.**

The generalized Riemann problem at a junction in $x=0$ with $N$ adjacent pipes with different flow models is defined through the set of equations

[TABLE]

where the states $\bar{U}_{0}^{(1)},...,\bar{U}_{0}^{(N)}$ are constant states in $\Omega$ and $\bar{\Pi}\in\mathbb{R}^{d}$ is a constant vector of dimension $d=N+N_{o}^{\mathcal{M}_{1}}$.

Definition 4.2**.**

A $\Phi$-solution to the Riemann problem (31) is a self-similar function $U(x,t):\mathbb{R}^{+}\times\mathbb{R}^{+}\rightarrow\Omega$ for which the following hold:

There exists a constant state $U_{\ast}(\bar{U}_{0})=\lim_{x\rightarrow 0^{+}}U(x,t)$ such that all components $U^{(j)}(x,t)$ coincide with the restriction to $x>0$ of the Lax solution to the standard Riemann problem for $x\in\mathbb{R}$,

[TABLE] 2. 2.

The state $U_{\ast}$ satisfies $\Phi(U_{\ast})=\bar{\Pi}$ for all $t>0$.

Riemann solution for isentropic Euler equations. The solution of the standard Riemann problem (32) for the isentropic models $\mathcal{M}_{2}$ and $\mathcal{M}_{3}$ with initial data $(U_{L},U_{R})$ for $x<0$ and $x>0$, respectively, can be described by a set of elementary waves such as rarefaction and shock waves.

These waves are parameterisations of the Rankine-Hugoniot jump condition and the Riemann invariants [19, 25]. Due to the construction of the network and the subsonic flow conditions, only $1$-waves can hit the junction, whereas $2$-waves leave the junction. These waves separate the solution in three states $(U_{L},U_{\ast}^{m},U_{R})$, $m\in\{\mathcal{M}_{2},\mathcal{M}_{3}\}$, see Fig. 2. The components of $U_{\ast}^{m}=(\rho_{\ast}^{m},q_{\ast}^{m})$ are determined by the following equations [6]:

[TABLE]

where

[TABLE]

Riemann solution for polytropic Euler equations. For the polytropic Euler equations, the set of waves is extended by a contact discontinuity which is located between the other two, see Fig. 2.

Here, $1$-waves enter the junction, $2$- and $3$-waves leave the junction while separating the solution in four states $(U_{L},U_{L\ast},U_{R\ast},U_{R})$. The velocity and the pressure are constant across the contact discontinuity, i.e., we have

[TABLE]

The four sought variables $(p_{\ast},u_{\ast},\rho_{L\ast},\rho_{R\ast})$ are again implicitly defined by means of parameterizations [25, Chapt. 4], [19, Chapt. 14.11]. It holds

[TABLE]

where for $k=L,R$

[TABLE]

with $\mu^{2}=(\gamma-1)/(\gamma+1)$ and $c_{k}^{2}=\gamma p_{k}/\rho_{k}$. The parameter $p_{\ast}$ is determined from the second equality in (38). The functions $\psi(p_{\ast},U_{k})$ and $\phi(p_{\ast},U_{k})$ are twice continuously differentiable at $p_{\ast}=p_{k}$. The total energy for the inner region is derived from $E_{k\ast}=p_{\ast}/(\gamma-1)+\rho_{k\ast}u_{\ast}^{2}/2$ for $k=L,R$.

Lax curves. By means of the Riemann solutions, we can set up the parameterisations of the $k$-waves, the so called $k$-Lax curves. For the isentropic models, $\mathcal{M}_{2}$ and $\mathcal{M}_{3}$, they are given by

[TABLE]

The Lax-curves for the model $\mathcal{M}_{1}$ read

[TABLE]

Coupling conditions for the $\Phi$-solution. Since in our network modelling all pipes are only outgoing from a junction with respect to $x$-coordinate, the sign of the velocity in incoming pipes (w.r.t. flow direction) has to be changed when switching from the standard to the generalised Riemann problem. This changes the parameterisations of all $\mathcal{L}_{1}$-curves. Indeed, $\mathcal{L}_{1}^{m}$ has to be replaced by $\mathcal{L}_{2}^{m}$ for $m=\mathcal{M}_{2},\mathcal{M}_{3}$, and $\mathcal{L}_{1}^{\mathcal{M}_{1}}$ by $\mathcal{L}_{3}^{\mathcal{M}_{1}}$. We also note that a contact discontinuity always travels with positive wave speed. This is a consequence of constant initial data, the spatial parametrisation of the pipes and the restriction to subsonic flow.

The coupling conditions in (31) can now be expressed in terms of the Lax-curves, see Fig. 3 for a schematic illustration. Let us set $U_{\ast}=U_{L\ast}$ for $\mathcal{M}_{1}$. Then, the sought state $U_{\ast}$ in Def. 4.2 satisfies

[TABLE]

with

[TABLE]

Given constant states $\bar{U}^{(j)}$, mass flux, enthalpy and entropy can be extracted from $U^{(j)}_{\ast}$ using the Lax curves:

[TABLE]

with $f_{j}\in\{q_{j},h_{j},s_{j}\}$, where $h_{j}$ and $s_{j}$ depend on the model chosen and are given by the expressions

[TABLE]

Without loss of generality, let $\mathbb{I}_{o}^{\mathcal{M}_{1}}=\{1,2,\ldots,n_{0}\}$, $\mathbb{I}_{i}^{\mathcal{M}_{1}}=\{n_{0}+1,\ldots,n_{1}\}$, $\mathbb{I}_{i}^{\mathcal{M}_{2}}=\{n_{1}+1,\ldots,n_{2}\}$, and $\mathbb{I}_{i}^{\mathcal{M}_{3}}=\{n_{2}+1,\ldots,n_{3}\}$. Accordingly, the free parameters are $(\sigma_{1},\ldots,\sigma_{N})$ and $(\tau_{1},\ldots,\tau_{n_{0}})$. The coupling conditions (M), (H), and (S) can now be written as

[TABLE]

with $s^{\ast}$ defined through

[TABLE]

We set $s_{j}=s^{\ast}$ for $j\in\mathbb{I}_{o}^{\mathcal{M}_{2}}\cup\mathbb{I}_{o}^{\mathcal{M}_{3}}$, which is always possible, if $s^{\ast}$ is well-defined by the coupling conditions (54). The regularity of the Lax curves ensures the property $\Phi\in C^{1}(\mathbb{R}^{N}\times\mathbb{R}^{n_{0}},\mathbb{R}^{d})$. It remains to show that (54) has a unique solution. Newton’s method can then be applied to determine the solution vector $(\sigma^{\ast},\tau^{\ast})$, which determines the states $U_{\ast}^{(j)}$. For the well-posedness of the generalised Riemann problem (31) with the coupling function $\Phi$ defined in (54), the following local result as a generalisation of [18, Theorem 2.1] can be given.

Theorem 4.1**.**

Let $N>\dim(\mathbb{I}_{i})>0$ and $\Phi$ as defined in (54). Assume constant initial data $\bar{U}^{(j)}\in D_{-}^{m},\;j\in\mathbb{I}_{i}^{m}$, and $\bar{U}^{(j)}\in D_{+}^{m},\;j\in\mathbb{I}_{o}^{m}$, $m\in\{\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3}\}$, with $\Phi(\bar{U})\!=\!0$ are given. Then there exist positive constants $\delta$ and $K$ such that for all initial states $\tilde{U}\in\Omega$ with $\sum_{j=1,\ldots,N}\|\tilde{U}^{(j)}-\bar{U}^{(j)}\|\!<\!\delta$, the Riemann problem (31) admits a unique $\Phi$-solution $U(x,t)=\mathcal{R}^{\Phi}(\tilde{U})$ satisfying $\Phi(U(0^{+},t))\!=\!0$ and

[TABLE]

Additionally, if $\nu$ is replaced by $\hat{\nu}$, where $\sum_{i=1,\ldots,N}\|\nu_{i}-\hat{\nu}_{i}\|\!<\!\delta$, and $\mathcal{R}^{\Phi}_{\hat{\nu}}(\tilde{U})$ is the corresponding $\Phi$-solution for the same initial state $\tilde{U}$, then

[TABLE]

with $\mathcal{R}^{\Phi}_{\nu}(\tilde{U})\!=\!\mathcal{R}^{\Phi}(\tilde{U})$.

Proof.

In the spirit of the implicit function theorem, it is sufficient to study the determinant of the Jacobian matrix $D_{(\sigma,\tau)}\Phi(\sigma_{0},\tau_{0})$ with the two argument vectors $\sigma_{0}=(\bar{p}_{1},\ldots,\bar{p}_{n_{1}},\bar{\rho}_{n_{1}+1},\ldots,\bar{\rho}_{n_{3}})$ and $\tau_{0}=0\in\mathbb{R}^{n_{0}}$. Note that $\Phi(\sigma_{0},\tau_{0})=0$. The Jacobian reads

[TABLE]

Here, we have used the notations $f_{\mu_{i}}=\partial_{\mu_{i}}f_{i}$, $\hat{q}_{\mu_{i}}=\|\nu_{i}\|\partial_{\mu_{i}}q_{i}$ for $f\!=\!h,s$ and $\mu=\sigma,\tau$, and $s^{\ast}_{\sigma_{i}}=\partial_{\sigma_{i}}s^{\ast}$. Note that $s^{\ast}_{\sigma_{j}}=0$ for $j>n_{3}$. From (51), we derive the following derivatives:

[TABLE]

where $\bar{c}_{j}=\sqrt{\gamma\bar{p}_{j}/\bar{\rho}_{j}}$ for $\mathcal{M}_{1}$ and $\bar{c}_{j}=\sqrt{\kappa_{j}\gamma\bar{\rho}_{j}^{\gamma-1}}$ for $\mathcal{M}_{2},\mathcal{M}_{3}$. None of the derivatives can vanish, except $\partial_{\sigma_{j}}s_{j}$ and $\partial_{\sigma_{j}}s^{\ast}$. Without loss of generality, we number the incoming pipes in such a way that $\bar{s}_{n_{0}+1}=\max_{j\in\mathbb{I}_{i}}\bar{s}_{j}$. Then $\bar{s}_{n_{0}+1}-\bar{s}^{\ast}\geq 0$, and therefore $s^{\ast}_{\sigma_{n_{0}+1}}\leq 0$ since $\bar{q}_{j}<0$ for $j\in\mathbb{I}_{i}$. Note that $\dim(\mathbb{I}_{i})>0$ has been assumed.

Case 1: $n_{0}>0$.

A closer inspection of the special structure of the Jacobian (58) reveals that it is regular if and only if all $3\times 3$–matrices

[TABLE]

are regular. Taking into account the sign of the derivatives, we find

[TABLE]

Case 2: $n_{0}=0$.

In this case, model $\mathcal{M}_{1}$ does not appear for outgoing pipes. Hence, the entropy mix $s^{\ast}$ is simply passed as constant entropy value to the lower order models applied in the outflow region. The parameters $\tau_{0}$ disappear and the Jacobian matrix has the simplified form

[TABLE]

We first note that $\hat{q}_{\sigma_{j}},h_{\sigma_{j}}\!>\!0$ for all $j=1,\ldots,N$. Consequently, the column vectors are linearly dependent if and only if the first vector can be written as the sum of the others. It would request that $\hat{q}_{\sigma_{1}}=\alpha_{2}\hat{q}_{\sigma_{2}}+\ldots+\alpha_{N}\hat{q}_{\sigma_{N}}$ with $\alpha_{j}=-h_{\sigma_{1}}/h_{\sigma_{j}}$, which contradicts $\hat{q}_{\sigma_{1}}>0$. This shows $\det(D_{\sigma}\Phi(\sigma_{0}))\neq 0$.

Now, the implicit function theorem ensures the existence of a $\delta\!>\!0$, a neighbourhood $\mathcal{U}(v_{0})$ of $v_{0}=(\sigma_{0},\tau_{0})$ in case 1 or $v_{0}=\sigma_{0}$ in case 2, and a function $\varphi:B(\bar{U},\delta)\rightarrow\mathcal{U}(v_{0})$ such that $\varphi(\bar{U})=v_{0}$ and $\Phi(v;U)=0$ if and only if $v=\varphi(U)$ for all $U\in B(\bar{U},\delta)$. The solution $U(x,t)$ can then be identified by the restriction to $x\in\mathbb{R}^{+}$ of the solution to the standard Riemann problem (32) with $\bar{U}_{0}\!=\!\tilde{U}$. The Lipschitz estimate (56) follows from the $C^{1}$-regularity of $\Phi$. Since $\Phi$ depends smoothly on $\|\nu_{i}\|$, the same arguments as above can be used to show (57). ∎

4.2 Compressor Coupling

We will now show that the compressor couplings, i.e., the coupling conditions (CM) and (CS) accomplished with either (CP1) or (CP2), are also well-defined. The generalized Riemann problem for the compressor and its self-similar $\Phi$-solution can be formulated using Definition 4.2 with $N\!=\!2$ and $\Phi$ given by the corresponding compressor model.

We have the following theorem:

Theorem 4.2**.**

Let $N\!=\!2$, $\|\nu_{1}\|\!=\!\|\nu_{2}\|$ and $\Phi$ defined through (CM), (CS), and either (CP1) or (CP2). Assume constant initial data $\bar{U}^{(1)}\in D_{-}^{i}$ and $\bar{U}^{(2)}\in D_{+}^{j}$ with $i,j\in\{\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3}\}$ and $\Phi(\bar{U})\!=\!\bar{\Pi}\!=\!(0,\bar{\Pi}_{2},0)^{T}$ with $\bar{\Pi}_{2}>0$ are given. Then there exist positive constants $\delta$ and $K$ such that for all states $\tilde{U}\in\Omega$ with $\|\tilde{U}^{(1)}-\bar{U}^{(1)}\|+\|\tilde{U}^{(2)}-\bar{U}^{(2)}\|\!<\!\delta$, the Riemann problem (31) admits a unique $\Phi$-solution $U(x,t)=\mathcal{R}^{\Phi}(\tilde{U})$ satisfying $\Phi(U(0^{+},t))\!=\!\bar{\Pi}$ and

[TABLE]

Proof.

We proceed as in Theorem 4.1 and study the determinants of the Jacobian matrices. Let us first start with (CP1) and the case where model $\mathcal{M}_{1}$ is taken at the outflow. Then the free parameter vector is $(\sigma_{1},\sigma_{2},\tau_{2})$ and the coupling conditions can be written in the form $\Phi(\sigma_{1},\sigma_{2},\tau_{2})-\bar{\Pi}=0$. More precisely, we have

[TABLE]

with $C\!=\!R\gamma/(\gamma-1)\!>\!0$. The derivatives taken at constant state values read

[TABLE]

and

[TABLE]

for $k=2,3$. The Jacobian of $\hat{\Phi}$ evaluated at $(\bar{p}_{1},\bar{p}_{2},0)$ if $1\in\mathbb{I}_{i}^{\mathcal{M}_{1}}$ or $(\bar{\rho}_{1},\bar{p}_{2},0)$ if $1\in\mathbb{I}_{i}^{\mathcal{M}_{2}}\cup\mathbb{I}_{i}^{\mathcal{M}_{3}}$ has the form

[TABLE]

where we used the notation $f_{\mu_{i}}=\partial_{\mu_{i}}f_{i}$ for $f=q,s$, and $\mu=\sigma,\tau$. A short calculation of all derivatives of the second coupling condition $\hat{\Phi}_{2}$ reveals

[TABLE]

Together with $s_{\tau_{2}}<0$, $q_{\sigma_{1}}>0$, and $q_{\sigma_{2}}>0$, this shows

[TABLE]

If one of the models $\mathcal{M}_{2}$ or $\mathcal{M}_{3}$ is used in the outflow region, the entropy equality becomes trivial and the coupling conditions reduce to

[TABLE]

Due to $q_{2}^{\prime}(\bar{\rho}_{2})\!=\!\lambda_{2}^{\mathcal{M}_{k}}>0$ and $p_{2}^{\prime}(\bar{\rho}_{2})\!=\!\bar{c}_{2}^{2}>0$, we finally conclude that

[TABLE]

It remains to study the case (CP2). Again we start with model $\mathcal{M}_{1}$ at the outflow. Then, the second component of the coupling conditions reads

[TABLE]

with $C=C_{p}R\gamma/(\gamma-1)>0$. We note $q_{2}>0$ and have $\bar{p}_{2}>\bar{p}_{1}$ due to the assumption $\bar{P}^{\ast}>0$. This gives the inequalities

[TABLE]

and eventually $\det(D\hat{\Phi})>0$. Simplifying the model in the outflow region to $\mathcal{M}_{2}$ or $\mathcal{M}_{3}$ leads to the coupling conditions

[TABLE]

The same arguments as above show $\det(D\hat{\Phi})<0$.

In all cases, the implicit function theorem guarantees the existence of a unique $\Phi$-solution to the Riemann problem (31) for initial values $\tilde{U}$ varying in a small neighborhood of $\bar{U}$. ∎

5 The Cauchy Problem at the Junction

We first introduce a few notations.

Definition 5.1**.**

Let

[TABLE]

For the extended variable ${\bf q}=(U,\Pi)$, a constant state $\bar{U}$ and a constant vector $\bar{\Pi}$, we consider the metric space

[TABLE]

equipped with the distance and total variation (TV)

[TABLE]

For positive $\delta\in[0,\bar{\delta}]$, we set $D_{\delta}({\bf q})\!=\!\{{\bf q}\in X:TV({\bf q})\leq\delta\}$ and introduce the set of varying states $\mathcal{U}_{\delta}(\bar{U})\!=\!\{U\in\bar{U}+{\bf L}^{1}(\mathbb{R}^{+};\Omega):TV(U)\leq\delta\}$ in the neighborhood of $\bar{U}$.

Let $G$ denote the vector of the right-hand side functions in (16) for all pipes and be defined through

[TABLE]

For the map $G:[0,T]\times\mathcal{U}_{\bar{\delta}}(\bar{Y})\rightarrow{\bf L}^{1}(\mathbb{R}^{+};\Omega)$, we assume that there exist positive constants $L_{1}$ and $L_{2}$ such that for all $t,s\in[0,T]$ the following inequalities are satisfied:

[TABLE]

This is the usual assumption on $G$, which also covers non-local terms [7, 8] as well as real applications [10].

Next we define the Cauchy problem at junctions, which corresponds to our special set of coupling conditions, and weak solutions.

Definition 5.2**.**

Let $N>\dim(\mathbb{I}_{i})>0$ and $\Phi$ defined through (M), (H), (S), or $N\!=\!2$, $\|\nu_{1}\|\!=\!\|\nu_{2}\|$ and $\Phi$ defined through (CM), (CS), and either (CP1) or (CP2). A weak solution $U=(U^{(1)},\ldots,U^{(N)})$ on $[0,T]$ to the Cauchy problem

[TABLE]

is a function $U\in{\bf C}^{0}([0,T];\bar{U}+{\bf L}^{1}(\mathbb{R}^{+};\Omega))$ such that $U(t)\in{\bf BV}(\mathbb{R}^{+};\Omega)$ for all $t\in[0,T]$, the initial conditions, $U(x,0)\!=\!U_{0}(x)$, and the condition at the junction, $\Phi(U(0^{+},t))=\Pi(t)$, for a.e. $t>0$ are satisfied. Further, for all $\varphi\in{\bf C}^{\infty}_{c}(\mathbb{R}^{+}\times(0,T);\mathbb{R})$ it holds

[TABLE]

The weak solution is entropic if for all non-negative $\varphi\in{\bf C}^{\infty}_{c}(\accentset{\circ}{\mathbb{R}}^{+}\times(0,T);\mathbb{R}^{+})$ and $i=1,\ldots,N$ it holds

[TABLE]

for all convex entropy-entropy flux pairs $(\eta_{m_{i}}(U^{(i)}),\psi_{m_{i}}(U^{(i)}))$ of the model $m_{i}\in\mathcal{M}$.

5.1 The homogeneous Cauchy problem

We first start with the homogeneous Cauchy problem, i.e, we set $G_{m_{i}}=0$ in (86) for all $i=1,\ldots,N$. The function $\Pi(t)$ is assumed to be constant, $\Pi(t)=\bar{\Pi}$, to cover the situations discussed in Theorem 4.1 and 4.2. In this case, a solution of the Cauchy problem can be constructed applying a proper modification of the wave front tracking method by Bressan [3] and its natural extension for networks introduced by Colombo and Mauri [12]. In terms of the standard Riemann semigroup, we have the following

Theorem 5.1**.**

Let $G_{m_{i}}\!=\!0$ for all $i=1,\ldots,N$, $\Pi(t)=\bar{\Pi}$, and the assumptions of Theorem 4.1 or 4.2 be satisfied with constant values $\bar{{\bf q}}=(\bar{U},\bar{\Pi})$. Then there exist positive constants $\delta$, $K$, a domain $D$, and a semigroup $S:\mathbb{R}^{+}\times D\rightarrow D$ such that

1. (1)

$\text{cl}_{{\bf L}_{1}}\mathcal{U}_{\delta}(\bar{U})\subseteq D$. 2. (2)

For all $U\in D$, $S_{0}(U)=U$ and $S_{s}S_{t}(U)=S_{s+t}(U)$. 3. (3)

For all $U\in D$, the map $t\rightarrow S_{t}(U)$ is a weak entropic solution to the homogeneous Cauchy problem (86) in the sense of Definition 5.2. 4. (4)

For $\hat{U},\tilde{U}\in D$ and $s,t\geq 0$ it holds

[TABLE] 5. (5)

If $U\in D$ is piecewise constant and $t>0$ sufficiently small, then $S_{t}(U)$ coincides with the juxtaposition of the solutions to Riemann problems centered at the points of jumps or at the junction.

Proof.

The proof is based on the wave front tracking algorithm, see [3]. The initial data $U_{0}$ is approximated by a piecewise constant function $\bar{U}_{0}$ with a finite number of discontinuities such that $TV(\bar{U}_{0})\leq\delta$ and $\|U_{0}-\bar{U}_{0}\|_{{\bf L}^{1}}\leq\varepsilon$. Next, we approximate the local Riemann problems at each point of a jump in the inner of a pipe by an approximate Riemann solver - the accurate solver from [3, Sect. 7.2]. Rarefaction waves are substituted by rarefaction fans as specified by the accurate solver. At junctions, we solve the Riemann problems by means of our solution procedures introduced above. If $\delta$ is sufficiently small, Theorems 4.1 and 4.2 ensure the existence and uniqueness of corresponding solutions. This construction can be continued up to the first collision among waves in a pipe or until a wave hits the junction. We then apply the simplified solver from [3, Sect. 7.2] and its natural extension to junctions from [12, Sect. 4.2] in order to keep the total number of waves finite. At a junction, waves with small strength are reflected into a non-physical wave so that no wave in any other pipe is produced. For $\lim_{k\rightarrow\infty}\varepsilon_{k}=0$, the algorithm produces a sequence of $\varepsilon$-solutions $\{U^{\varepsilon_{k}}\}_{k\in\mathbb{N}}$ in the sense of Bressan [3, Lemma 7.1]. We will now show that these solutions fulfill Helly’s embedding theorem [3, Theorem 2.3].

Fix a wave front tracking approximate solution $U^{\varepsilon}$. Then we have to show the following three properties:

[TABLE]

for some $C,M,L>0$. Let $\mathcal{J}_{i}$ and $\mathcal{A}_{i}$ denote the set of discontinuities and the set of approaching wave fronts in the $i$-th pipe for every $t>0$, respectively, see [3, Sect. 7.3]. The strengths of the waves of the first, second and third family are denoted by $v_{j,\alpha}$, $j=1,2,3$ and $\alpha\in\mathcal{J}_{i}\cup\mathcal{A}_{i}$ with $v=(\sigma,\tau)$ for outgoing waves of the most accurate model $\mathcal{M}_{1}$ and $v=\sigma$ otherwise. The notation $v_{\alpha}$ is used when the wave type is not explicitly specified. We now define Glimm-type functionals

[TABLE]

with a suitable constant $\hat{K}_{J}>0$, which will be specified later, and the functionals

[TABLE]

The interaction of waves for all models considered above are well understood. Estimates for the interaction functional $Q$ are given, e.g., in [8, Lemma 4.1] and [12, Sect. 4.2.]. Here, we focus on the interaction at junctions. In what follows, all wave fronts before interaction time $t$, at which a shock or rarefaction fan hits the junction, are denoted by $v^{-}$ and wave fronts after interaction time $t$ resulting from our Riemann solver are denoted by $v^{+}$. We first recall a simple consequence of the implicit function theorem applied in the proofs of Theorems 4.1 and 4.2. With the function $\varphi:B(\bar{U},\delta)\rightarrow\mathcal{U}(v_{0})$, where $\varphi(U)=v$ and $\varphi(\bar{U})=v_{0}$, there exists a constant $C>0$ such that for all $\tilde{U}\in B(\bar{U},\delta)$ with $\varphi(\tilde{U})=\tilde{v}$ it holds

[TABLE]

Identifying the states before and after interaction by $U^{-}$ and $U^{+}$ with parameters $v^{-}$ and $v^{+}$ and setting $v_{0}=0$, which is always possible by a change of the coordinate system, we get

[TABLE]

with a positive constant $\tilde{K}$. From this, we conclude with the general estimate

[TABLE]

We are now ready to choose appropriate constants $C_{i}(v_{\alpha})$ in (91). We set

[TABLE]

We will now show the boundedness of the total variation of $U^{\varepsilon}(x,t)$. Changes of the functional $V$ as function of $t$ can be expressed in the form

[TABLE]

with

[TABLE]

Thus, due to estimate (94), we have $\triangle V_{k}(t)\leq-K_{J}V_{k}(t-)$, $k=1,2,3$, and therefore also $\triangle V(t)\leq-K_{J}V(t-)$. Consequently, the map $t\rightarrow V(t)$ is non-increasing over time at the junction. For the map $t\rightarrow Q(t)$, we proceed analogously to [6, 12] and deduce

[TABLE]

from which we get the estimate

[TABLE]

Let now $\hat{K}_{J}$ be such that $\hat{K}_{J}V(0)-K_{J}<0$. Then $\triangle Y(t)\leq 0$, showing that the map $t\rightarrow Y(t)$ is also non-increasing. Observe that the TV-norm is equivalent to $V$, i.e., there exists a constant $C_{1}>0$ such that $TV(U)\leq C_{1}\,V(U)$ and $V(U)\leq C_{1}\,TV(U)$ for all $U\in\mathcal{U}_{\delta}(\bar{U})$. Thus, we can estimate

[TABLE]

where we have used that $TV(U^{\varepsilon}(\cdot,0))\leq\delta$. Boundedness of $U^{\varepsilon}$ follows from the fact that all $\varepsilon$-solutions are piecewise constant with a finite number of discontinuities and have bounded variation. The stability with respect to time - the third condition in (89) - has been proven in [3, Sect. 7] independently of the underlying geometry. As a consequence, an application of Helly’s compactness theorem yields the convergence of a subsequence $\{U^{\varepsilon_{l}}\}_{l\in\mathbb{N}}$, where the limit function preserves the properties of the $\varepsilon$-functions. The limit orbits $t\rightarrow S_{t}(U)$ are solutions in the sense of Definition 5.2, see [3, Lemma 7.1]. This proves the existence of a semigroup $S:\mathbb{R}^{+}\times D\rightarrow D$ and a domain $D$ such that $\text{cl}_{{\bf L}^{1}}\mathcal{U}_{\delta}(\bar{U})\subseteq D$. Hence, statements (1), (2), (3), and (5) clearly hold.

We pass now to the ${\bf L}^{1}$-stability of the $\varepsilon$-solutions with respect to the initial data, i.e., statement (4). For this aim, we will use stability functionals from [11, 12]. Let $U^{\varepsilon}$, $V^{\varepsilon}\in\mathcal{U}_{\delta}(\bar{U})$ be two piecewise constant $\varepsilon$-solutions and denote by $q_{i}^{(l)}(x)$ the shock sizes implicitly given as solutions of

[TABLE]

where $\mathbb{I}^{\mathcal{M}_{j}}$ is the index set of all pipes with model $\mathcal{M}_{j}$ and $\mathcal{L}_{l}$, $l=1,2,3$, are the Lax-curves for the shock waves (see also [3, (8.4)]). For a compact notation, we set $q_{i}^{(3)}(x)=0$ for $i\in\mathbb{I}^{\mathcal{M}_{2}}\cup\mathbb{I}^{\mathcal{M}_{3}}$. Now define the functional

[TABLE]

with weights

[TABLE]

The constant $K_{s}$ is defined by

[TABLE]

with $W_{max}=\max\{\max_{i,m_{i}=1}W_{i}^{(3)},\max_{i,m_{i}=2,3}W_{i}^{(2)}\}$ and $K_{J}$ from (94). Further, $Y$ is the functional defined in (90), while the functions $A_{i}^{(j)}$ are defined in dependence of the models $\mathcal{M}_{m_{i}}$ on the $i$-th pipe by $A_{i}^{(j)}:=A_{i,j}^{\mathcal{M}_{m_{i}}}$ with

[TABLE]

and

[TABLE]

These definitions are standard and can be found, e.g., in [3, (8.9)], [12, p. 566], and [11]. Now we fix $\kappa_{1}$ and $\kappa_{2}$ in such a way that $\delta>0$ in the definition of $\mathcal{U}_{\delta}$ can be chosen to uniformly satisfy $1\leq W_{i}^{(j)}(x)\leq\kappa$ for every $i=1,\ldots,N$ and $j=1,2,3$. Indeed, since $K_{s}\geq 1$ and $A_{i}^{(j)}\geq 0$, the lower bound is obvious. By definition, we have $A_{i}^{(j)}\leq\kappa_{0}(Y(U^{\varepsilon}(\cdot,t))+Y(V^{\varepsilon}(\cdot,t)))$ with a certain constant $\kappa_{0}>0$ and consequently $W_{i}^{(j)}\leq K(1+2\delta(\kappa_{0}\kappa_{1}+\kappa_{1}\kappa_{2}))$, where $K=1$ for $j=2,3$ and $K=K_{s}$ for $j=1$. Hence, the functional $\Theta$ is equivalent to the ${\bf L}^{1}$ distance, i.e.,

[TABLE]

with a positive constant $C_{1}$. Applying the same calculations as in the proof of [3, Theorem 8.2] shows that at any time $t>0$ when an interaction happens neither in $U^{\varepsilon}$ nor $V^{\varepsilon}$, we have

[TABLE]

with a positive constant $C_{2}$ that depends only on a bound on the total variation of the initial data. At an interaction time $t>0$, we get from above $\triangle(Y(U^{\varepsilon}(\cdot,t))+Y(V^{\varepsilon}(\cdot,t)))<0$ and hence, by choosing $\kappa_{2}$ large enough, we obtain

[TABLE]

Thus, $\Theta(U^{\varepsilon}(\cdot,t),V^{\varepsilon}(\cdot,t))-\Theta(U^{\varepsilon}(\cdot,s),V^{\varepsilon}(\cdot,s))\leq C_{2}(t-s)$ for all $0\leq s\leq t$. This concludes the proof of statement (4) by standard arguments given in [3, Sect. 8.3]. ∎

5.2 The inhomogeneous Cauchy problem

Including the source terms $G_{m_{i}}$ in (86) for all $i=1,\ldots,N$, we have the following result for the well-posedness of the Cauchy problem:

Theorem 5.2**.**

Let the assumptions of Theorem 4.1 or 4.2 be satisfied with constant values $\bar{{\bf q}}=(\bar{U},\bar{\Pi})$. Let $\mathcal{T}_{t}$ be the right translation defined by $(\mathcal{T}_{t}\Pi)(s)=\Pi(t+s)$. Then there exist positive constants $\delta$, $\delta^{\prime}$, $K$, domains $D^{t}$ for $t\in[0,T]$, and a map $\mathcal{E}(s,t_{0}):D^{t_{0}}\rightarrow D_{\delta}$ with $t_{0}\in[0,T]$ and $s\in[0,T-t_{0}]$ such that

1. (1)

$D_{\delta^{\prime}}(\bar{{\bf q}})\subseteq D^{t}\subseteq D_{\delta}(\bar{{\bf q}})$* for all $t\in[0,T]$.* 2. (2)

$\mathcal{E}(0,t_{0}){\bf q}={\bf q}$* for all $t_{0}\in[0,T],\;{\bf q}\in D^{t}$.* 3. (3)

$\mathcal{E}(s,t_{0})D^{t_{0}}\subset D^{t_{0}+s}$* for all $t_{0}\in[0,T],\;s\in[0,T-t_{0}]$.* 4. (4)

For all $t_{0}\in[0,T]$, $s_{1},s_{2}\geq 0$ with $s_{1}+s_{2}\in[0,T-t_{0}]$

[TABLE] 5. (5)

For all $(U_{0},\Pi)\in D^{t_{0}}$, the map $t\rightarrow\mathcal{E}(t,t_{0})(U_{0},\Pi)=(U(t),\mathcal{T}_{t}\Pi)$ is the entropic solution to the Cauchy problem (86) in the sense of Definition 5.2. 6. (6)

For all $t_{0}\in[0,T]$ and $(U_{0},\Pi)\in D^{t_{0}}$

[TABLE]

where $(U(t),\mathcal{T}_{t}\Pi)=\mathcal{E}(t,t_{0})(U_{0},\Pi)$ and $S_{t}$ denotes the semigroup generated from (86) with $G=0$. 7. (7)

For all $t_{0}\in[0,T],\;s\in[0,T-t_{0}]$ and ${\bf q},\tilde{{\bf q}}\in D^{t_{0}}$

[TABLE]

Proof.

The proof of this theorem is based on the operator splitting technique introduced by Colombo and Guerra [9] in the framework of differential equations in metric spaces and the techniques applied by Colombo, Guerra, Herty, and Schleper [10]. In the metric space $X$ defined in (82), the following map is considered:

[TABLE]

It approximates the solution of the inhomogeneous Cauchy problem in the time interval $[t_{o},t_{o}+t]$. Let ${\bf q}=(U,\Pi)$, $\tilde{{\bf q}}=(\tilde{U},\tilde{\Pi})$ and extent the stability function to

[TABLE]

This extension keeps all properties of the original function $\Theta$ and therefore there exists again a positive constant $C_{1}$ such that for all ${\bf q}$, $\tilde{{\bf q}}\in D_{\delta}$ it holds

[TABLE]

with $d_{X}$ defined in (83). Then, the same calculation as in [10, Lemma 4.7] yields

[TABLE]

for all $t_{o}\in[0,T]$, for all ${\bf q},\tilde{{\bf q}}\in D$ (the domain of the semigroup $S$ from Theorem 5.1), $\varepsilon>0$ sufficiently small and a positive constant $C$. This also gives the Lipschitz dependence of $F$ from $t$ and ${\bf q}$, see [10, Proposition 4.8]. Following [9, Definition 2.2], we next construct the Euler $\varepsilon$-polygonal $F^{\varepsilon}$ generated by $F$ on the interval $[t_{o},t_{o}+t]$,

[TABLE]

where $k=\max\{h\in\mathbb{N}:h\varepsilon<t_{o}+t\}$. In such a way, the solution of the inhomogeneous Cauchy problem is approximated by the solution of the homogeneous Cauchy problem, corrected by an Euler approximation of the source term after every time step of length $\varepsilon$. It follows from [10, Proposition 4.9] that this Euler approximation is a local flow in the sense of [9] and hence allows to apply [9, Theorem 2.6]. This ensures the existence of a unique limit semigroup $\mathcal{E}$ generated by $F$ and being first order tangent to $F$ in the sense of the tangency condition (6). With this key observation, the remaining properties stated in the theorem directly follow from standard results, see, e.g., the proof of Theorem 2.3 in [10]. ∎

Acknowledgments

This work was supported by the German Research Foundation within the collaborative research center TRR154 “Mathematical Modeling, Simulation and Optimization Using the Example of Gas Networks” (DFG-SFB TRR154/2-2018, TP B01).