NUMERICAL NETWORK MODELS AND ENTROPY PRINCIPLES FOR ISOTHERMAL JUNCTION FLOW

We numerically explore network models which are derived for the isothermal Euler equations. Previously we proved the existence and uniqueness of solutions to the generalized Riemann problem at a junction under the conditions of monotone momentum related coupling constant and equal crosssectional areas for all connected pipe sections. In the present paper we extend this proof to the case of pipe sections of different cross-sectional areas. We describe a numerical implementation of the network models, where the flow in each pipe section is calculated using a classical high-resolution Roe scheme. We propose a numerical treatment of the boundary conditions at the pipe-junction interface, consistent with the coupling conditions. In particular, mass is exactly conserved across the junction. Numerical results are provided for two different network configurations and for three different network models. Mechanical energy considerations are applied in order to evaluate the results in terms of physical soundness. Analytical predictions for junctions connecting three pipe sections are verified for both network configurations. Long term behaviour of physical and unphysical solutions are presented and compared, and the impact of having pipes with different cross-sectional area is shown.

1. Introduction. Models for networks of hyperbolic conservation laws are for example used to describe traffic flow [6,15] and fluid flow in junctions [1,2,3,7,10,13,16]. Gas transport in pipe networks is the most common fluid flow application, but systems of water distribution, tunnel aerodynamics, mine ventilation and gas dynamics in engines and engine manifolds have also been considered [16]. With the extension of the theory to two-phase flows, at present developed for the drift-flux model [3,4], one could as well potentially investigate flow distribution issues in e.g. heat exchangers. More advanced models would however be needed in order to account for phase change in the fluid and heat transfer to the flows.
1.1. The network models. A network model for fluid flow in pipes consists of a set of hyperbolic equations describing the flow in each pipe together with a set of coupling conditions. The hyperbolic conservation laws are given by

