Well-balanced scheme for gas-flow in pipeline networks

1.   Introduction

The study of mathematical models for gas flow in pipe networks has recently gained interest in the mathematical community, see e.g. [3,4,9,12,27]. While in the engineering literature [31] the topic has been discussed some decades ago, a complete mathematical theory has only emerged recently, see e.g. [16] for the Euler system on networks, [14] for the p–system on networks and [8] for a recent review article on general mathematical models on networks. Depending on the scale of phenomena of interest, different mathematical models for gas flow might be useful. A complete hierarchy of fluid–dynamic models has been developed and discussed in [9]. Therein, typical flow rates and pressure conditions are given and it is shown that a steady state algebraic model can be sufficient to describe average states in a gas network. Models based on an asymptotic expansion of the pressure may lead to further improvements in case of typical, slowly varying, temporal flow patterns [23,18]. If a finer resolution of the spatial and temporal dynamics is required, the isothermal Euler equations (1) provide a suitable model [3]. Here we focus on schemes which capture both steady states and small temporal and spatial perturbations.

Schemes which preserve a steady state exactly are called well-balanced, and their development is a lively topic in the field of hyperbolic balance laws, see e.g. the monograph [32]. Usually, these schemes use specific knowledge of an equilibrium state. As a consequence, well-balanced schemes for still water such as [1,29], or for moving water, [30], or for wet-dry fronts such as [6,10], which all approximate solutions to the shallow water equations, use different discretization techniques. A unified approach to well-balancing in one space dimension was recently proposed by Chertock, Herty and Özcan [11], who integrate the source term in space and subtract it from the numerical flux. This results in a new set of so-called equilibrium variables. The equilibrium variables are then reconstructed, and their values at the cell interfaces are used to compute numerical fluxes, which involve a new equilibrium limiter. Chertock et al. demonstrated the feasibility of this approach using a second order scheme for subsonic gas flow in a pipe with wall friction.

To the best of our knowledge, there is currently no well-balanced scheme for hyperbolic flows on networks. Here spurious oscillations may not only be caused by an imbalance of numerical fluxes and source terms, but also by discretization errors at junctions and compressors. In the present paper, we extend the equilibrium variable approach and develop a second order well-balanced scheme on a network. We also simplify the numerical fluxes introduced in [11] by removing the equilibrium limiter. High-order schemes for gas networks have been introduced only recently in [2,7,28]. A challenging question would be to extend our new well-balancing method to these more accurate schemes.

In the following we introduce the mathematical model for the temporal and spatial dynamics of gas flow in pipe networks. For simplicity, we study a single node $ x = x_o $ where $ M $ pipes meet. The flow within each pipe $ i = 1 \ldots M $ is governed by the isothermal Euler equations

with conservative variables $ U_i $, flux $ F(U_i) $ and source $ S(U_i) $ given by

Here $ \rho_i, q_i $ and $ p(\rho_i) $ are the density, momentum, and pressure of the gas, $ f_{g, i} $ is the friction factor and $ D_i $ the diameter of pipe $ i $. The pressure of the gas is given by Gamma-law,

For the numerical tests we focus on the special case of isothermal pressure,

with speed of sound $ a>0 $. However the well-balanced scheme discussed in the following sections, including the proofs of Theorems 3.1 and 4.1 are valid for any pressure given by Gamma-law (3). We complete (1) with initial conditions within and boundary conditions at the ends of the pipes. The boundary conditions at a node of multiple pipes are implicitly given by coupling conditions [3,35], which take the form of $ M $ nonlinear algebraic equations involving traces $ U_i $ of the conserved variables at node $ x_o $. We write them in the general form

In [12,13,35] general conditions are identified which guarantee a well–posed problem for initial data with suitable small total variation. If at a node the conservation of mass and equality of adjacent pressures is assumed, then existence, uniqueness and continuous dependence on the initial data for the p–system was shown in [15].

For pipes with a junction, a flow which is steady within each pipe, and fulfills the coupling conditions at the junction is called a steady state [22]. The coupling conditions are further detailed in Sections 2 and 3.

2.   Review of coupling conditions for the p–system

Our construction of well-balanced schemes is based upon analytical results for the isothermal Euler equations which have already been established (see e.g. [3,15] and the references therein). For completeness, and to fix the notation, we would like to give a brief summary.

When we study the coupling condition, we set the source terms $ S(U_i) $ to zero. This is based on the heuristic assumption that wall friction can be neglected at the instance of interaction at the node. It can also be justified rigorously for the semi-discrete scheme (38). The eigenvalues for the homogeneous $ 2\times 2 $ system (1) are $ \lambda_1 = \frac{q}{\rho}-a $ and $ \lambda_2 = \frac{q}{\rho}+a $. We assume that all states are subsonic, i.e.,

This assumption is satisfied for typical gas flow conditions in high–pressure gas transmission systems. We denote the set of all incoming (respectively outgoing) pipes by $ I^- $ (respectively $ I^+ $). For $ i \in I^- $, we parametrize the incoming pipes by $ x \in \Omega_i : = (-\infty, x_o). $ Similarly, we parametrize outgoing pipes $ j \in I^+ $ by $ x \in \Omega_j : = (x_o, \infty). $ Let us fix a time $ t_o \ge 0 $. It is important to note that there are $ 2M $ different traces at the node $ (x_o, t_o) $. We denote them by,

Note that the limits are exchanged in (8) and (7) (respectively (10) and (9)), so $ U_{{i}}^{{o}} $ and $ U_{{j}}^{{o}} $ are limits along the $ x $-axis. We call them the old traces. The states $ U_{{i}}^{{*}} $ and $ U_{{j}}^{{*}} $ are limits along the $ t $ axis, and we call them the new traces (see Figure 1).

The construction of the coupling conditions starts by connecting, within each pipe, the old with the new trace by a Lax curve entering the pipe. For an incoming pipe, this will be a curve $ \overline{U}_i(\sigma_i) $ of the first family with left state $ U_l : = U_{{i}}^{{o}} $. According to [17], this curve for $ \gamma = 1 $ is given by

Analogously, for an outgoing pipe, we use curves of the second family with right state $ U_r : = U_{{j}}^{{o}} $, which are given by

We refer to [17] for case $ \gamma>1 $. The Lax curves 1-R and 1-S(respectively 2-R and 2-S) have $ C^2 $ continuity at the point $ U_l $(respectively $ U_r $).

The parameters $ \sigma_i, i = 1 \ldots M $ will be determined from the $ M $ coupling conditions

Now we set

By construction, the new traces satisfy the coupling conditions (5). Note that the number of unknowns in (13) are halved compared to (5) as each conservative variable, $ U_i\in \mathbb R^2 $ can be written in terms of single parameter $ \sigma_i \in \mathbb R $.

So far, we have reviewed the general framework which was established and applied in [3,15,33]. In the following we focus on a particular coupling condition for which we design a well-balanced scheme. Let $ A_i = \frac{\pi}{4}D_i^2 $ be the area of the cross section of pipe $ i $. The default coupling condition requires that the total incoming mass flux at each node $ x_o $ equals the total outgoing mass flux,

since mass should not be accumulated or lost at the junction. Various approaches have been studied in order to model the other $ (M -1) $ coupling conditions. A seemingly obvious choice would be conservation of momentum. However as momentum is a vector quantity it is difficult to describe the conservation of momentum in a one-dimensional model as the junctions of the pipe are three-dimensional. Multi-dimensional approaches considering a 2D node for a 1D flow have been discussed in [24,5]. Coupling conditions based on enthalpy have also been studied in [3,27,35,34]. In general these coupling conditions can be represented in the form,

In [35], Reigstad et al. showed that the coupling conditions are well-posed if the function $ h $ is monotonic with respect to the parameter $ \sigma $ of the Lax curves,

Our analytical results are obtained under the general coupling condition (16) (respectively (16)), together with the monotonicity condition (17). In the numerical experiments, we use the following pressure-coupling condition: we require that the traces of the pressures should take one and the same value $ p^* = p^*(t_o) $ at the $ t $-axis,

This condition is common in the engineering literature, and for a certain regime it is indeed a suitable approximation of the two-dimensional situation [24].

Thus the coupling function (5) at a junction reads