GUNHILD A. REIGSTAD
where in each pipe section k, we seek the solution U k (x, t) for t ∈ R + and x ∈ R + . It is assumed that the junction-pipe interface is positioned at x = 0. The boundary conditions at each junction-pipe interface are determined by a set of coupling conditions. The conditions relate the flows in each of the N pipe sections and are closely linked to the applied conservation law (1). The coupling conditions we will be concerned with in this paper, derived for the isothermal Euler equations, are presented below in Definition 1.1. Following the approach of [8,9], this definition of the solution to the generalized Riemann problem at the junction is set up allowing for pipe-sections of different cross-sectional areas. The coupling conditions are denoted RP1 and RP2. Definition 1.1. A solution to the generalized Riemann problem at a junction, is a set of self-similar functions U k (x, t) such that RP0: For all k ∈ {1, . . . , N }, there exists a state such that U k (x, t) is given by the restriction to x ∈ R + of the Lax solution to the standard Riemann problem for x ∈ R: RP1: Mass is conserved at the junction: RP2: There is a unique, scalar momentum related coupling constant at the junction: Equation (4) is a constructed half-Riemann problem where U * k is calculated using the coupling conditions (5) and (6). The two statesŪ k and U * k are assumed to be connected by a wave with non-negative speed. For the isothermal Euler equations at subsonic conditions this implies that the states are connected by waves of the second family. A k is the cross-sectional area of pipe k, ρ * k and v * k are density and velocity at the boundary, respectively.
To select entropic solutions we follow the approach of [7,22]: where ρ 0 is some reference density and a is the constant fluid speed of sound, see (10).
well assumed that no vacuum states would appear. The problem was set up as a half-Riemann problem where the solution was found as a maximisation problem on the total mass flux at the junction. Supply and demand functions were defined and used as an important part of this problem formulation. The additional restriction needed in order to complete the problem formulation was discussed in [2]. Two options were considered: equal pressure at the vertex and subsonic flow on all outgoing pipes. The choice of equal pressure was used in [1]. The models in [1,2] were numerically implemented, and results from a set of test cases were presented.
A comparison between the network model proposed in [1] and results from a two-dimensional simulation model was performed by Herty and Seaïd [14]. A flow merge and a flow split configuration in a 90 • T-junction were used as test cases. It was found that in both cases the trend between the results were similar. However, for the flow merge case, deviations were seen. Therefore the approach of using empirical geometry dependent pressure loss coefficients was said to be supported for this configuration.
Colombo and Garavello [7] modified the mathematical description of network models for gas flow in pipe networks and established a network model for the psystem. Flow directions in the pipe sections were no longer given a priori, and the maximisation problem was replaced by an entropy condition (Eq. (7)) in addition to the mass-and momentum related coupling conditions (Eq. (5) and (6)). The selection of momentum related coupling constant was investigated by looking at the well posedness of the Riemann problem at a junction connecting three pipes. The analysis showed that equal momentum flux gave a problem that continuously depended on the initial conditions. This was not the case when pressure was used as coupling constant.
For the network model with momentum flux as coupling constant, Colombo and Garavello [7] investigated the existence and uniqueness of solutions, and proved that it was global in time and local in the subsonic state space.
The Cauchy problem at a junction for the p-system was investigated by Colombo and Garavello [8]. Well posedness was proved for initial states that are perturbations of subsonic states that are stationary, entropic solutions to the generalized Riemann problem.
Numerical results for network models based on the p-system are presented in [9]. In particular, the pressure is modelled by the γ-law: with the parameters γ = 1.4, p 0 = 1 and ρ 0 = 1. Three different junctions were considered. They consisted of two, three and four pipes connected at the junction, respectively. In all cases the pipes connected at the junction are of different cross-sectional area. In the last case the two momentum related coupling constants considered, pressure and momentum flux, were seen to give qualitatively different results for the chosen initial data. Brouwer et al. [5] used a friction dominated pipe model together with suitable coupling conditions in order to simulate a gas pipeline case from the Norwegian continental shelf. The pipeline network is characterised by supply and demand restrictions at a given set of locations, and the results are presented in terms of stationary pressure distribution and throughput. The simulation results were compared to previous presented results on the same case. However, due to missing case information the comparison was rather limited.
The analysis performed in our previous work [22] was motivated by the usage of two different momentum related coupling constants; pressure and momentum flux. Network models with the two coupling constants were considered for the special case where all pipe sections have equal cross-sectional areas. Under this assumption, existence and uniqueness of entropic solutions was found to be global in time and local in the subsonic state space: Proposition 1. (Proposition 9 in [22]) Consider the generalized Riemann problem at a junction of N pipes with equal cross-sectional areas. With pressure or momentum flux as coupling constant (RP2), there exists a unique solution satisfying RP0-RP2 provided that the initial data belongs to the subsonic region as defined in [22]. There does not exist solutions that satisfy RP3 (entropic solutions) for all initial data in this subsonic region.
Using the concept of ideal flow in a junction as starting point, the Bernoulli invariant was suggested as momentum related coupling constant. The invariant is derived from the conservation equation for total mechanical energy in smooth flows and is constant along streamlines. These properties are both valid for ideal flow, as such flow should be reversible and have uninterrupted streamlines.
The subsequent analysis of the network model containing Bernoulli invariant as coupling constant allowed us to prove existence and uniqueness globally both in time and in the subsonic state space: Proposition 2. (Proposition 10 in [22]) Consider the generalized Riemann problem at a junction of N pipes with equal cross-sectional areas. With Bernoulli invariant as coupling constant (RP2) there exists a unique entropic solution satisfying RP0-RP3 provided that the initial data belongs to the subsonic region as defined in [22].

1.3.
Overview. The aim of the present work is twofold: to present an appropriate numerical method for junction flow and to explore the impact of the different momentum related coupling constants numerically. As network theory is a fairly young field of research, there does not seem to be any established benchmarks which could serve as standard test cases. Therefore, we have here constructed two basic network configurations for our numerical investigations. For both of them we are able to numerically evaluate the entropy condition (7). This allows us to analyse if a solution is physical or not, even at conditions where we are not able to compare our results to analytical ones.
Before we describe the numerical method and the evaluation of results, two different issues are investigated analytically. First, the previous work performed in [22] was based on the assumption of pipe-sections of equal cross-sectional areas. In Section 2 we consider problems that are not restricted by this assumption. For monotone momentum related coupling constants, we provide a proof for the existence and uniqueness of solutions to the generalized Riemann problem in the sense of Definition 1.1.
Second, the standard Riemann problem is a well known test case with analytic solution. Formulated as a network problem, the case consists of two pipe sections of equal cross-sectional area connected at a junction. An analytical investigation of such a network configuration is presented in Section 3. The solutions of the network model and the standard Riemann problem are compared. Two conditions on the momentum related coupling constants, which ensure that the two solutions coincide, are identified. Both pressure, momentum flux and Bernoulli invariant as coupling constant are seen to fulfil the conditions. Thus the corresponding models should all predict the analytical solution to the standard Riemann problem.
The implemented network models are presented in Section 4. The flow in each pipe section is calculated using an approximate Riemann solver of Roe, presented in Section 4.1 and Section 4.2. Section 4.3 shows how the conservative flux function at the pipe junction interface has been evaluated and set such that mass is conserved at the junction. Finally, the numerical procedure used to calculate U * k for each pipe section is presented in Section 4.4.
Numerical results for the first test case, a closed network consisting of three pipe sections connected at two junctions, is presented in Section 5.1. Using the total energy in the network as indicator, we show how both pressure and momentum flux as coupling constant lead to unphysical solutions for different initial conditions. We also show the long term development of the total energy as waves interact in pipe sections and at the junctions.
Steady state conditions in an open network model is presented in Section 5.2. The network consists of a symmetric branch structure where the flow in each pipe section is divided equally between two pipes at each junction. Simulations on a model where all pipes are of equal cross-sectional area verify the presence of unphysical solutions for pressure and momentum flux as coupling constant. Using the same initial conditions, but changing the cross-sectional areas such that the velocities in all pipe sections are equal, we observe no unphysical behaviour.
1.4. The isothermal Euler equations. In this work, we follow the approach of [1,2,7,14,17,19] and consider one dimensional, compressible, pipe flow governed by the isothermal Euler equations. These consist of the isentropic Euler equations together with a specific pressure law Here ρ and v are the fluid density and velocity, respectively, p(ρ) is the pressure and a is the constant fluid speed of sound. A more general formulation of (10) was considered in [7]. By an appropriate extension of the theory presented in [22], it should be possible to extend our numerical method to this more general formulation.