We now turn to a junction which models a compressor between the incoming pipe $ i = 1 $ and the outgoing pipe $ i = 2 $, both of the same diameter. The coupling conditions for the new traces are

Here $ C\!R \ge 1 $ is the compression ratio. It is usually time–dependent, and we consider it to be a given, external quantity. Thus the coupling function for the compressor becomes

Summarizing, the analytical problem at the nodes is to connect the old traces $ U_{{i}}^{{o}} $ within each pipe to a new trace $ U_{{i}}^{{*}} $ along the incoming Lax curve in such a way that the new traces satisfy the coupling conditions across the node. It was proven in [20,12,21] that this problem has a unique solution.

If the old traces are subsonic, and their variation is small enough, then the new traces will be subsonic as well. The new traces serve as initial data in the Riemann solver which determines the numerical flux.

In this paper, we show the analytical results of well-balancing for the case of multiple pipes meeting at the junction. However the scheme is valid for the case of compressor as well, as can be seen from the numerical results in Section 5.2.

3.   Coupling conditions in terms of equilibrium variables

The difficulty in preserving steady states is that the divergence of the conservative fluxes is approximated by a flux-difference, while the source is usually integrated by a quadrature over the cell. If this is not tuned carefully, the equilibrium state is not maintained, and spurious oscillations may be created. Chertock, Herty and Özcan [11] resolved this difficulty for one-dimensional balance laws by integrating the source term and hence writing it in conservative form. They applied this approach to the Cauchy problem for $ 2\times 2 $ balance laws. Here we extend their method to a node in a network. Equation (1) can be stated as

where the flux variable,

and fluxes $ K_i, L_i $ and an integrated source term $ R_i $ is given by

The point $ \tilde{x_i} $ belongs to $ \overline{\Omega}_i $ and is arbitrary but fixed. Later on, we choose $ \tilde{x_i} = x_o $ for all $ i $. We call $ (K, L) $ the equilibrium variables, since they are constant for steady states. Let us consider the general pressure law (3). Given the integrated source term $ R_i $ and the equilibrium variables $ V_i = (K_i, L_i) $ we now solve equation (24) for the conservative variables $ (\rho_i, q_i) $. Omitting the subscripts of the pipes in the following calculation, we rewrite (24) as

In other words, our task is to find the non-negative real roots of

with

For $ \rho\ge0 $ the function $ \varphi $ is convex, has roots $ \rho_0 = 0 $ and $ \rho_1 = b^{1/\gamma} $ and a critical point in

with corresponding minimum value

An elementary calculation shows that if $ K, L $ and $ R $ belong to a sonic state $ (\rho_s, u_s) $, then $ \rho_* = \rho_s $. Thus, if $ c = \varphi_* $, then $ \rho_* $ is a double root of (25). Else, if $ 0 \geq c > \varphi_* $, we have two real roots $ \rho_\pm $, with

We find the supersonic root $ \rho_- $ by a Newton iteration with initial value $ \rho_0 $, and the subsonic root $ \rho_+ $ using the initial value $ \rho_1 $. By convexity of $ \varphi $, the convergence is quadratic. In particular, if $ \gamma = 1 $, then the subsonic root in the $ i^\text{th} $ pipe is given by

and hence

Note that $ R_i $ appears as a parameter in (27). Similar to the discussion in the previous section the conditions are stated for the traces of the equilibrium variables at $ x_o. $ The dependence on $ x_o $ is omitted for readability.

or more generally

where $ H $ could be any any quantity such as momentum flux or Bernoulli invariant as discussed in Section 2. Similarly, the coupling condition for a compressor in terms of $ K $ and $ L $ reads

For subsonic states, the coupling conditions (29), (30), (32) are equivalent to the coupling conditions (15), (18), (20).

For steady states, $ K_i $ and $ L_i $ are constant within each pipe, and the coupling conditions are fulfilled at the junction. Evaluating the coupling conditions in terms of the equilibrium variables $ V = (K, L) $ and the parameter $ R $ is an essential ingredient of the well-balancing. Therefore, we rewrite the Lax-curves in terms of the equilibrium variables:

Similarly for the admissible boundary states on the pipes $ i \in I^+ $ with given value $ R_r $ and $ u_r = q_r/\rho_r, $ we obtain

The equilibrium variables satisfying the coupling condition (29) and (30) are given by

Note that all variables defined along the Lax-curves also depend on the old traces and the integrated source terms as parameters, e.g. $ \overline{V}_{i}(\sigma_i) = \overline{V}_{i}(\sigma_i; V_i^o, R_i^o) $. For given datum $ U_l, U_r $ we depict the parameterized wave curves for incoming and outgoing pipes in the phase space of $ K $ and $ L $, respectively, in Figure 2. From the figure we observe that in the subsonic region, the 1–Lax curve is monotonically decreasing and 2–Lax curves is monotonically increasing. The following theorem proves that the coupling conditions stated in the variables $ K $ and $ L $ locally have a unique solution.

Theorem 3.1. Consider a nodal point with $ |I^-| \geq 1 $ incoming and $ |I^+| \geq 1 $ outgoing adjacent pipes. Suppose that the initial data $ \widehat{U_i} = (\hat{\rho}_{i}, \hat{q}_{i}) $, $ i\in I^\pm $ are subsonic on each pipe and fulfill the coupling conditions (15) and (16). Let $ {\hat V}_i = (\widehat{K_{{i}}^{{}}}, \widehat{L_{{i}}^{{}}}), i\in I^\pm $ be the corresponding equilibrium variables, with integrated source terms $ \widehat{R}_i $.

Then there exists an open neighborhood $ \mathcal{V} \subset \mathbb{R}^{ 2M \times M } $ of $ (\widehat{V}, \widehat{R}): = (\widehat{V}_i, \widehat{R}_i)_{i \in I^\pm} $ such that for any old trace $ (V^o, R^o)\in \mathcal{V} $ there exists a unique new trace $ V^* $ such that $ (V^*, R^o) \in \mathcal{V} $ fulfill the coupling conditions (29) and (30). Moreover, $ V_i^* $ is connected to $ V_i^o $ by an incoming Lax curve along the respective pipe. The neighborhood can be chosen sufficiently small, such that the corresponding states are subsonic.

Proof. Denote by $ M = |I^-| + |I^+| $ the total number of connected pipes. For $ V: = (V_i)_{i \in I^\pm} : = (K_i, L_i)_{i \in I^\pm} $ the coupling conditions (29) and (30) are given by the function $ {\Psi}: \mathbb{R}^{ 2 M }\times \mathbb{R}^{ 2 M } \times \mathbb{R}^M \to \mathbb{R}^{M }. $

where $ H $ is a given coupling condition of the form (16). By assumption we have $ \Psi({\hat V}, (\widehat{V}, \widehat{R})) = 0. $ Now, we define

where

and the components of $ \overline{V}_i(\sigma_i) $ are given by equation (33) and equation (34), respectively. Further, $ \sigma = (\sigma_i)_{i \in I^\pm} $. Next, we compute the determinant of $ D_\sigma \overline{\Psi}(\sigma, (V^o, R^o)) $ at $ \sigma = 0. $ We have for $ i \in I^- $ and $ j \in I^+ $

and therefore

From equations (33) and (34), we obtain at $ \sigma_i = 0 $

Hence,

Thus for any monotonic coupling condition, we see that $ \text{det}(D_\sigma\overline{\Psi})\ne 0 $. The coupling conditions such as pressure, momentum flux, Bernoulli invariant mentioned in Section 2 are all monotonic at the point $ \sigma_i = 0 $.

For example, for the pressure coupling condition (using $ H(K_i, L_i, R_i) = P_i $ in (31)), we get

Hence, and therefore

Similarly for coupling condition given by continuity of momentum flux, $ H(K_i, L_i, R_i) = L_i-R_i $ and hence

Thus, by the implicit function theorem there exists an open neighborhood $ \mathcal{V} $ of $ \widehat{V} $ such that for all initial data $ V^o \in \mathcal{V} $ there exists $ \sigma^* $ such that $ V^*: = \bar V(\sigma^*) $ fulfills the coupling conditions, i.e.

Since the corresponding state of $ \widehat{V} $ in conservative variables is strictly subsonic we may assume, by possibility decreasing the size of $ \mathcal{V}, $ that also the conservative variables corresponding to $ (V^*, R^o) $ are subsonic.