2.
Existence and uniqueness of solutions to the generalized Riemann problem at a junction. A proof of uniqueness and existence of solutions to the generealized Riemann problem, as stated in Definition 1.1, was presented in [22] for the special case of equal cross-sectional areas in all pipe sections. This proof is valid for monotone momentum related coupling constants as defined below: for both shock waves and rarefaction waves connecting the two states U * k andŪ k .
In Definition 2.1, M * k = v * k/a is the Mach number and H * k is the coupling constant value for pipe section k. The constant is a function of ρ * k or M * k alone, because both variables may be expressed as a function of the other variable and ofŪ k , which is a constant state, using the appropriate wave equation [22].
Here we will extend the proof of existence and uniqueness to the case of pipe sections of different cross-sectional areas. That is, we will show that the conditions RP0 -RP2 have a unique solution that exists for all subsonic initial data when the momentum related coupling constant is monotone. The proof will extensively be based on the results found in [22].
First we state the extended definition of subsonic initial conditions: Assume that a set {Ū k } of initial data is given. Assume that this set satisfies the conditions and Such a set of initial data is said to belong to the subsonic region.
R2 and S2 denote rarefaction and shock wave of the second family, respectively. The notation indicates which of the wave equations that is used in H * k (M * k ). Condition (1) is the obvious constraint of subsonic initial conditions in each pipe section. The two additional conditions set the limits for where the solutions are subsonic.
Then we may state the following: Proposition 3. Assume that the initial dataŪ k belongs to the subsonic region in the sense of Definition 2.2 and that the momentum related coupling constant H is monotone in the sense of Definition 2.1. Then there exists a unique set of subsonic solutions satisfying RP0-RP2.
Proof. Proposition 1 in [22] proves the uniqueness of a state U * k that satisfies RP0 for the conditions 1. The initial stateŪ k is subsonic, 2. The momentum related coupling constant is monotone in the sense of Definition 2.1 3. The coupling constant value is given and satisfies The proof of the proposition is independent of the assumption on the cross-sectional area of each pipe section, and is hence valid for the generalized Riemann problem in the sense of Definition 1.1. As the last condition is satisfied for a set of initial data that belongs to the subsonic region, all three conditions are met. The uniqueness of a state U * k satisfying RP0 is thus proved. Figure 1. A junction with 2 connected pipes Proposition 3 in [22] proves the uniqueness of the set of solutions {U * k } that satisfies RP0-RP2, given that the states U * k uniquely satisfies RP0 and that the set of solutions exist. Using the rarefaction and shock wave equations, it is proved that This implies that monotone coupling constants are as well monotone functions of the [ρM ] * k variable. As the cross sectional area of each pipe section is a constant, we may therefore conclude that the total mass function J (H), as defined in Equation (12), is a monotonically increasing function of a monotone coupling constant, H, in the sense of Definition 2.1. Thus Proposition 3 in [22] holds for the case of pipe sections with different cross-sectional areas.
As in [22], the definition of the subsonic region in Definition 2.2 guarantees the existence of the unique set of solutions.
3. The standard Riemann problem test case. In addition to the monotonicity constraint (11), two requirements on the momentum related coupling constant are found by looking at the standard Riemann problem test case. In this test case two pipe sections of equal cross-sectional area are connected at a junction along the same axis, see Figure 1. Thus, the solution predicted by the network model must coincide with that of a Standard Riemann problem. The initial conditions in two such pipe sections are given asŪ 1 andŪ 2 . Their boundary conditions at the pipe junction interfaces, Eq. (3), are U * 1 and U * 2 , respectively. The corresponding standard Riemann problem would have ρ L =ρ 1 , v L = −v 1 and U R =Ū 2 with a resulting intermediate state U m . In the following we begin by identifying a set of {(ρ * k , v * k ), k = 1, 2}, which correctly describes the intermediate state. and That is: Proof. The solutions to the standard Riemann problem is characterised by rarefaction and shock waves of the 1 st and 2 nd family (denoted R1, R2, S1 and S2): In the network model, the boundary conditions are found from the equations: Equation (25) to (28) are identical to Equation (21) to (24) if Equation (17) and (18) hold. Consequently Equation (19) and (20) are fulfilled.
Next, we consider the momentum related coupling constant of the network model and identify two constraints: Lemma 3.2. Consider isothermal flow. The coupling conditions of a network model for two pipe sections of equal cross-sectional area reduces to the following set of equations: Equation (17) and (18) form a unique solution to this set if and only if Proof. The existence of a unique relation between ρ * 1 and ρ * 2 is guaranteed by the monotonicity constraint (32). Symmetry in the ρ * k v * k variable is required for Equation (17) to be a solution of the Equations (29) and (30). If one or both of the conditions are not met, Equation (17) and (18)  Proof. Lemma 3.1 proves that if the network model predicts equal densities at the pipe-junction interfaces and velocities that are equal in absolute value and of opposite sign, then the solution will coincide with that of the standard Riemann problem. Lemma 3.2 proves that the network model will predict this kind of densities and velocities if and only if the two constraints on the momentum related coupling constant are met.
A simple analysis, omitted here, show that pressure and Bernoulli invariant as momentum related coupling constant all fulfil the conditions set in Lemma 3.2.
Hence they all correctly reproduce the solutions to the standard Riemann problem. Observe however that in this case, the entropy condition, RP3, is trivially satisfied. Therefore, studies involving junctions with more than two pipe sections are required if we wish to gain insight into the general entropy consistency of the various coupling constants.
4. Numerical implementation. In this Section, we describe a numerical implementation of the network model. The approximate Riemann solver of Roe has been chosen, as it allows for easy extension to second-order accuracy for smooth solutions using the wave-limiter approach of LeVeque [18]. This approach is described in more detail in Section 4.2. In Sections 4.3 and 4.4, we describe the boundary treatment needed to ensure that the coupling conditions are satisfied.
4.1. The approximate Riemann solver of Roe. The set of equations governing the flow in each pipe is solved using the approximate Riemann solver of Roe [23]. The approximation is made by replacing the nonlinear problem (1) by a linearised problem at each cell interface [18]: The Roe matrix, A i−1/2 , must fulfil three conditions : Condition 1. A i−1/2 must be diagonalisable and have real eigenvalues The conditions are set in order to ensure a hyperbolic system (Condition 1) that is consistent with the original conservation law (Condition 2). In addition W p will be an eigenvector of A i−1/2 if the states Q i−1 and Q i are connected by a single wave, W p = Q i − Q i−1 , in the true Riemann solution (Condition 3). The last condition also guarantees that the selected scheme (39) is conservative [18].
For the isothermal Euler equations, the following matrix fulfils the conditions [24, Eq. (11.42) and Eq. (11.43)]: 4.2. Numerical algorithm. Equation (39) describes the numerical algorithm, written on conservation form. The algorithm is based on Godunov's method, applied for non-linear systems. It contains both a conservative Roe flux, F , and a correction term, F , where the latter extends the algorithm to a high-resolutions method [18].
The conservative Roe flux is found using the Roe matrix A i−1/2 : The matrix | A i−1/2 | is defined as [18]: where andλ i−1/2 and R i−1/2 are the eigenvalues and right eigenvector matrix of A i−1/2 , respectively. The high-resolution correction terms are defined by Equation (44), where the limited wave W p i−1/2 is given by Equation (45). In this equation, each wave of the solution is modified by a limiter function φ(θ).
In our implementation, the monotonized central-difference limiter function (MC limiter) was chosen:

GUNHILD A. REIGSTAD
The ratio θ is found by: where the index I is defined as: The waves, W p i−1/2 , may be calculated from the eigenvectors of the Roe matrix: wherer p i−1/2 is the right eigenvector of A i−1/2 belonging to family p and The wave speeds needed in Equation (44) and Equation (48) are defined as the eigenvalues belonging to the Roe matrix: 4.3. The flux function at the junction boundary. In our numerical model, all boundary conditions are managed using the ghost cell approach. We have chosen to extend the computational domain by one ghost cell at each end of a pipe section. At a boundary towards a junction, the boundary condition found from the network model, U * k , is used to determine the ghost cell variables. By doing so, some special care must be taken about axis directions and the calculation of the flux across the boundary.
We denote the condition in cell i in the k-th pipe section at time t = n∆t by Q n i,k and the ghost cell by Q n GC,k . The ghost cell variables are updated according to Equation (52).
The variable j k is used as an indicator of the relation between the global axis, i x,g , used in the numerical model and the local axis, i x,l , used by the network model at each junction: When calculating the numerical flux across the pipe-junction boundary, it is necessary to ensure that the mass conservation coupling condition (5) is fulfilled for each junction. The selected approach is outlined below in Proposition 5.
Proposition 5. Denote the internal grid cell which is closest to the junction boundary, Q n SJ,k . The numerical network model will conserve the total mass if the upwind Godunov flux (54) is used at pipe junction interfaces: where F tot (Q n GC,k , Q n SJ,k ) = F (Q n GC,k , Q n SJ,k ) + F (Q n GC,k , Q n SJ,k ).
Proof. The numerical mass flux at a pipe-junction interface may be written as: where V [h] is used as notation for the h-th element of the vector V . Numerical conservation of mass at a junction is then expressed as: Using the selected flux (54), Equation (57) becomes: For a junction connecting N pipes this is equal to: which satisfies Equation (5).
Similarly, at the right hand side of the domain the correction term is: We may use Condition 3, presented in Section 4.1, and simplify Equation (60) and (61). Noticing that A i−1/2 may be written as: and inserting the expression for | A i−1/2 | found in Equation (41), Equation (60) may be written as: Similarly, Equation (61) may be written as: 4.4. Calculation of boundary conditions at pipe junction interfaces. The implemented procedure calculating the boundary conditions is in accordance with the results in [22]. We consider the three different momentum related coupling constants; pressure (33), momentum flux (34) and Bernoulli invariant (35), which are all monotone in the sense of Definition 2.1 [22]. The inverted functions, density and velocity as function of the momentum related coupling constant, will depend on the appropriate wave equation. They are denoted (ρ) R ,(M ) R and (ρ) S ,(M ) S when rarefaction and shock wave equations are used, respectively. For a given coupling constant value,H, the density, Mach number and mass flux at a pipe-junction interface may be calculated: Numerically, the inverted function values are found by iteration on the function Here the appropriate wave equation is used to describe the relation between the conditions of the ghost cell, Q n GC,k = [ρ * k , ρ * k v * k ] T , and of the nearby internal cell, Q n SJ,k = [ρ k ,ρ kvk ] T , which is regarded as initial condition of the pipe section. A Newton-Raphson algorithm is used to solve (66), utilising that the derivative dH * k/dM * k is known. If the algorithm fails, a bisection method is used on the interval M * k ∈ [−1, 1]. The total mass flux function is defined by Equation (12). Using this equation together with the inverted mass flux function (65), we calculate J (H − ) and J (H + ), see Equation (13) and (14). If the set of initial states, {Q n SJ,k }, belong to the subsonic region in the sense of Definition 2.2, a solution may be calculated. This is done with the aid of a bisection method on the intervalH ∈ [H − , H + ].
5. Numerical results. The impact of the selected momentum related coupling constant will be explored using two different network configurations. The first configuration, a closed system consisting of three pipes and two junctions, will be used to show how the selected constant affects the energy of the system. The second configuration, an open system of 16 pipes connected in a network, will be used to show how the coupling constant influences the steady state flow situation in the network.
In all cases, the fluid speed of sound is set to a = 300 m/s.