Using the same technique, we can also treat the compressor, because the pressure changes monotonically along the Lax curves.

Corollary 1. Consider a nodal point with $ |I^-| \geq 1 $ incoming and $ |I^+| \geq 1 $ outgoing pipes or a compressor connected to 1 incoming and 1 outgoing pipe. Suppose we have an equilibrium state with subsonic initial data $ \widehat{U_i} = (\hat{\rho}_{i}, \hat{q}_{i}) $, $ i\in I^\pm $ in each pipe and fulfilling the pressure coupling conditions (15) and (18) for a junction or (20) for compressor. Let $ {\hat V}_i = (\widehat{K_{{i}}^{{}}}, \widehat{L_{{i}}^{{}}}), i\in I^\pm $ be the corresponding equilibrium variables, with integrated source terms $ \widehat{R}_i $. Then there exists an open neighborhood $ \mathcal{V} \subset \mathbb{R}^{ 2M \times M } $ of $ (\widehat{V}, \widehat{R}): = (\widehat{V}_i, \widehat{R}_i)_{i \in I^\pm} $in the subsonic regime, such that for any old trace $ (V^o, R^o)\in \mathcal{V} $ there exists a unique new trace $ V^* $ such that $ (V^*, R^o) \in \mathcal{V} $ fulfill the coupling conditions (29) and (30) for a junction or (32) for a compressor.

Proof. The corollary follows from the proof of Theorem 3.1. One can also check that the pressure for a general Gamma law is also monotonic, i.e.,

In case of a compressor the pressure of the gas is multiplied by the compression ratio $ C\!R>0 $ and hence it does not affect the monotonicity of the coupling condition.

Remark 1. Note that we have omitted the source term when computing the solution along the Lax curves. This is justified by the semi-discrete formulation of the finite volume scheme in the next section, which implies that we are evaluating the coupling condition at a set of measure zero in space-time. Since the source term is bounded, it does not contribute to the integral over the cells.

Remark 2. Another possibility is to reformulate system (22) as a system of three equations in $ (K_i, L_i, R_i) $ with the equation for $ R_i $ given by

Therefore, the corresponding hyperbolic field has a zero eigenvalue in an independent subspace. This yields a characteristic boundary at each adjacent pipe. Hence in phase space, at each pipe $ i $ any value $ \tilde{R}_i $ can be connected along a wave curve to $ R_i $. The central wave leads to a contact discontinuity of zero velocity at the nodal point. Hence, the trace of $ R_i $ at $ x = x_o $ is independent of $ \tilde{R}_i. $

4.   A well-balanced central-upwind scheme for nodal dynamics

We compute the evolution of the conservative variables using the second–order central upwind scheme [26,25]. The computational domain $ \Omega_i $ is discretized in cells $ [x_{j-\frac{1}2}, x_{j+\frac{1}2}] $ of size $ \Delta x $ and centered at $ x_j = \bar{x} + (j-\frac{1}2)\Delta x $ for $ j = 1, \ldots , N $. We choose $ \bar{x} $ such that $ x_N = x_o $ for $ i \in I^- $ and $ x_0 = x_o $ for $ i\in I^+. $ For simplicity the same number of cells $ N $ for all adjacent pipes will be used. The approximated cell averages at fixed time $ t $ are computed as

The evolution of conservative variables, density and momentum using central upwind scheme [25,26] reads

where $ \mathcal{V}_i^{j-1/2}, \mathcal{V}_i^{j+1/2} $ are the fluxes across the left and right interface of cell $ j $, respectively. At the junction, the flux is the new trace of the equilibrium variable,

The new traces $ V_i^* $ are constructed with the help of Theorem 3.1 based on the old traces $ V_i^{N, E}, i\in I^- $ and $ V_i^{1, W}, i\in I^+ $. The point values of $ K $ and $ L $ at the cell interfaces, i.e., $ K_{{i}}^{{j}, {E}}, K_{{i}}^{{j}, {W}}, L_{{i}}^{{j}, {E}}, L_{{i}}^{{j}, {W}}, $ are computed using piecewise linear reconstruction of $ K_{{i}}^{{j}} $ and $ L_{{i}}^{{j}} $ calculated using the cell averages $ (\rho_{{i}}^{{j}}, q_{{i}}^{{j}}) $ using equation (24). The values $ R_i $ are computed using a second–order quadrature rule applied to the integral starting for example at $ \tilde{x} = x_o $ with $ R_i = 0 $ at each pipe according to the following equations