5.1.
Simulations on a closed network. An outline of the closed network configuration is shown in Figure 2. Three pipe sections, labelled S1 to S3, each of length L = 50 m, are connected by two junctions. All pipe sections are assumed to have equal cross-sectional areas. The global axis direction, i g , is set in the direction from junction J1 to junction J2, as indicated in the figure. It should be noted that for junction J2, the global axis and the local axis, i l , that is applied by the network theory are of opposite directions, see Figure 2b. Initially, each pipe is filled with stagnant fluid of uniform pressure, with different pressure levels in each pipe. Two different sets of initial pressures have been selected. These are shown in Table 1.    Figure 3 and Figure 4. The simulations were run until T = 0.06 s. For simplicity, only results for the first pipe-section are presented, and only for pressure as momentum related coupling constant. The results for the other coupling constants and pipe sections are however similar. The analytical pressure-and velocity profiles, shown in the figures 3a, 3b, 4a and 4b are calculated using the wave equations for rarefaction and shock waves of the second family. If U * k andŪ k are connected by a rarefaction wave, Equation (67) and (68) describe the density and velocity profiles in the pipe section.
In these equations, λ 2 is the eigenvalue of the second family: Equation (71) and (72) describe the profiles when the two states are connected by a shock wave.
Comparing analytical and numerical results, it is seen that only the lowest grid resolution leads to a significant deviation.  The numerical profiles are seen to converge towards a solution as the grid size is reduced. The exact solution, that is, the limit which the numerical results should converge to, is shown as "Analytic 1". The profile is derived as described in Observation 1 to Observation 3 and Equation (85), and relies on the condition that the two initial waves entering each pipe section have not yet interacted. This condition is satisfied at T = 0.06 s. The "Analytic 2" profile is derived in order to distinguish between physical and unphysical solutions. We will describe this further later in this section. Figure 2 has an interface towards a junction in both ends. Denote the boundary conditions at these two pipejunction interfaces of pipe section k, as U b k,J1 and U b k,J2 . Note that both conditions relate to the global axis, i g . Then, for the closed configuration shown in Figure 2 the following holds:

Observation 1. Each of the pipe sections shown in
Proof. Figure 2 shows that the two junctions with connected pipes have the same initial conditions. Hence, However, as the directions of the global axis and the local axis used by the network theory are opposite at the second junction, J2, the velocities are related by Equation (74).
Observation 2. Denote boundary conditions as in Observation 1 and consider the closed configuration shown in Figure 2. A pipe section containing two rarefaction waves that has not yet interacted has the following energy-content as a function of time: Proof. Energy and energy flux is an entropy entropy-flux pair for the isothermal Euler equations. [11, p. 212] Thus the entropy condition for a pipe section is given by Equation (77), with energy, E, defined by Equation (78).
In the present case, the pipe section only contains rarefaction waves. For such a smooth solution Equation (77) shows that energy is conserved: Integration along the pipe section from the interface towards J1 to the interface towards J2 gives Using the results in Observation 1 and integrating from 0 to t gives the result in equation (76).
Observation 3. Denote boundary conditions as in Observation 1, the initial condition in the pipe asŪ k and consider the closed configuration shown in Figure 2. A pipe section containing two shock waves that have not yet interacted has the following energy-content as a function of time: where L is the length of the pipe section.
Proof. Denote the speed of the shock waves emerging from the interfaces towards J1 and J2 by s k,J1 and s k,J2 , respectively. Due to symmetry, the relation between the speeds may be found as: The final expression in (82) is obtained when taking into account that the initial velocity is zero, see Table 1. At a given time t, the two shock waves divide the pipe section into three regions. Closest to the pipe-junction interface the conditions are given by U b k,J1 and U b k,J2 , respectively. The mid-region is given byŪ k . Total energy is found by summation of energy in each region: where energy, E, is defined by Equation (78). Inserting the results from Observation 1 and Equation (82), we may write: Equation (81) is found by inserting the expressions forĒ k , E b k,J1 and s k,J1 (82) into Equation (84).
The total energy may then be found by summation: where R2 is the subset of pipes containing rarefaction waves, and S2 is the subset containing shock waves. The "Analytic 1" profile is compared with the numerical results (see Figure 3c and Figure 4c), in order to determine the appropriate grid size. The results show that only at the smallest grid size, ∆x = 5.0 × 10 −3 m, the deviation from the analytical profile is negligible. All further results presented for this configuration are therefore obtained using this grid size.
The second curve, denoted "Analytic 2", is used to determine if the numerically obtained solution is entropic in terms of Definition 1.2. These solutions are characterised by the absence of energy production in the junctions. In order to identify energy production or dissipation in the junctions, we need to account for the energy development in each of the pipe sections. Looking at the entropy condition in Equation (77) we see that pipe sections containing shock waves will have a reduced energy content as a function of time. This energy loss may be calculated as shown in Equation (86). By subtracting this loss from the total initial energy, as shown in Equation (87), we may calculate the total energy of the system for conditions with no energy dissipation or production in the junctions.
Entropic solutions are then characterised by: Otherwise there is energy production in the junctions.
We use this criterion (88) to investigate the numerical results presented in Figure 3c and Figure 4c. It is clearly seen that the network model with pressure as momentum related coupling constant yields an entropic solution only for the second set of initial data. For the first set of initial data, the model predicts a condition with energy production in the junctions.
Next, we present results for all three momentum related coupling constants. Figure 5 and Figure 7 show the pressure and velocity profiles in each of the pipe sections at T = 0.06 s. The corresponding energy profiles are shown in Figure 6 and Figure 8. As expected, only the choice of Bernoulli invariant as momentum related coupling constant leads to energy conservation at the junctions. The duality between the two momentum related constants, pressure and momentum flux, is as well as expected [22]. For the first set of initial data, pressure as constant results in energy production at the junctions, while equal momentum flux gives energy dissipation. For the second set of initial data this is the other way around.

Development of total energy as a function of time.
A second set of simulations were run until T = 1.00 s. Within this time-frame the first set of waves entering each pipe section interact, the resulting waves propagate out of the pipe sections and into the junctions, new waves are created and re-enter each pipe section and so forth. The resulting energy profiles for the two selected cases are shown in Figure 9. As seen, there is a net reduction of total energy for all three momentum related coupling constants. However, the profiles for pressure and momentum flux as coupling constant are not decreasing monotonically. This is due to the direction of the pipe flows. For instance we may consider the profile for pressure as coupling constant in Figure 9a. From T = 0.00 s to T 0.17 s there is an increase in total energy. At T 0.17 s there is a significant change in the profile and the total energy is decreasing. At T = 0.00 s the first set of waves enter the pipe sections. Relative to the local axis at the junctions, two pipe sections have boundary conditions with negative flow velocity, one pipe section has positive boundary velocity. As seen in Figure 3c, this implies energy production for a network model having pressure as coupling constant. At T 0.17 s the second set of waves enter the pipe sections. In this set, the flow velocities are of opposite sign compared to the first set. Thus, two pipe sections have boundary conditions with positive flow velocity, one pipe section has negative boundary velocity. For the selected coupling constant this implies energy dissipation at the junctions, which explains the significant change from increasing to decreasing total energy seen in Figure 9. The source is connected at the inlet of pipe section S1 and holds a constant pressure of 2.0 bar. At T = 0.0 s, a membrane separating the source from the pipe is broken, and fluid is entering pipe section S1.
At T = 10.0 s, a stationary flow through the network is established. The flow conditions are evaluated using an integral analysis within a set of fixed boundaries. For any boundary the following holds if mass and energy is conserved: where S i and S o are the sets of pipe sections with flow direction into and out of the fixed boundary, respectively. In our analysis we use four different boundaries, indicated as CV 1 to CV 4 in Figure 10. For each boundary we calculate the deviances: based on densities, ρ k and velocities, v k that are averaged over all grid cells of each pipe section. First we assume that the cross-sectional areas of all pipe sections are equal. The corresponding results are presented in Table 2. As seen, M ratio deviates from zero. This is most probably due to minor deviations from steady state in one or more of the pipe sections. However, the effect of energy production and dissipation at the junctions is clearly seen from E ratio . There is a significant difference between M ratio and E ratio for the boundaries CV 1 , CV 2 and CV 3 for pressure and momentum flux as coupling constant. In the case of pressure as constant, there is energy dissipation at the junctions J1 to J7. In the case of momentum flux as constant, energy is produced in the same junctions.
Looking at the last boundary, CV 4 , we observe that the two deviances M ratio and E ratio are approximately equal for all three momentum related coupling constants. At this boundary the mass balance (89), the symmetry property of the momentum related coupling constant (31) and the application of the coupling constant (6) on loss-less pipe flow result in the following set of equations: test case of a standard Riemann problem. Hence, for the CV 4 boundary, the three different momentum related coupling constants all yield entropic solutions. The steady state pressure and velocity profiles for the pipe-sections are showed in Figure 11.
Next we change the assumption of pipe sections of equal cross-sectional areas to the following: Meaning that at each junction, except J8, the cross sectional area of each of the outgoing pipe sections is half of the cross sectional area of the ingoing pipe.
The corresponding results are shown in Table 3. For this network there is no significant difference between the three momentum related constants and the deviances are negligible. An analysis on a junction connecting three pipe sections, similarly as in [22], shows why.
Junction J1 with the connected pipe sections S1, S2 and S3 is representative for the junctions J1 to J7. These are the junctions impacting the results for CV 1  Figure 9. Energy profiles for the three different momentum related coupling constants. Note that only Bernoulli invariant as momentum related coupling constant,H = ln(ρ) + 1 /2M 2 , yields a monotonically decreasing profile, as explained in Section 5.1.2.
to CV 3 . Due to the symmetry in the branches of the network, the velocities in the pipe sections S2 and S3 are equal. Observation 4 and Observation 5 shows that network models using pressure or momentum flux as coupling constant yields energy conservation in the junctions due to the selected cross-sectional area ratios (95)-(96).