The equilibrium variables $ W \in \{ K, L \} $ at the right and left boundaries of cell $ j $ can be calculated as,

with numerical derivatives

$ \theta\in[1, 2] $ and minmod limiter

For interior interfaces, we may use any conservative numerical flux functions whose numerical diffusion vanishes at equilibrium states. Here we choose the central upwind flux,

where $ a_{i, \pm}^{j+1/2} $ are the maximum and minimum eigenvalues of the Jacobian, i.e.,

and $ \alpha_i^{j+1/2} $ is the local diffusion computed as $ \alpha_i^{j+1/2} = \frac{a_{i, +}^{j+1/2}a_{i, -}^{j+1/2}}{a_{i, +}^{j+1/2}-a_{i, -}^{j+1/2}}. $ The conservative variables $ U_{{i}}^{{j+1}, {W}}, U_{{i}}^{{j}, {E}} $ can be computed from the corresponding equilibrium variables $ V_{{i}}^{{j+1}, {W}}, V_{{i}}^{{j}, {E}} $ and the integral term $ R_i^{j+1/2} $ using equation (27).

Remark 3. Note that in [11] an additional limiter was introduced to suppress the numerical viscosity in (44) at equilibrium and assure well-balancing:in particular, the numerical viscosity $ \alpha_i^{j+1/2} $ was multiplied with a factor

where $ \phi $ was given by $ \frac{|K_{j+1}-K_j|}{\Delta x}\frac{|\Omega|}{max\{K_j, K_{j+1}\}} $ for the mass equation and analogously with $ K $ replaced by $ L $ for the momentum equation.

The following theorem shows that well-balancing is already assured by the continuity of the integrated source terms $ R_i^{j+1/2} $ at equilibrium.

Theorem 4.1. The numerical scheme given by (38) and flux defined by (44) preserves the steady state across a node of $ M $ adjacent pipes and coupling conditions given by (29) and (31).

Proof. Consider steady state $ (\widehat{V}, \widehat{R}): = (\widehat{V}_i, \widehat{R}_i)_{i \in I^\pm} $. Then all numerical derivatives defined in (42) vanish at equilibrium.

Since the equilibrium variables $ V_{{i}}^{{j+1}, {W}} = V_{{i}}^{{j}, {E}} = \widehat{V}_i $ and the integral term $ R_i^{j+1/2} $ are constant across the cell interface, we see that the conservative variables at the cell interfaces are also continuous at equilibrium i.e. $ U_{{i}}^{{j+1}, {W}} = U_{{i}}^{{j}, {E}} $. Thus the numerical fluxes are given by $ \mathcal{V}_i^{j+1/2} = \widehat{V}_i $ for all $ j = 2, \ldots , N-1 $.

At the nodal point the flux variables satisfy the coupling conditions (29) and (30). Then, the boundary data $ K_{{i}}^{{*}} $ and $ L_{{i}}^{{*}} $ are obtained according to Theorem 3.1. Since the states are unique we obtain $ K_{{i}}^{{j}} = K_{{i}}^{{*}} = \widehat{K}_i\:\text{ and }\: L_{{i}}^{{j}} = L_{{i}}^{{*}} = \widehat{L}_i $ and hence the boundary fluxes for each pipe at the junction are $ \mathcal{V}_i^{N+1/2} = \mathcal{V}_i^{N-1/2} = (\widehat{K}_i, \widehat{L}_i)^T $ for incoming pipes $ i \in I^- $ and $ \mathcal{V}_j^{1/2} = \mathcal{V}_j^{3/2} = (\widehat{K}_j, \widehat{L}_j)^T $ for outgoing pipes $ j \in I^+. $ Hence, the scheme is well-balanced across the node.

Some remarks are in order. The algorithm uses the same time step for all adjacent pipes. This is not necessary but simplifies the computation of the coupling condition. Also, the algorithm is second–order in the pipe but it may reduce to first order at the coupling condition. The algorithm can be extended to second–order across the nodal point using techniques presented in [2]. However, note that the steady state is constant and therefore the scheme preserves the steady state to any order across the nodal point.

5.   Numerical tests

In this section, we test the well-balanced(WB) scheme with numerical examples for steady state and near steady state flows. The results of this WB method have been compared with a second order non well-balanced method(NWB). The NWB scheme is given by,

where $ \mathcal{F} $ is the HLL flux given by,

where the flux terms are as defined in (2) and $ \mathcal{S}_i^j $ is the source term given in (2) at the point $ U_i^j $. The coupling conditions (15), (18) are used to calculate the density and momentum at a node. The coupling conditions at the node are solved with Newton's method for both WB and NWB schemes with initial guess for Newton's iteration given by the solution in the pipe at the node, $ V_i^{N, E}/U_i^{N, E} \forall i\in I^- $ and $ V_i^{1, W}/U_i^{1, W} \forall i\in I^+ $ for WB/ NWB schemes respectively.

We run the tests at CFL number = 0.4 and minmod parameter, $ \theta = 1 $. All the pipes in the examples have been considered to be of same diameter and friction factor, $ \frac{f_g}{2D} = 1 $ and the speed of sound for the gas, $ a = 1 $. We test the scheme for several well-balanced flows at junctions and compressors, as well as perturbations of such steady states.

5.1.   Steady state at a node

In the first example, we study the WB scheme for steady state at a node with three types of pipe combinations–1 incoming and 1 outgoing pipes; 1 incoming and 2 outgoing pipe; and 2 incoming and 1 outgoing pipe. The initial conditions are selected in such a way that the node is at steady state with the equilibrium variables constant in each pipe and satisfying the coupling conditions at the node.

The initial condition for first case with 1 incoming and 1 outgoing pipe are $ K_1 = K_2 = 0.15 $ and $ p^* = 0.332 $ corresponding to $ L_1 = L_2 = 0.4 $. Similarly for the second case, of 1 incoming and 2 outgoing pipe, $ K_1 = 0.15, K_2 = K_3 = 0.075 $ and $ p^* = 0.332 $ or $ L_1 = 0.4, L_2 = L_3 = 0.3492 $; and for 2 incoming and 1 outgoing pipes, $ K_3 = 0.15, K_1 = K_2 = 0.075 $ and $ p^* = 0.332 $ or $ L_3 = 0.4, L_1 = L_2 = 0.3492 $.

The L-1 error for the three cases are given in Table 1,

As can be seen from the results in Table 1, the L-1 error $ ||K-\widehat K||\text{ and }||L-\widehat L|| $ is accurate up to machine precision using the WB scheme, whereas it is of the order of $ 10^{-7} $ to $ 10^{-8} $ with the NWB scheme. We can also note that the coupling conditions in terms of (K, L) converge quickly using Newton's method and do not affect the well-balancing property of the scheme at the node.

5.2.   Steady state with a compressor

In the second example, we study the well-balancing at steady state for compressor connecting two pipes with compression ratios $ C\!R = 1.5, 2, 2.5 $. The initial conditions are selected in a way that the compressor is at steady state for time, T = 0. The momentum in the two pipes is given by $ K_1 = K_2 = 0.15 $ and pressure is given by $ p_1^* = 0.332 $ and $ p_2^* = CRp_1^* $. The L1 errors using the WB and NWB scheme are given in Table 2.

Similar to the first example, we see that the L1 errors using the WB scheme are accurate up to machine precision. Also the coupling conditions for the compressor do not affect the well-balancing of the scheme.

5.3.   Perturbations to steady state for a node

From the first two examples, we can see that the WB scheme preserves steady state. We now compare the WB and NWB scheme for initial conditions given by perturbations to momentum at steady state. The initial conditions for the perturbed state are given by,

where $ \widehat{K}_i $ and $ \widehat{L}_i $ are constant steady state equilibrium variables in the two pipes and $ \eta_i $ is the magnitude of perturbation at the node.

At first, we consider a node connecting two pipes. The equilibrium variables for this case are given by, $ \widehat{K}_i = 0.15 $ and $ \widehat{L}_i = 0.4 $. At first we consider perturbation of $ \eta_i = 10^{-3} $ at the junction. The momentum at time T = 0.2 are as shown in Figure 3,