Observation 4.
A network model using pressure as momentum related coupling constant (33) yields energy conservation in the junctions of Figure 10 for pipe sections with cross-sectional areas as described by the Equations (95) and (96).
Proof. We consider junction J1 with the connected pipe sections S1, S2 and S3. Mass conservation at the junction (5) may be expressed as: Due to symmetry in the network, From the final expression in (100) we see that energy is conserved in the junction for v * 1 = 0 and for A 2 = 1 /2A 1 . As seen from Equation (96), the last condition for energy conservation is fulfilled.
Observation 5. A network model using momentum flux as momentum related coupling constant (34) yields energy conservation in the junctions of Figure 10 for pipe sections with cross-sectional areas as described in the Equations (95) and (96).
Proof. We consider junction J1 with the connected pipe sections S1, S2 and S3. Mass conservation at the junction (5) may be expressed as: The entropy condition (7) becomes: Energy conservation,Q = 0, is obtained for two different conditions; z 1 = 0 and b 1 (M * 1 ) = b 2 (M * 2 ). The first condition implies that M * 1 = 0. The second condition is fulfilled for M * 1 = M * 2 and M * 1 = −M * 2 . In our case, the physically relevant solution is M * 1 = −M * 2 . From Equation (102) we see that this corresponds to the condition A 1 = 2A 2 , which is fulfilled by Equation (96).
Observation 4 and Observation 5 show that the velocities in each of the pipe sections are equal, if we refer to the global axis directed from the source to the outlet of pipe section 16. This is also the situation for the network model using Bernoulli invariant as coupling constant. Looking at the expressions for momentum flux and Bernoulli invariant as coupling constant, Equation (34) and (35), we see that this implies equal density, and thus equal pressure, too. For the present case the resulting pressure is p = 1.7305 bar. The velocity is v = 43.468 m/s. 6. Summary. We investigate solutions to the generalized Riemann problem for the isothermal Euler equations analytically and numerically. Restricted to monotone coupling constants and pipe sections of equal cross-sectional area, a proof of existence and uniqueness of such solutions was given in our previous work [22]. In the present paper we have expanded this proof to the case of pipe sections of different cross-sectional areas.
We have also compared the solution of a network model that describes a standard Riemann problem to the analytical solution of such a problem. The analysis lead to two constraints on the momentum related coupling constant; it must be symmetric in the ρv variable and monotone in the ρ variable.
A classical high-resolution Roe scheme is used to simulate the flow in each pipe section. For the pipe-junction interfaces we have proposed a numerical method which ensures conservation of mass at the junctions. In our implementation, boundary conditions are generally managed using the ghost cell approach. At pipejunction interfaces the ghost cell variables are thus set equal to the boundary condition predicted by the network model, U * k . Conservation of mass at each junction is obtained by using the upwind Godunov flux as numerical flux function at the pipe-junction boundaries.
Numerical results are presented for two different network configurations; a closed system consisting of three pipe sections connected at two junctions and an open system consisting of pipe sections in a symmetrical branch-tree structure. The results are analysed with respect to entropy consistency. That is, the evaluation is performed by applying the entropy condition in order to determine if a solution is physically reasonable or not. This is done in order to evaluate the three different options for momentum related coupling constant used in the network model.
We present long-term simulation results, which show the impact of unphysical solutions predicted by the network models. For the closed system, this is seen as an increase in total energy at certain time intervals. For the open system, steady state energy fluxes are used to identify energy production in junctions.
Analytical predictions of entropy consistency for junctions connecting three pipe sections of equal cross-sectional areas were provided in [22]. Numerical results for such junctions, derived in the present paper, are in accordance with these predictions. They show that the options of pressure and momentum flux yields unphysical solutions for certain initial data. Only Bernoulli invariant as momentum related coupling constant yields entropic solutions for all the presented test cases.