We can see from the results that both the WB and NWB schemes provide similar solutions for the perturbation of order $ 10^{-3} $ at the node. We now reduce this perturbation to $ \eta_i = 10^{-6} $.

From Figure 4 we can see that the NWB scheme develops oscillations for N = 100 when the perturbation is of order $ 10^{-6} $. The perturbation is resolved better for a finer grid with N = 500 per pipe. However in the case of WB method, the scheme is able to capture the perturbations well even for a coarser grid of N = 100 per pipe.

We now do a similar test for a node connected to 1 incoming and 2 outgoing pipes. The equilibrium states are given by, $ \widehat{K}_1 = 0.15, \widehat{K}_2 = \widehat{K}_3 = 0.075 $ and $ \widehat{L}_1 = 0.4, \widehat{L}_2 = \widehat{L}_3 = 0.3492 $. We also run this simulation for two perturbations of order $ 10^{-3} $ and $ 10^{-6} $ up to a time T = 0.2. Figure 5 shows the result for momentum with $ \eta_1^* = 10^{-3}, \eta_2^* = \eta_3^* = 0.5\text{x}10^{-3} $ and Figure 6 for $ \eta_1^* = 10^{-6}, \eta_2^* = \eta_3^* = 0.5\text{x}10^{-6} $ respectively.

We see from the results that even for the perturbations of order $ 10^{-3} $, the results from NWB scheme are unstable when there is a sharp increase in momentum. The results of NWB scheme are even more oscillatory when the perturbations are of order $ 10^{-6} $. Further even with a finer resolution, we can see a spike in the region where there is a jump in momentum. However, these issues are resolved with the WB scheme. The results of WB scheme with coarser grid are a bit more diffusive than the finer grid, but there are no instabilities arising in the solution.

5.4.   Shock at the node

So far we have tested the scheme for smooth near-equilibrium solutions. Now we consider a discontinuity at a junction of 1 incoming and 2 outgoing pipes, i.e. $ \rho_1 = 5, \rho_2 = 4, \rho_3 = 3 $ and $ q_1 = q_2 = q_3 = 1 $. This results in a rarefaction and two shocks moving into the pipes, so the solution is far from equilibrium. The boundary condition is given by momentum at the inlet, $ q_{in} = 1 $ and density at outlet, $ \rho_{2, out} = 4, \rho_{3, out} = 3 $. The solution at time $ T = 0.1, 0.25 $ with WB and NWB scheme are shown in Figure 7 and Figure 8 respectively.

From the above test, we can see that the scheme is able to capture shocks. The solution shows a 1-wave moving towards left into the first pipe from the junction and a 2-wave towards right into the second and third pipes. Similar solutions were obtained by Egger[19] using a conservative mixed finite element method with a stagnation enthalpy coupling condition. In Figures 7 and 8, we use the classical, non-well -balanced scheme (using $ U $-variables) as reference and observe that the solutions are nearly identical. We can also note that the solution satisfies the coupling conditions at the node i.e. the pressure is constant at the node and sum of momentum of second and third pipe is equal to the momentum of first pipe. Note that the Newton iteration failed for this example when using the $ \mathcal{H} $-limiter (46).

6.   Conclusion

In this paper we have extended the recent, equilibrium variable based approach to well-balancing, developed by Chertock, Herty and Özcan[11] for one-dimensional systems, to a network of gas pipelines with friction. In particular we studied intersections of pipes at a node and compressors within a pipeline network. We prove well-posedness and well-balancing of the new scheme. For compressors and for junctions of three pipes, numerical experiments demonstrate that equilibria are resolved up to machine accuracy. Most interestingly, near equilibrium flows are resolved robustly and accurately, even in cases where a standard non-balanced scheme fails. For flows away from equilibrium, including shocks emanating from a junction, the scheme is as good as a standard scheme using conservative instead of equilibrium variables.

Acknowledgments

This work has been supported by HE5386/13–15, BMBF ENets Project, and DFG Research Training Group 2326: Energy, Entropy, and Dissipative Dynamics, RWTH Aachen University.