An asymptotic-preserving scheme for isentropic flow in pipe networks

Michael T. Redle; Michael Herty; Michael T. Redle; Michael Herty

doi:10.3934/nhm.2025013

Networks and Heterogeneous Media

2025, Volume 20, Issue 1: 254-285. doi: 10.3934/nhm.2025013

Previous Article Next Article

Research article

An asymptotic-preserving scheme for isentropic flow in pipe networks

Michael T. Redle ^{1
,
,},
Michael Herty ²

1.
Chair of Applied and Computational Mathematics, RWTH Aachen University, 52062 Aachen, Germany
2.
Chair of Numerical Analysis, Institute for Applied Mathematics (IGPM), RWTH Aachen University, 52062 Aachen, Germany

Received: 15 November 2024 Revised: 04 March 2025 Accepted: 13 March 2025 Published: 25 March 2025

We considered the simulation of isentropic flow in pipelines and pipe networks. Standard operating conditions in pipe networks suggested an emphasis to simulate low Mach and high friction regimes—however, the system was stiff in these regimes and conventional explicit approximation techniques proved quite costly and often impractical. To combat these inefficiencies, we developed a novel asymptotic-preserving scheme that was uniformly consistent and stable for all Mach regimes. The proposed method for a single pipeline followed the flux splitting suggested in Haack et al., in which the flux was separated into stiff and non-stiff portions then discretized in time using an implicit-explicit approach. The non-stiff part was advanced in time by an explicit hyperbolic solver; we opted for the second-order central-upwind finite volume scheme. The stiff portion is advanced in time implicitly using an approach based on Rosenbrock-type Runge-Kutta methods, which ultimately reduced this implicit stage to a discretization of a linear elliptic equation. To extend to full pipe networks, the scheme on a single pipeline was paired with coupling conditions defined at pipe-to-pipe intersections to ensure a mathematically well-posed problem. We showed that the coupling conditions remained well-posed at the low Mach/high friction limit—which, when used to define the ghost cells of each pipeline, resulted in a method that was accurate across these intersections in all regimes. The proposed method was tested on several numerical examples and produced accurate, non-oscillatory results with run times independent of the Mach number.

Keywords:

Citation: Michael T. Redle, Michael Herty. An asymptotic-preserving scheme for isentropic flow in pipe networks[J]. Networks and Heterogeneous Media, 2025, 20(1): 254-285. doi: 10.3934/nhm.2025013

Related Papers:

[1]	Gunhild A. Reigstad . Numerical network models and entropy principles for isothermal junction flow. Networks and Heterogeneous Media, 2014, 9(1): 65-95. doi: 10.3934/nhm.2014.9.65
[2]	Yogiraj Mantri, Michael Herty, Sebastian Noelle . Well-balanced scheme for gas-flow in pipeline networks. Networks and Heterogeneous Media, 2019, 14(4): 659-676. doi: 10.3934/nhm.2019026
[3]	Martin Gugat, Falk M. Hante, Markus Hirsch-Dick, Günter Leugering . Stationary states in gas networks. Networks and Heterogeneous Media, 2015, 10(2): 295-320. doi: 10.3934/nhm.2015.10.295
[4]	Michael Herty . Modeling, simulation and optimization of gas networks with compressors. Networks and Heterogeneous Media, 2007, 2(1): 81-97. doi: 10.3934/nhm.2007.2.81
[5]	Raimund Bürger, Christophe Chalons, Rafael Ordoñez, Luis Miguel Villada . A multiclass Lighthill-Whitham-Richards traffic model with a discontinuous velocity function. Networks and Heterogeneous Media, 2021, 16(2): 187-219. doi: 10.3934/nhm.2021004
[6]	Luis Almeida, Federica Bubba, Benoît Perthame, Camille Pouchol . Energy and implicit discretization of the Fokker-Planck and Keller-Segel type equations. Networks and Heterogeneous Media, 2019, 14(1): 23-41. doi: 10.3934/nhm.2019002
[7]	Maya Briani, Emiliano Cristiani . An easy-to-use algorithm for simulating traffic flow on networks: Theoretical study. Networks and Heterogeneous Media, 2014, 9(3): 519-552. doi: 10.3934/nhm.2014.9.519
[8]	Tingting Ma, Yayun Fu, Yuehua He, Wenjie Yang . A linearly implicit energy-preserving exponential time differencing scheme for the fractional nonlinear Schrödinger equation. Networks and Heterogeneous Media, 2023, 18(3): 1105-1117. doi: 10.3934/nhm.2023048
[9]	Mapundi K. Banda, Michael Herty, Axel Klar . Gas flow in pipeline networks. Networks and Heterogeneous Media, 2006, 1(1): 41-56. doi: 10.3934/nhm.2006.1.41
[10]	Junjie Wang, Yaping Zhang, Liangliang Zhai . Structure-preserving scheme for one dimension and two dimension fractional KGS equations. Networks and Heterogeneous Media, 2023, 18(1): 463-493. doi: 10.3934/nhm.2023019

Abstract

1. Introduction

This paper focuses on the development of a novel numerical method for gas flow in pipelines and pipe networks. To describe the gas transport, we use the isothermal/isentropic Euler equations with a source to account for the friction along the pipe walls:

$\begin{equation} \begin{aligned} & \rho_t + {\left( {\rho u} \right)}_x = 0, \\ &{\left( {\rho u} \right)}_t + {\left( {\rho u^2 +p} \right)}_x = -\frac{ \kappa}{2D} \rho u {\left| {u} \right|}, \end{aligned} \end{equation}$

(1.1)

where $x$ is the (one-dimensional) spatial variable, $t$ denotes time, $\rho$ is the fluid density, $u$ denotes the fluid velocity, $p = \rho^\gamma$ is the pressure under the isentropic assumption, in which $\gamma$ is the ratio of specific heats, and $\kappa$ denotes the Fanning friction coefficient. In the one-dimensional approach, which is quite common when modeling gas flows in pipes due to the ratio of pipe length $L$ and cross section $D$ (see e.g., ^[1,2]), the unknowns $\rho$ and $u$ are averaged over the (assumed constant) cross section.

The standard operating conditions hold the gas flow at moderate velocities (see e.g., ), ^[3], where they note a reference Mach number around 0.001; and, thus, we must account for the low Mach number or high friction regimes of system (1.1). However, the limiting solution as this parameter goes to zero brings about a number of difficulties—the most notable of which is that the underlying system becomes very stiff. This makes designing efficient and accurate methods to approximate the solution quite challenging. For example, conventional explicit schemes applied to system (1.1) would greatly over-resolve the solution in time, as the wave speeds of the system are proportional to the inverse of the reference Mach number—further implying that the Courant–Friedrichs–Lewy (CFL) stability restriction is proportional to this small parameter; see, e.g., ^[4,5,6]. Hence, explicit schemes may prove quite computationally expensive, especially in the low Mach/high friction regimes. Moreover, in the context of pipe networks, separate pipe segments may be in different regimes. Thus, to avoid (ⅰ) the high computational cost of explicit schemes in low Mach regimes; and (ⅱ) the potentially-expensive nonlinear solvers of implicit schemes in segments not near the asymptotic state, one may require different methods dependent on the location in the network.

One alternative to avoid this over-resolution issue (and, thus, the potential need of separate simulation codes dependent on location) is asymptotic-preserving (AP) schemes. By the definition in ^[7,8], a scheme is AP if the discretization of the continuous system remains consistent and stable regardless of the value taken for the underlying singular parameter (denoted by ${\varepsilon}$ in this paper). Furthermore, AP schemes allow for the space and time discretization to be independent of ${\varepsilon}$ ; i.e., the CFL condition for stability does not depend on the small parameter in the system.

Due to the potential speedup when discretizing with AP methods, they have attracted a lot of attention in the numerical and engineering community. AP schemes have been extensively studied for the kinetic equations; see, e.g., ^{[9,10,11,12,13,14,15]} and references therein. More recently, AP methods for the kinetic equations have been further extended to use AP Monte Carlo methods to reduce numerical diffusion effects ^[16,17], and to use AP neural networks to solve the underlying partial differential equation (PDE) system ^[18,19]. There has also been extensive development of AP methods in the low Mach regime of the isentropic Euler system ^{[20,21,22,23,24,25,26,27]}, compressible Euler equations ^{[28,29,30,31,32,33]}, and Navier-Stokes systems ^{[34,35,36,37]}, as well as in the low Froude regime of the standard shallow water equations ^{[38,39,40,41,42,43]} and variants with additional sources ^{[44,45,46,47]}.

The number of studies significantly shrinks, though, when it comes to the consideration of nonlinear friction terms that remain in the limiting solution. When this nonlinear friction term, such as that in system (1.1), remains in the low Mach/Froude regime, it adds the difficulty that this source must be discretized implicitly in some manner. If discretized naively, this requires a nonlinear solve, which would preferably be avoided. To our knowledge, only the studies in ^[44,48] consider some nonlinear friction term, in which only ^[44] avoids a nonlinear solve via their proposed semi-implicit hybrid finite volume/finite element method.

Furthermore, the number of studies of AP schemes in the extension to pipe networks is also quite limited. To our knowledge, there are only two such developments. In ^[49], they present an AP hybrid discontinuous Galerkin method on networks for the linear convection-diffusion equation. The work in ^[48] proposes a finite element AP method applied to the barotropic Euler equations with a nonlinear friction term on pipe networks. However, the aforementioned method does opt for an implicit time discretization reliant on a nonlinear solver.

In this paper, we propose an AP scheme for the isentropic Euler equations with a nonlinear friction source on pipe networks. The method on a single pipeline uses the hyperbolic flux splitting proposed in ^[26,47] to split the stiff and non-stiff parts of the system, which are then treated through an implicit-explicit approach. The non-stiff portion of the system is advanced in time explicitly and discretized in space using the second-order central-upwind (CU) finite volume scheme developed in ^[50,51]. The stiff part of the system is advanced in time implicitly using an approach related to Rosenbrock-type Runge-Kutta methods seen in, e.g., ^[52,53]. This, in combination with the selected flux splitting of ^[26,47], ultimately reduces this implicit stage to an elliptic equation that can then be solved linearly after discretizing in space via standard central difference. The proposed scheme for a single pipeline is extended to pipe networks by defining coupling conditions, such as those in ^[54], at pipe-to-pipe intersections to ensure a mathematically well-posed problem. We show that the coupling conditions remain well-posed in the low Mach/high friction regimes and use them to define the ghost cell values of each pipeline. The resulting first-order in time, second-order in space method is tested on several numerical examples of pipe networks and produces accurate and non-oscillatory results, in addition to running significantly faster than the analogous fully explicit scheme in the low Mach/high friction regimes.

This paper is organized as follows. In §2, we discuss the asymptotics of the isentropic Euler equations (1.1) on pipe networks and associated numerical difficulties. We then develop an AP scheme for the underlying system in §3. We present the performance of the proposed scheme on three examples in §4, and make some concluding remarks in §5.

2. Dimensional analysis

To look at the asymptotics of system (1.1), we non-dimensionalize the system by introducing characteristic time $t_0$ , characteristic density $\rho_0$ , characteristic velocity $w_0$ , and characteristic pressure $p_0$ , along with using the pipe length $L$ as the characteristic length. Therefore, the dimensionless quantities for system (1.1) read

$\begin{equation*} \widehat{x} = \frac{x}{L}, \qquad \widehat{t} = \frac{t}{t_0}, \qquad \widehat{\rho} = \frac{\rho}{\rho_0}, \qquad \widehat{u} = \frac{u}{w_0}, \qquad \widehat{p} = \frac{p}{p_0}. \end{equation*}$

Taking $w_0 = L/t_0$ , substituting these quantities into system (1.1), and dropping the hat notation for the sake of simplicity, we obtain the dimensionless isentropic Euler equations with a friction source term:

$\begin{equation*} \begin{aligned} & \rho_t + {\left( {\rho u} \right)}_x = 0, \\ &{\left( {\rho u} \right)}_t + {\left( {\rho u^2 + \frac{1}{ \text{Ma}^2}p} \right)}_x = -\frac{\kappa}{2\delta} \rho u {\left| {u} \right|}, \end{aligned} \end{equation*}$

where

$\text{Ma} = w_0\sqrt{\frac{\rho_0}{p_0}}, \qquad \delta = \frac{D}{L},$

are the reference Mach number and the ratio of the cross section to the pipe length, respectively. We then choose to take the reference Mach number $\text{Ma} = {\varepsilon}$ , and follow the suggestions of ^[3,48,55] in taking $\delta = {\varepsilon}^2/C_\delta$ , where $C_\delta$ is a constant, resulting in the parameterized system

$\begin{equation} \begin{aligned} & \rho_t + {\left( {\rho u} \right)}_x = 0, \\ &{\left( {\rho u} \right)}_t + {\left( {\rho u^2 + \frac{1}{{\varepsilon}^2}p} \right)}_x = -\frac{C_\delta \kappa}{2{\varepsilon}^2} \rho u {\left| {u} \right|}, \end{aligned} \end{equation}$

(2.1)

which can otherwise be written in the following vector form:

$\begin{equation} {\mathit{\boldsymbol{U}}}_t + {\mathit{\boldsymbol{F}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)}_x = {\mathit{\boldsymbol{S}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)}, \quad {\mathit{\boldsymbol{F}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)} = \begin{pmatrix} \rho u \\[4pt] \rho u^2 + p/{\varepsilon}^2 \end{pmatrix}, \quad {\mathit{\boldsymbol{S}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)} = \begin{pmatrix} 0 \\[4pt] -\frac{C_\delta \kappa}{2{\varepsilon}^2} \rho u {\left| {u} \right|} \end{pmatrix}, \end{equation}$

(2.2)

where ${\mathit{\boldsymbol{U}}}: = (\rho, \rho u)^\top$ , ${\mathit{\boldsymbol{F}}}({\mathit{\boldsymbol{U}}})$ denotes the flux, and ${\mathit{\boldsymbol{S}}}({\mathit{\boldsymbol{U}}})$ is the source due to friction along the pipe walls.

Remark 2.1. Note that in the limit ${\varepsilon} \rightarrow 0$ , it is clear that system (2.1) approaches a state where the spatial derivative of pressure balances the friction source due to the pipe walls. This follows a classical limit commonly seen in the modeling of pipelines and pipe networks; see, e.g., ^[3].

2.1. Continuous extension to pipe networks

To simulate an entire pipe network, we must of course consider appropriate boundary conditions at each pipe entrance and exit—the most complicated of which are at pipe-to-pipe intersections. These are obtained by defining some coupling conditions at the pipe-to-pipe intersections or junctions; see, e.g., ^[54,56]. In this section, we briefly describe some suitable condition options that may be prescribed at pipe junctions. The following coupling conditions presented are commonly used examples for the mathematical framework of pipeline networks; see, for example, ^[57] and references therein.

Let us denote the entrance and exit of pipe $k = 1, \ldots, K$ as $x_{\rm i}^{(k)}$ and $x_{\rm f}^{(k)}$ , respectively, and the variables in pipe $k$ as ${\mathit{\boldsymbol{U}}}^{(k)} = (\rho^{(k)}, (\rho u)^{(k)})^\top$ . Consider a single junction in which pipes with indices $1, \ldots, m$ denote the ingoing pipelines and those with indices $m+1, \ldots, K$ are the outgoing pipelines, as depicted in Figure 1. Then, at the junction, one must have the coupling condition for:

Figure 1. Illustration of a pipe network junction with

$m$ ingoing pipelines and

$K-m$ outgoing pipelines.

DownLoad: Full-Size Img PowerPoint

Conservation of momentum:

$\begin{equation} \sum\limits_{k = 1}^m (\rho u)^{(k)}(x_{\rm f}^{(k)}, t) = \sum\limits_{\ell = m+1}^K (\rho u)^{(\ell)}(x_{\rm i}^{(\ell)}, t), \end{equation}$

(2.3)

paired with one of the following:

$({\text{a}})$ Equal pressures (see, e.g., ^[54]):

$\begin{equation} p{\left( {\rho^{(k)}(x_{\rm f}^{(k)}, t)} \right)} = p{\left( {\rho^{(\ell)}(x_{\rm i}^{(\ell)}, t)} \right)} \qquad \forall k = 1, \dots, m, \qquad \ell = m+1, \dots, K; \end{equation}$

(2.4)

$({\text{b}})$ Equal momentums (see, e.g., ^[58,59]):

$\begin{equation} \begin{aligned} &{\left( {\rho^{(k)}(u^{(k)})^2} \right)}(x_{\rm f}^{(k)}, t) + \frac{1}{{\varepsilon}^2}p{\left( {\rho^{(k)}(x_{\rm f}^{(k)}, t)} \right)} = {\left( {\rho^{(\ell)}(u^{(\ell)})^2} \right)}(x_{\rm i}^{(\ell)}, t) + \frac{1}{{\varepsilon}^2}p{\left( {\rho^{(\ell)}(x_{\rm i}^{(\ell)}, t)} \right)}\\[6pt] & \quad\quad \quad \quad\quad \quad \quad \quad\quad \quad \quad \quad \quad \quad \quad\quad \quad \quad \quad \quad \quad \quad \quad \forall k = 1, \dots, m, \qquad \ell = m+1, \dots, K; \end{aligned} \end{equation}$

(2.5)

$({\text{c}})$ Geometry or flow-dependent pressure loss (see, e.g., ^[1,60]):

$\begin{equation} \frac{1}{{\varepsilon}^2}p{\left( {\rho^{(k)}(x_{\rm f}^{(k)}, t)} \right)} = \frac{1}{{\varepsilon}^2}p{\left( {\rho^{(\ell)}(x_{\rm i}^{(\ell)}, t)} \right)} - h_{k, \ell}, \qquad \forall k = 1, \dots, m, \qquad \ell = m+1, \dots, K, \end{equation}$

(2.6)

where $h_{k, \ell}$ denotes the pressure loss at the junction.

Notice that the ${\varepsilon}$ -dependence remains in the coupling conditions, as it arose from the non-dimensionalization of system (1.1).

The isentropic Euler system (2.1) paired with (2.3) and one of systems (2.4) and (2.5) has been proven a well-posed problem under the condition that the initial data (ⅰ) is not in a vacuum state ( $\rho^{(k)}(x, t = 0) > 0\ \forall k$ ); (ⅱ) has a flow direction that does not change ( $u^{(k)}(x, t = 0) \geq 0\ \forall k$ ); and (ⅲ) is under subsonic conditions ( $u^{(k)}(x, t = 0) \leq c\ \forall k$ ) ^[54,59]. These coupling conditions are implemented into the boundary conditions of each pipe $k$ ; we present how this is done for the proposed scheme in §3.5.1, following the rich literature in this field, e.g., ^{[61,62,63,64,65,66,67,68,69]}.

2.2. The low mach/high friction limit

To investigate the asymptotic behavior of system (2.1) inside each pipeline, we take the asymptotic expansion of variables $\rho$ and $u$ . Note that for simplicity, we removed the superscript $(k)$ introduced in the previous section for now, as the coupling conditions are only introduced on the boundaries. Thus, the associated asymptotic expansions of our unknowns are

$\begin{equation*} \begin{aligned} \rho & = \rho^{(0)} + {\varepsilon}^2 \rho^{(2)} + \cdots, \\[6pt] u & = u^{(0)} + {\varepsilon}^2 u^{(2)} + \cdots. \end{aligned} \end{equation*}$

Consequentially, we can obtain the asymptotic expansion for pressure $p = \rho^\gamma$ using Taylor series:

$\begin{equation} p = {\left( {\rho^{(0)}} \right)}^\gamma + {\varepsilon}^2 \gamma {\left( {\rho^{(0)} } \right)}^{\gamma-1}\rho^{(2)}+ \cdots. \end{equation}$

(2.7)

Note that the ${\varepsilon}^1$ terms are skipped since there are no $\mathcal{O}({\varepsilon}^{-1})$ terms appearing in system (2.1). Substituting the expansions of $\rho, \ u, \ p$ into system (2.1) and collecting like powers of ${\varepsilon}$ , we obtain the asymptotic behavior of the isentropic Euler equations with a friction source:

$\begin{equation} \begin{aligned} \mathcal{O}({\varepsilon}^{-2}):\quad & {\left[ {{\left( {\rho^{(0)}} \right)}^\gamma} \right]}_x = -\frac{C_\delta \kappa}{2}\rho^{(0)} u^{(0)} {\left| {u^{(0)}} \right|}, \\ \mathcal{O}(1):\quad & \rho^{(0)}_t + {\left( {\rho^{(0)}u^{(0)}} \right)}_x = 0, \\ &{\left( {\rho^{(0)}u^{(0)}} \right)}_t + {\left[ {\rho^{(0)}{\left( {u^{(0)}} \right)}^2 + \gamma {\left( {\rho^{(0)}} \right)}^{\gamma -1} \rho^{(2)}} \right]}_x \\ & = -\frac{C_\delta \kappa}{2} {\left( {\rho^{(2)} u^{(0)} {\left| {u^{(0)}} \right|} + \rho^{(0)} u^{(2)} {\left| {u^{(0)}} \right|} + \rho^{(0)} u^{(0)} {\left| {u^{(2)}} \right|}} \right)}. \end{aligned} \end{equation}$

(2.8)

Here, note that the $\mathcal{O}({\varepsilon}^{-2})$ asymptotic is exactly that discussed in Remark 2.1.

Similarly, we consider the ${\varepsilon} \rightarrow 0$ case on the coupling conditions at the pipe intersections. It is clear that (ⅰ) the condition related to conservation of momentum will remain as it is ${\varepsilon}$ -independent; and (ⅱ) regardless of your selection of coupling conditions (2.4), (2.5), or (2.6), all simplify to requiring equal pressures at the junction at the low Mach/high friction limit. Note that the different limiting equations have already been presented, e.g., in ^[3]. The main aspect here is on their numerical treatment.

2.3. Numerical Difficulties in Low Mach/High Friction regimes

By computing the eigenvalues of the Jacobian $\partial {\mathit{\boldsymbol{F}}}/ \partial {\mathit{\boldsymbol{U}}}$ , one can find that the wave speeds of system (2.2) are

${\left\{ { u \pm \frac{1}{{\varepsilon}}\sqrt{p'(\rho)}} \right\}}.$

In turn, this implies that if one was to solve system (2.1) using some standard explicit method with a uniform mesh spatial discretization with cell size ${{\Delta x}}$ , the corresponding time-step restriction due to the CFL condition would be

$\begin{equation} {{\Delta t}}_{\text{ex}} \leq \nu \frac{{{\Delta x}}}{ {\max\limits_x}{\left\{ {{\left| {u} \right|} + \frac{1}{{\varepsilon}}\sqrt{p'(\rho)}} \right\}}} = {\mathcal{O}({\varepsilon} {{\Delta x}})}, \end{equation}$

(2.9)

where $0 < \nu \leq 1$ denotes the CFL number. On top of this, explicit schemes typically have numerical diffusion of ${\mathcal{O}{\left({({{\Delta x}})^p/{\varepsilon}} \right)}}$ ^[26], where $p$ is the order of the method, further implying that one would need to select ${{\Delta x}} = {\mathcal{O}({\varepsilon}^{1/p})}$ to combat excessive numerical diffusion in the results. Therefore, to obtain unsmeared results, the time-step restriction for explicit schemes would need to be ${{\Delta t}}_{\text{ex}} = \mathcal{O}({\varepsilon}^{1+1/p})$ . In other words, explicit schemes are quite inefficient in the low Mach and high friction regimes due to the significant computational cost.

An alternative to avoid the heavy time-step restriction seen in system (2.9) would instead be to advance in time using an implicit method. However, this also could prove quite costly, as the nonlinearity of system (2.1) implies a dependence on some nonlinear iterative solver of an $N\times N$ system of equations, where $N$ is the number of cells in the spatial discretization.

Thus, we want a scheme that removes this ${\varepsilon}$ -dependence in the time-step restriction, converges to the asymptotics in system (2.8), and maintains the correct coupling conditions at the ${\varepsilon} \rightarrow 0$ limit. The scheme proposed in the following section aims to meet these desired properties on the discrete level.

3. Asymptotic-preserving scheme

To form an AP scheme for the system (2.1) on a single pipeline, we follow the hyperbolic flux splitting idea from ^[26,47]. To do so, we separate the slow and fast dynamics into two parts, resulting in the corresponding split system:

$\begin{equation} \begin{aligned} & \rho_t + \alpha{\left( {\rho u} \right)}_x + {\left( {1-\alpha} \right)} {\left( {\rho u} \right)}_x = 0, \\ &{\left( {\rho u} \right)}_t + {\left( {\rho u^2 + \frac{p - a(t) \rho}{{\varepsilon}^2}} \right)}_x + \frac{a(t)}{{\varepsilon}^2}\rho_x = -\frac{C_\delta \kappa}{2{\varepsilon}^2} \rho u {\left| {u} \right|}. \end{aligned} \end{equation}$

(3.1)

Equivalently, the split form Eq (3.1) can be written into an updated vector form, which reads

$\begin{equation} {\mathit{\boldsymbol{U}}}_t + \widetilde{{\mathit{\boldsymbol{F}}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)}_x + \widehat{{\mathit{\boldsymbol{F}}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)}_x = {\mathit{\boldsymbol{S}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)}, \end{equation}$

(3.2)

where

$\begin{equation} \widetilde{{\mathit{\boldsymbol{F}}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)} = \begin{pmatrix} \alpha \rho u \\ {\rho u^2 + \frac{p - a(t) \rho}{{\varepsilon}^2}} \end{pmatrix}, \qquad \widehat{{\mathit{\boldsymbol{F}}}}{\left( {{\mathit{\boldsymbol{U}}}} \right)} = \begin{pmatrix} {\left( {1-\alpha} \right)}\rho u \\ {\frac{ a(t) \rho}{{\varepsilon}^2}} \end{pmatrix}, \end{equation}$

(3.3)

are the slow (non-stiff) and fast (stiff) dynamics parts of the fluxes, respectively, and ${\mathit{\boldsymbol{S}}}({\mathit{\boldsymbol{U}}})$ is defined in Eq (2.2).

To guarantee the non-stiff subsystem ${\mathit{\boldsymbol{U}}}_t + \widetilde{{\mathit{\boldsymbol{F}}}}{\left({{\mathit{\boldsymbol{U}}}} \right)}_x = {{\bf{0}}}$ is indeed non-stiff and hyperbolic, $\alpha$ and $a(t)$ must be chosen appropriately to remove the $1/{\varepsilon}^2$ dependence on the wave speeds. Computed by finding the eigenvalues of the Jacobian $\partial \widetilde{{\mathit{\boldsymbol{F}}}}/ \partial {\mathit{\boldsymbol{U}}}$ , the wave speeds of this subsystem are

$\begin{equation} {\left\{ {u \pm \sqrt{{\left( {1-\alpha} \right)}u^2 + \frac{\alpha{\left[ {p'(\rho)-a(t)} \right]}}{{\varepsilon}^2}}} \right\}}. \end{equation}$

(3.4)

Thus, to ensure the eigenvalues of the non-stiff subsystem are real and are ${\mathcal{O}} (1)$ instead of ${\mathcal{O}} ({\varepsilon}^{-2})$ , as in the original system (2.1), we take

$\begin{equation} \alpha = {\varepsilon}^b \qquad \text{and} \qquad a(t) = \min\limits_{x} p'(\rho), \end{equation}$

(3.5)

and $b \geq 2$ . In turn, the new wave speeds in system (3.4) now avoid the dissipation and time-step issues seen in original system (2.1), allowing the slow dynamics to be discretized using any appropriate hyperbolic solver. We describe the hyperbolic solver we use in §3.2.

Remark 3.1. Note that, in comparison to the work of ^[26,47], we have less freedom in the selection of this parameter $b$ arising in Eq (3.5). This is due to the fact that the asymptotic behavior shown in Eq (2.8) allows for a nonconstant $\rho^{(0)}$ when ${\varepsilon} \rightarrow 0$ , meaning ${\left[{p'(\rho)-a(t)} \right]}$ is ${\mathcal{O}}(1).$

3.1. Time discretization

To discretize the split systems (3.2) and (3.3) in a way that relaxes the time-step stability restriction, we use the implicit-explicit (IMEX) approach; that is, we approximate the non-stiff flux terms $\widetilde{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})_x$ explicitly, and use an implicit approximation for the stiff flux terms $\widehat{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})_x$ and the friction source ${\mathit{\boldsymbol{S}}}({\mathit{\boldsymbol{U}}}).$ Since the source term ${\mathit{\boldsymbol{S}}} ({\mathit{\boldsymbol{U}}})$ is nonlinear, we opt for a time discretization related to Rosenbrock-type Runge-Kutta methods, which will still allow the system to be solved without the need of nonlinear solvers; see, e.g., ^[52,53] and references therein. To this end, we use the following first-order IMEX time discretization:

$\begin{align} &\frac{\rho^{n+1} - \rho^n}{{{\Delta t}}} + \alpha {\left( {\rho u} \right)}^n_x + {\left( {1-\alpha} \right)} {\left( {\rho u} \right)}^{n+1}_x = 0, \end{align}$

(3.6)

$\begin{align} &\frac{{\left( {\rho u} \right)}^{n+1} - {\left( {\rho u} \right)}^n}{{{\Delta t}}} + {\left( {\rho u^2 + \frac{p -a^n \rho}{{\varepsilon}^2}} \right)}^n_x + \frac{a^n}{{\varepsilon}^2}\rho^{n+1}_x = -\frac{C_\delta \kappa}{2{\varepsilon}^2} {\left| {u^n} \right|}{\left( {\rho u} \right)}^{n+1}. \end{align}$

(3.7)

The influence of the Rosenbrock-type method appears in the discretization of the source of system (3.7), in which only the $\rho u$ term is evaluated at time $t^{n+1}$ . For simplicity and notation purposes, let us define

$\begin{equation*} {\mathit{\boldsymbol{R}}}^n : = (R^{\rho, n}, R^{\rho u, n})^\top = -\widetilde{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})_x. \end{equation*}$

One can then solve Eqs (3.6) and (3.7) for $\rho^{n+1}$ and $(\rho u)^{n+1}$ , respectively, to obtain

$\begin{align} \rho^{n+1} & = \rho^n + {{\Delta t}} R^{\rho, n} - {{\Delta t}} (1-\alpha){\left( {\rho u} \right)}^{n+1}_x, \end{align}$

(3.8)

$\begin{align} [6pt] {\left( {\rho u} \right)}^{n+1} & = \frac{1}{\Psi^n(x)}{\left[ {{\left( {\rho u} \right)}^n + {{\Delta t}} R^{\rho u, n} - \frac{a^n {{\Delta t}}}{{\varepsilon}^2}\rho^{n+1}_x} \right]}, \end{align}$

(3.9)

where

$\begin{equation} \Psi^n(x) : = 1 + {{\Delta t}} \cdot \frac{C_\delta \kappa}{2{\varepsilon}^2} {\left| {u^n} \right|} \end{equation}$

(3.10)

is always greater than 1. Following then that of ^[23], we then differentiate Eq (3.9) with respect to $x$ and substitute this into Eq (3.8) to obtain an elliptic equation for $\rho^{n+1}:$

$\begin{equation} \rho^{n+1} - \frac{({{\Delta t}})^2}{{\varepsilon}^2}a^n(1-\alpha){\left( {\frac{\rho^{n+1}_x}{\Psi^n}} \right)}_x = \rho^n + {{\Delta t}} R^{\rho, n} -{{\Delta t}} (1-\alpha){\left[ {\frac{{\left( {\rho u} \right)}^n + {{\Delta t}} R^{\rho u, n}}{\Psi^n}} \right]}_x. \end{equation}$

(3.11)

Assuming all values at time $t^n$ are known, this implies that the choice of the Rosenbrock-type method in Eq (3.7) results in just a linear equation for the unknown $\rho^{n+1}$ . Furthermore, the time discretization in Eq (3.7) in combination with an appropriate spatial discretization, such as that further detailed in §3.2 and §3.3, will result in Eq (3.11) being a linear system that can be solved for $\rho^{n+1}$ . This solution of $\rho^{n+1}$ is then substituted into Eq (3.9) to obtain $(\rho u)^{n+1}$ .

3.2. Spatial discretization of non-stiff flux terms

We use the CU finite volume scheme to discretize the non-stiff flux terms $\widetilde{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})$ . For each pipe in the network, we split the domain into $N$ finite volume cells $C_j = [x_{{j-\frac{1}{2}}}, x_{{j+\frac{1}{2}}}]$ , where ${{\Delta x}} = x_{{j+\frac{1}{2}}} - x_{{j-\frac{1}{2}}}$ . Assume that at time level $t = t^n$ , the solution is realized in terms of its cell averages

$\begin{equation*} \bar{{\mathit{\boldsymbol{U}}}}_j^n = \frac{1}{{{\Delta x}}}\int_{C_j} {\mathit{\boldsymbol{U}}}(x, t^n)\ dx. \end{equation*}$

Then, we approximate the contribution of the non-stiff flux terms ${\mathit{\boldsymbol{R}}} : = -\widetilde{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})_x$ in each cell using

$\begin{equation} {\boldsymbol{{{\mathcal{{R}}}}}}^n_j = -\frac{\widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{{j+\frac{1}{2}}} - \widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{{j-\frac{1}{2}}}}{{{\Delta x}}}, \end{equation}$

(3.12)

where $\widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{j \pm \frac{1}{2}}$ are computed using the CU fluxes; see, e.g., ^[50,51,70]:

$\begin{equation} \widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{{j+\frac{1}{2}}} = \frac{ s^+_{{j+\frac{1}{2}}} \widetilde{{\mathit{\boldsymbol{F}}}}{\left( {{\mathit{\boldsymbol{U}}}^-_{{j+\frac{1}{2}}}} \right)} - s^-_{{j+\frac{1}{2}}} \widetilde{{\mathit{\boldsymbol{F}}}}{\left( {{\mathit{\boldsymbol{U}}}^+_{{j+\frac{1}{2}}}} \right)}}{s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} +\frac{s^+_{{j+\frac{1}{2}}} s^-_{{j+\frac{1}{2}}}}{s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}}{\left( {{\mathit{\boldsymbol{U}}}^+_{{j+\frac{1}{2}}} - {\mathit{\boldsymbol{U}}}^-_{{j+\frac{1}{2}}}} \right)}. \end{equation}$

(3.13)

Here, ${\mathit{\boldsymbol{U}}}^\pm_{{j+\frac{1}{2}}}$ are one-sided point values of ${\mathit{\boldsymbol{U}}}$ at the cell interface, computed at time $t = t^n$ using the piecewise linear reconstruction

$\begin{equation*} \widetilde{\mathit{\boldsymbol{U}}}^n(x) : = \bar{{\mathit{\boldsymbol{U}}}}^n_j + {\left( {{\mathit{\boldsymbol{U}}}_x} \right)}^n_j{\left( {x - x_j} \right)}, \quad x \in C_j. \end{equation*}$

Therefore, the values of the reconstruction at the interface $x_{{j+\frac{1}{2}}}$ are

$\begin{equation} {\mathit{\boldsymbol{U}}}^-_{{j+\frac{1}{2}}} = \bar{{\mathit{\boldsymbol{U}}}}^n_j + \frac{{{\Delta x}}}{2}{\left( {{\mathit{\boldsymbol{U}}}_x} \right)}^n_j, \qquad {\mathit{\boldsymbol{U}}}^+_{{j+\frac{1}{2}}} = \bar{{\mathit{\boldsymbol{U}}}}^n_{j+1} - \frac{{{\Delta x}}}{2}{\left( {{\mathit{\boldsymbol{U}}}_x} \right)}^n_{j+1}, \end{equation}$

(3.14)

where the slopes ${\left({{\mathit{\boldsymbol{U}}}_x} \right)}^n_j$ are computed using a nonlinear limiter to avoid oscillations. In this paper, we use the generalized minmod limiter (see, e.g., ^[71,72,73]):

$\begin{equation*} {\left( {{\mathit{\boldsymbol{U}}}_x} \right)}^n_j = { \text{minmod}}{\left( { \theta \frac{\bar{{\mathit{\boldsymbol{U}}}}^{n}_{j+1}-\bar{{\mathit{\boldsymbol{U}}}}^{n}_{j}}{{{\Delta x}}}, \ \frac{\bar{{\mathit{\boldsymbol{U}}}}^{n}_{j+1}-\bar{{\mathit{\boldsymbol{U}}}}^{n}_{j-1}}{2{{\Delta x}}}, \ \theta \frac{\bar{{\mathit{\boldsymbol{U}}}}^{n}_{j}-\bar{{\mathit{\boldsymbol{U}}}}^{n}_{j-1}}{{{\Delta x}}} } \right)}, \end{equation*}$

where the minmod function is

$\begin{equation*} \mbox{minmod}(z_1, z_2, \ldots) = \left\{\begin{aligned} &\min(z_1, z_2, \cdots)&&\mbox{if}\; z_i > 0, \; \forall i, \\ &\max(z_1, z_2, \cdots)&&\mbox{if}\; z_i < 0, \; \forall i, \\ &0&&\mbox{otherwise}, \end{aligned}\right. \end{equation*}$

and is applied component-wise. Lastly, $s^\pm_{{j+\frac{1}{2}}}$ of Eq (3.13) denote the one-sided speeds of propagation, estimated at time $t = t^n$ using the largest and smallest eigenvalues of the Jacobian $\partial \widetilde{{\mathit{\boldsymbol{F}}}}/ \partial {\mathit{\boldsymbol{U}}}$ seen in Eq (3.4):

$\begin{equation} \begin{aligned} s^+_{{j+\frac{1}{2}}} = & \max\left\{u^-_{{j+\frac{1}{2}}} + \sqrt{{\left( {1-\alpha} \right)}{\left( {u^-_{{j+\frac{1}{2}}}} \right)}^2 + \frac{\alpha}{{\varepsilon}^2}{\left[ { p'{\left( {\rho^-_{{j+\frac{1}{2}}}} \right)}-a^n} \right]}}, \right.\\ & \quad \left.u^+_{{j+\frac{1}{2}}} + \sqrt{{\left( {1-\alpha} \right)}{\left( {u^+_{{j+\frac{1}{2}}}} \right)}^2 + \frac{\alpha}{{\varepsilon}^2}{\left[ { p'{\left( {\rho^+_{{j+\frac{1}{2}}}} \right)}-a^n} \right]}}, 0\right\}, \\[6pt] s^-_{{j+\frac{1}{2}}} = & \min\left\{u^-_{{j+\frac{1}{2}}} - \sqrt{{\left( {1-\alpha} \right)}{\left( {u^-_{{j+\frac{1}{2}}}} \right)}^2 + \frac{\alpha}{{\varepsilon}^2}{\left[ { p'{\left( {\rho^-_{{j+\frac{1}{2}}}} \right)}-a^n} \right]}}, \right.\\ & \quad \left.u^+_{{j+\frac{1}{2}}} - \sqrt{{\left( {1-\alpha} \right)}{\left( {u^+_{{j+\frac{1}{2}}}} \right)}^2 + \frac{\alpha}{{\varepsilon}^2}{\left[ { p'{\left( {\rho^+_{{j+\frac{1}{2}}}} \right)}-a^n} \right]}}, 0\right\}. \end{aligned} \end{equation}$

(3.15)

3.3. Fully discrete method

To obtain the fully discrete AP scheme for cell $C_j$ , we take the elliptic equation for $\rho^{n+1}$ in (3.11) and the equation for $(\rho u)^{n+1}$ in (3.9) and (ⅰ) compute non-stiff flux terms ${\mathit{\boldsymbol{R}}}^n = -\widetilde{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})_x$ using the CU numerical fluxes described in §3.2; (ⅱ) discretize the first spatial derivative terms using standard second-order central difference; and (ⅲ) use a standard finite difference approximation for second derivatives on the term ${\left({\rho_x^{n+1}/\Psi^n} \right)}_x$ . Thus, the fully discrete AP scheme for $\bar \rho_j^{n+1}$ reads

$\begin{equation} \begin{aligned} &\bar{\rho}^{n+1}_j - \frac{({{\Delta t}})^2}{({{\Delta x}})^2{\varepsilon}^2}a^n(1-\alpha){\left[ {\phi^n_{{j+\frac{1}{2}}} \bar{\rho}^{n+1}_{j+1} - {\left( {\phi^n_{{j+\frac{1}{2}}} + \phi^n_{{j-\frac{1}{2}}}} \right)}\bar{\rho}^{n+1}_j + \phi^n_{{j-\frac{1}{2}}} \bar{\rho}^{n+1}_{j-1}} \right]}\\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad = \bar{\rho}^n_j+ {{\Delta t}} {\mathcal{{R}}}^{\rho, n}_j - {{\Delta t}} (1-\alpha) D_x {\left[ {\frac{1}{\Psi^n_j}{\left( {{\left( {\bar{\rho u}} \right)}^n_j + {{\Delta t}} {\mathcal{{R}}}^{\rho u, n}_j} \right)} } \right]}, \end{aligned} \end{equation}$

(3.16)

where ${\boldsymbol{{{\mathcal{{R}}}}}}^n_j = ({\mathcal{R}}^{\rho, n}_j, {\mathcal{R}}^{\rho u, n}_j)^\top$ is defined in Eqs (3.12) and (3.13), $\Psi_j^n = \Psi^n(x_j)$ from Eq (3.10),

$D_x\xi_j = \frac{\xi_{j+1}-\xi_{j-1}}{2{{\Delta x}}}$

is the discrete central difference operator, and

$\phi^n_{{j+\frac{1}{2}}} = \frac{1}{2}{\left[ {\frac{1}{\Psi^n_j} + \frac{1}{\Psi^n_{j+1}}} \right]}$

is used for the evaluation of $1/\Psi_{{j+\frac{1}{2}}}^n$ to keep the second-order spatial approximation in the elliptic solve. Note that the definition of $\Psi_j^n$ in Eq (3.10) results in $0 < \phi^n_{{j+\frac{1}{2}}} \leq 1$ , in turn implying that the corresponding matrix to be formed on the righthand side of Eq (3.16) is non-singular (by the Gershgorin theorem).

After solving the linear system in (3.16) for ${\left\{ {\bar \rho_j^{n+1}} \right\}}$ , it can be directly substituted into the spatially discretized version of Eq (3.9) to obtain the fully discrete approximation for ${\left({\bar{\rho u}} \right)}^{n+1}_j$ :

$\begin{equation} {\left( {\bar{\rho u}} \right)}^{n+1}_j = \frac{1}{\Psi^n_j}{\left[ {{\left( {\bar{\rho u}} \right)}^n_j + {{\Delta t}} {\mathcal{{R}}}^{\rho u, n}_j - \frac{a^n {{\Delta t}}}{{\varepsilon}^2}D_x\bar\rho^{n+1}_j} \right]}. \end{equation}$

(3.17)

3.4. Numerical diffusion and stability

In this section, we describe how the proposed method addresses the numerical difficulties of explicit schemes previously discussed in §2.3.

The stability of the proposed AP scheme is ensured by Lemma 3.1 of ^[26], which states that if each piece (the implicit and explicit portions) of the split method is stable, then the full method is also stable. For the fast (stiff) dynamics, one can show that implicit discretization in time, seen in Eqs (3.6) and (3.7), is stable for all $\rho^n$ and $u^n$ . Thus, stability of the proposed AP scheme is controlled by the time discretization of the slow (non-stiff) dynamics. The CFL condition for the slow dynamics can be found using the eigenvalues from Eq (3.4), resulting in the following time-step restriction:

$\begin{equation} {{\Delta t}}_{\text{AP}} \leq \nu \frac{{{\Delta x}}}{ {\max\limits_x}{\left\{ {{\left| {u} \right|} + \sqrt{{\left( {1-\alpha} \right)}u^2 + \frac{\alpha}{{\varepsilon}^2}{\left[ {p'(\rho)-a(t)} \right]}}} \right\}}}, \end{equation}$

(3.18)

in which, due to the selection of $\alpha$ and $a(t)$ made in Eq (3.5), the denominator is independent of ${\varepsilon}$ .

To ensure ${{\Delta t}}_{\text{AP}}$ is completely ${\varepsilon}$ -independent, we must also confirm the numerical diffusion of the proposed scheme does not have a ${\varepsilon}^{-1}$ dependence or worse. To do this, we analyze the leading order of numerical diffusion of the CU spatial discretization described in §3.2, and substitute this into the full discretization in Eqs (3.16) and (3.17). We start by rewriting the CU flux in Eq (3.13) in the form

$\begin{equation*} \widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{{j+\frac{1}{2}}} = \frac{\widetilde{{\boldsymbol{{F}}}}(\bar{{\mathit{\boldsymbol{U}}}}_j^n)+\widetilde{{\boldsymbol{{F}}}}(\bar{{\mathit{\boldsymbol{U}}}}_{j+1}^n)}{2} + {\boldsymbol{{{\mathcal{{D}}}}}}^n_{{j+\frac{1}{2}}}, \end{equation*}$

where

$\begin{equation} \begin{aligned} {\boldsymbol{{{\mathcal{{D}}}}}}^n_{{j+\frac{1}{2}}} = & \frac{s^+_{{j+\frac{1}{2}}}}{2{\left( {s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} \right)}}{\left[ {2 \widetilde{{\boldsymbol{{F}}}}({\mathit{\boldsymbol{U}}}^-_{{j+\frac{1}{2}}})-\widetilde{{\boldsymbol{{F}}}}(\bar{{\mathit{\boldsymbol{U}}}}_j^n)-\widetilde{{\boldsymbol{{F}}}}(\bar{{\mathit{\boldsymbol{U}}}}_{j+1}^n)} \right]}\\[4pt] &-\frac{s^-_{{j+\frac{1}{2}}}}{2{\left( {s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} \right)}}{\left[ {2 \widetilde{{\boldsymbol{{F}}}}({\mathit{\boldsymbol{U}}}^+_{{j+\frac{1}{2}}})-\widetilde{{\boldsymbol{{F}}}}(\bar{{\mathit{\boldsymbol{U}}}}_j^n)-\widetilde{{\boldsymbol{{F}}}}(\bar{{\mathit{\boldsymbol{U}}}}_{j+1}^n)} \right]}+\frac{s^+_{{j+\frac{1}{2}}} s^-_{{j+\frac{1}{2}}}}{s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}}{\left( {{\mathit{\boldsymbol{U}}}^+_{{j+\frac{1}{2}}} - {\mathit{\boldsymbol{U}}}^-_{{j+\frac{1}{2}}}} \right)}, \end{aligned} \end{equation}$

(3.19)

is the numerical diffusion term of the CU discretization, with ${\boldsymbol{{{\mathcal{{D}}}}}}_{{j+\frac{1}{2}}}^n : = ({\mathcal{D}}^{\rho, n}_{{j+\frac{1}{2}}}, {\mathcal{D}}^{\rho u, n}_{{j+\frac{1}{2}}})^\top$ to denote the different components. Here, $s^\pm_{{j+\frac{1}{2}}}$ are defined in Eq (3.15) and ${\mathit{\boldsymbol{U}}}^\pm_{{j+\frac{1}{2}}}$ are defined in Eq (3.14). To compute the leading order of the numerical diffusion in Eq (3.19), we introduce the following asymptotic expansions:

$\begin{align*} \bar \rho^n_j & = \rho^{(0), n}_j + {\varepsilon}^2 \rho^{(1), n}_j + \cdots, \\[6pt] \rho^\pm_{{j+\frac{1}{2}}} & = \rho^{(0), \pm}_{{j+\frac{1}{2}}} + {\varepsilon}^2 \rho^{(1), \pm}_{{j+\frac{1}{2}}} + \cdots. \end{align*}$

The asymptotic expansion of $u_j^n$ and $u^\pm_{{j+\frac{1}{2}}}$ follows analogously, and the discrete asymptotic expansion of the pressure term follows the same form as that in Eq (2.7). In addition, we expand $a^n$ in a way that follows its definition in Eq (3.5):

$a^n = \min\limits_x \gamma {\left( {\rho_j^{(0), n}} \right)}^{\gamma -1} + {\varepsilon}^2 a^{(2), n} + \cdots,$

where $a^{(2), n}$ is a constant in space at time $t^n$ . Using these expansions, we can obtain the leading order of the numerical diffusion (3.19) for $\rho$ :

$\begin{equation} \begin{aligned} {\mathcal{{D}}}^{\rho, n}_{{j+\frac{1}{2}}} = & \frac{\alpha s^+_{{j+\frac{1}{2}}}}{2{\left( {s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} \right)}}{\left[ {2 \rho^{(0), -}_{{j+\frac{1}{2}}} u^{(0), -}_{{j+\frac{1}{2}}}-\rho^{(0), n}_j u^{(0), n}_j-\rho^{(0), n}_{j+1} u^{(0), n}_{j+1}} \right]}\\[4pt] &-\frac{\alpha s^-_{{j+\frac{1}{2}}}}{2{\left( {s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} \right)}}{\left[ {2 \rho^{(0), +}_{{j+\frac{1}{2}}} u^{(0), +}_{{j+\frac{1}{2}}}-\rho^{(0), n}_j u^{(0), n}_j-\rho^{(0), n}_{j+1} u^{(0), n}_{j+1}} \right]}\\ &+\frac{s^+_{{j+\frac{1}{2}}} s^-_{{j+\frac{1}{2}}}}{s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}}{\left( {\rho^{(0), +}_{{j+\frac{1}{2}}} - \rho^{(0), -}_{{j+\frac{1}{2}}}} \right)}. \end{aligned} \end{equation}$

(3.20)

Since $\alpha = {\varepsilon}^b, \ b \geq 2$ from Eq (3.5), this is ${\mathcal{O}}(1)$ by the last term of Eq (3.20), as $\rho$ is generally not constant in space. Similarly, one can compute the leading order of numerical diffusion (3.19) for $\rho u$ , which comes from the ${\mathcal{O}}({\varepsilon}^{-2})$ portion of the flux $\widetilde{{\mathit{\boldsymbol{F}}}}({\mathit{\boldsymbol{U}}})$ seen in Eq (3.2). However, since $a^n$ and its expansion are related to $p'(\rho)$ and $p(\rho) = \rho^\gamma$ , it is clear that the expansions to compute the diffusion will not result in any cancellation of the ${\varepsilon}^{-2}$ term. Thus, ${\mathcal{{D}}}^{\rho u, n}_{{j+\frac{1}{2}}} = {\mathcal{O}}({\varepsilon}^{-2})$ .

Finally, substituting the CU fluxes into Eq (3.12), we can obtain the following approximations to ${\boldsymbol{{{\mathcal{R}}}}}_j^n$ :

$\begin{align} {\mathcal{R_j}}^{\rho, n} & = - D_x\widetilde{F}^\rho({\mathit{\boldsymbol{U}}}_j^n) + {\mathcal{O}}(({{\Delta x}})^2), \end{align}$

(3.21)

$\begin{align} [6pt] {\mathcal{R_j}}^{\rho u, n} & = - D_x\widetilde{F}^{\rho u}({\mathit{\boldsymbol{U}}}_j^n) + {\mathcal{O}}{\left( {\frac{({{\Delta x}})^2}{{\varepsilon}^2}} \right)}, \end{align}$

(3.22)

where $\widetilde{{\mathit{\boldsymbol{F}}}} : = (\widetilde{F}^\rho, \widetilde{F}^{\rho u})^\top$ and $D_x$ again represents the central difference operator. Here, the ${\mathcal{O}}(({{\Delta x}})^2)$ term in Eq (3.21) and the ${\mathcal{O}}({\varepsilon}^{-2}({{\Delta x}})^2)$ term in Eq (3.22) come from the leading orders from the diffusion ${\mathcal{{D}}}^{\rho, n}_{{j+\frac{1}{2}}}$ and ${\mathcal{{D}}}^{\rho u, n}_{{j+\frac{1}{2}}}$ , respectively, and the $({{\Delta x}})^2$ piece arises since these diffusion terms are introduced via the second-order CU discretization.

Normally, this ${\mathcal{O}}({\varepsilon}^{-2}({{\Delta x}})^2)$ diffusion term is concerning. However, if we substitute Eq (3.22) into the fully discrete AP scheme seen in Eqs (3.16) and (3.17), we see that this ${\mathcal{O}}({\varepsilon}^{-2}({{\Delta x}})^2)$ diffusion term is always paired with a division by $\Psi_j^n$ . Then, since $\Psi_j^n$ defined in Eq (3.10) is ${\mathcal{O}}({\varepsilon}^{-2})$ , the dependence of ${\varepsilon}$ in the leading order diffusion term cancels. This ${\varepsilon}$ -independence within the numerical diffusion, in combination with the denominator of Eq (3.18) being independent of ${\varepsilon}$ , implies that the time-step restriction of the proposed AP scheme is indeed independent of ${\varepsilon}$ . Furthermore, it is sufficient to take ${{\Delta t}}_{\rm AP} = {\mathcal{O}}({{\Delta x}})$ (as opposed to the ${\mathcal{O}}({\varepsilon}{{\Delta x}})$ restriction on explicit schemes) to enforce stability.

3.5. Boundary conditions on pipe networks

3.5.1. Coupling conditions at pipe intersections

The algorithm described in §3.1–§3.3 is applied to every pipe within the network of interest. To extend to a full network, we must implement coupling conditions (2.3), with one choice of Eq (2.4), (2.5), or (2.6), on the boundaries of each pipe that meets at the intersections. To obtain the boundary values corresponding to the junction, we solve the so-called half-Riemann problem: Assume a set of coupling conditions from §2.1 and given constant initial data $({\mathit{\boldsymbol{U}}}^0)^{(k)}$ for each pipe $k$ . We then seek the values ${\mathit{\boldsymbol{U}}}^*$ satisfying the coupling conditions such that the half-Riemann problem ^[74]

$\begin{equation} \begin{aligned} & {\mathit{\boldsymbol{U}}}^{(k)}_t + {\mathit{\boldsymbol{F}}}({\mathit{\boldsymbol{U}}}^{(k)})_x = {{\bf{0}}}, \\[6pt] &{\mathit{\boldsymbol{U}}}(x, 0) = \left\{ \begin{aligned} &{\mathit{\boldsymbol{U}}}^*, && \text{if } x \leq 0, \\ & ({\mathit{\boldsymbol{U}}}^0)^{(k)}, && \text{if } x > 0, \end{aligned}\right. \end{aligned} \end{equation}$

(3.23)

admits a self-similar solution in which all generated waves have nonnegative speeds. Here, ${\mathit{\boldsymbol{U}}}$ and ${\mathit{\boldsymbol{F}}}({\mathit{\boldsymbol{U}}})$ are defined in Eq (2.2). Then, if considering coupling conditions (2.3) and (2.4), this amounts to solving the following nonlinear system for $(\rho^*)^{(k)}$ :

$\begin{equation} \begin{aligned} &\sum\limits_{k = 1}^K L_2^-{\left( {(\rho^*)^{(k)}; (\rho^n)^{(k)}, L_1^+{\left( {(\rho^n)^{(k)};({\mathit{\boldsymbol{U}}}^0)^{(k)}} \right)} } \right)} = 0, \\[6pt] &p{\left( {(\rho^*)^{(k)}} \right)} = p{\left( {(\rho^*)^{(\ell)}} \right)} , \qquad k \neq \ell, \end{aligned} \end{equation}$

(3.24)

where the known states are $({\mathit{\boldsymbol{U}}}^0)^{(k)}$ and $({\mathit{\boldsymbol{U}}}^n)^{(k)}$ , and the functions $L_{1, 2}^\pm$ denote the forward and reversed 1-Lax curve and 2-Lax curve for the isothermal Euler equations in Eq (2.1) with $p = \rho^\gamma$ (see, similarly, ^[59]):

$\begin{equation} \begin{aligned} L_1^+{\left( {\rho; \widehat{\rho}, \widehat{q}} \right)} & = \left\{ \begin{aligned} & \frac{\rho}{\widehat{\rho}} \ \widehat{q} - \frac{2}{\gamma - 1} \rho \cdot \frac{1}{{\varepsilon}}{\left( {\sqrt{p'(\rho)}-\sqrt{p'(\widehat{\rho})}} \right)}, && \text{if } \rho < \widehat{\rho}, \\[6pt] &\frac{\rho}{\widehat{\rho}} \ \widehat{q} - \frac{1}{{\varepsilon}}\sqrt{\frac{\rho}{\widehat{\rho}}{\left( {\rho - \widehat{\rho}} \right)}{\left[ {p(\rho)-p(\widehat{\rho})} \right]} }, && \text{if } \rho > \widehat{\rho}, \end{aligned}\ \right.\\[8pt] L_2^-{\left( {\rho; \widehat{\rho}, \widehat{q}} \right)} & = \left\{ \begin{aligned} & \frac{\rho}{\widehat{\rho}} \ \widehat{q} - \frac{2}{\gamma - 1} \rho \cdot \frac{1}{{\varepsilon}}{\left( {\sqrt{p'(\widehat{\rho})}-\sqrt{p'(\rho)}} \right)}, && \text{if } \rho < \widehat{\rho}, \\[4pt] &\frac{\rho}{\widehat{\rho}} \ \widehat{q} + \frac{1}{{\varepsilon}}\sqrt{\frac{\rho}{\widehat{\rho}}{\left( {\rho - \widehat{\rho}} \right)}{\left[ {p(\rho)-p(\widehat{\rho})} \right]} }, && \text{if } \rho > \widehat{\rho}. \end{aligned}\ \right. \end{aligned} \end{equation}$

(3.25)

The construction of system (3.24) would follow analogously if we instead used the coupling conditions in (2.5) or (2.6).

The known states $({\mathit{\boldsymbol{U}}}^n)^{(k)}$ seen in system (3.24) are selected by taking the cell value nearest the junction at time $t = t^n$ ; i.e., we take $({\mathit{\boldsymbol{U}}}^n)^{(k)} = ({\mathit{\boldsymbol{U}}}^n_{N})^{(k)}$ for ingoing pipelines and $({\mathit{\boldsymbol{U}}}^n)^{(k)} = ({\mathit{\boldsymbol{U}}}^n_{1})^{(k)}$ for outgoing pipelines. The nonlinear system (3.24) is solved via Newton's method with initial guess ${\mathit{\boldsymbol{U}}}^* = ({\mathit{\boldsymbol{U}}}^n)^{(k)}$ , in which the Jacobian within the Newton iteration is invertible if the initial guess is subsonic ^[59]. Note the $1/{\varepsilon}$ dependence in the 1-Lax curve and 2-Lax curve only arises due to considering the non-dimensionalized system (2.2) and does not destroy the nonlinear solve of system (3.24). This is further verified in Lemma 3.2. For the numerical experiments conducted in this paper, we observed convergence to a tolerance of $10^{-8}$ within $\sim$ 2–3 Newton iterations.

Once $(\rho^*)^{(k)}$ is computed, we can directly calculate $(\rho^* u^*)^{(k)}$ using

$(\rho^* u^*)^{(k)} = L_2^-{\left( {(\rho^*)^{(k)}; (\rho^n)^{(k)}, L_1^+{\left( {(\rho^n)^{(k)};({\mathit{\boldsymbol{U}}}^0)^{(k)}} \right)} } \right)}.$

The results of the nonlinear solve $(\rho^*)^{(k)}$ and $(\rho^* u^*)^{(k)}$ are then prescribed as the ghost cell values in the spatial discretization described in §3.2.

Lemma 3.2. At the limiting case ${\varepsilon} \rightarrow 0$ , the nonlinear system (3.24) with Lax curve definitions (3.25) remains a well-posed problem.

Proof. Since the second condition of system (3.24) is ${\varepsilon}$ -independent, we only need to focus on the first condition. For brevity, let us first rewrite the Lax-curve definitions in system (3.25) as

$L_1^+{\left( {\rho; \widehat{\rho}, \widehat{q}} \right)} = \left\{ \begin{aligned} & \frac{\rho}{\widehat{\rho}} \ \widehat{q} - \frac{1}{{\varepsilon}}g(\rho, \widehat{\rho}), && \text{if } \rho < \widehat{\rho}, \\[4pt] &\frac{\rho}{\widehat{\rho}} \ \widehat{q} - \frac{1}{{\varepsilon}}h(\rho, \widehat{\rho}), && \text{if } \rho > \widehat{\rho}, \end{aligned}\ \right. \quad \text{and}\quad L_2^-{\left( {\rho; \widehat{\rho}, \widehat{q}} \right)} = \left\{ \begin{aligned} & \frac{\rho}{\widehat{\rho}} \ \widehat{q} + \frac{1}{{\varepsilon}}g(\rho, \widehat{\rho}), && \text{if } \rho < \widehat{\rho}, \\[4pt] &\frac{\rho}{\widehat{\rho}} \ \widehat{q} + \frac{1}{{\varepsilon}}h(\rho, \widehat{\rho}), && \text{if } \rho > \widehat{\rho}, \end{aligned}\ \right.$

where $g(\rho, \widehat{\rho})$ and $h(\rho, \widehat{\rho})$ are taken such that these are equivalent to system (3.25). In this proof, we will only consider the case in which $\rho > \widehat{\rho}$ for both $L_1^+$ and $L_2^-$ , and the proofs for the other three pairings follow analogously. Following this assumption, Lax curves for this case read

$L_1^+{\left( {\rho; \widehat{\rho}, \widehat{q}} \right)} = \frac{\rho}{\widehat{\rho}} \ \widehat{q} - \frac{1}{{\varepsilon}}h(\rho, \widehat{\rho}) \qquad \text{and} \qquad L_2^-{\left( {\rho; \widehat{\rho}, \widehat{q}} \right)} = \frac{\rho}{\widehat{\rho}} \ \widehat{q} + \frac{1}{{\varepsilon}}h(\rho, \widehat{\rho}).$

Note that in the context of the nonlinear system (3.24), the case in which $\rho > \widehat{\rho}$ for both $L_1^+$ and $L_2^-$ is equivalent to having $(\rho^0)^{(k)} < (\rho^n)^{(k)}$ and $(\rho^n)^{(k)} < (\rho^*)^{(k)}$ . We can then use this specific case of the Lax curves and substitute directly into system (3.24) to obtain the nonlinear system:

$\begin{equation} \begin{aligned} &\sum\limits_{k = 1}^K \frac{(\rho^*)^{(k)}}{(\rho^n)^{(k)}} {\left[ {\frac{(\rho^n)^{(k)}}{(\rho^0)^{(k)}} (q^0)^{(k)} -\frac{1}{{\varepsilon}}h{\left( {(\rho^n)^{(k)}, (\rho^0)^{(k)}} \right)} } \right]} +\frac{1}{{\varepsilon}}h{\left( {(\rho^*)^{(k)}, (\rho^n)^{(k)}} \right)} = 0, \\[6pt] &p{\left( {(\rho^*)^{(k)}} \right)} = p{\left( {(\rho^*)^{(\ell)}} \right)} , \qquad k \neq \ell. \end{aligned} \end{equation}$

(3.26)

This system with ${\varepsilon} > 0$ has already been proven to be a well-posed system in ^[59]. We then multiply the first equation of system (3.26) by ${\varepsilon}$ and take the limit ${\varepsilon} \rightarrow 0$ to obtain an equivalent nonlinear system in which we wish to solve for $(\rho^*)^{(k)}$ , which reads

$\begin{equation} \begin{aligned} &\sum\limits_{k = 1}^K h{\left( {(\rho^*)^{(k)}, (\rho^n)^{(k)}} \right)} -\frac{(\rho^*)^{(k)}}{(\rho^n)^{(k)}} h{\left( {(\rho^n)^{(k)}, (\rho^0)^{(k)}} \right)} = 0, \\ &p{\left( {(\rho^*)^{(k)}} \right)} = p{\left( {(\rho^*)^{(\ell)}} \right)} , \qquad k \neq \ell. \end{aligned} \end{equation}$

(3.27)

This is equivalent to that of system (3.26) with $(q^0)^{(k)}$ being taken as zero. Therefore, this is also already proven as a well-posed problem by the work in ^[59]. Hence, the system (3.24) is well-posed at the limiting case ${\varepsilon} \rightarrow 0$ . □

Remark 3.3. As stated in §2.2, the other coupling conditions (2.5) and (2.6) reduce to the pressure balance coupling condition in (2.4) at the limiting case ${\varepsilon} \rightarrow 0$ . Thus, Lemma 3.2 is additionally valid for coupling conditions (2.5) and (2.6), as their corresponding nonlinear systems would also reduce to Eq (3.27) at the low Mach/high friction limit.

Remark 3.4. There exists also other approaches to discretize the coupling conditions. Here, we mention the possibility to do so fully implicit in time ^[75], using finite-element based approaches ^[76], finite difference discretizations ^[77,78], linearization approaches ^[79], or central scheme type discretizations ^[80]. As the previous lemma shows that the treatment using half-Riemann problems is, however, sufficient to guarantee the asymptotic-preserving property of the scheme, we therefore do not investigate those alternative approaches.

3.5.2. Additional boundary condition requirements

In addition to defining the ghost cell values for $\rho$ and $\rho u$ at the single point where the pipes meet, the algorithm described in §3.1–§3.3 requires defined values for (ⅰ) the reconstructed values on both sides of the cell interfaces $x_{\frac{1}{2}}^{(k)}$ and $x_{N+\frac{1}{2}}^{(k)}$ ; and (ⅱ) the central difference derivative of numerical fluxes near pipeline boundaries.

The requirement for reconstructed values at the domain boundaries $x_{\frac{1}{2}}^{(k)}$ and $x_{N+\frac{1}{2}}^{(k)}$ of each pipeline is standard and directly needed within the numerical flux evaluations; see Eq (3.13). These evaluations are trivial at the inlets and outlets of the pipe network, as they typically involve Dirichlet or Neumann boundary conditions. Again, the difficulty lies at the junctions of the network. Since we do not have enough information to form a piecewise-linear reconstruction (via Eq (3.14)) of the ghost cells sitting at the junctions, we instead opt for the first-order reconstruction; that is, at pipe intersections, we take

$\begin{equation*} ({\mathit{\boldsymbol{U}}}^\pm_{N+\frac{1}{2}})^{(k)} = {\mathit{\boldsymbol{U}}}^* \qquad \text{and} \qquad ({\mathit{\boldsymbol{U}}}^\pm_{\frac{1}{2}})^{(k)} = {\mathit{\boldsymbol{U}}}^*, \end{equation*}$

for ingoing and outgoing pipelines, respectively. Here, ${\mathit{\boldsymbol{U}}}^*$ denotes the solution to nonlinear system (3.24).

The other boundary condition issue lies in the central difference within the fully discrete method to calculate $\bar \rho^{n+1}_j$ . As seen in Eq (3.16), we require a central difference of ${\mathcal{R_j}}^{\rho u, n}$ , which by Eq (3.12) implies the need to evaluate the unknowns $\widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{-\frac{1}{2}}$ and $\widetilde{{\boldsymbol{{{\mathcal{{F}}}}}}}^n_{N+\frac{3}{2}}$ for each pipe. This would require even more ghost cell evaluations than the reconstructions on the boundaries. Thus, we instead assume a zero-extrapolation of ${\mathcal{R_j}}^{\rho u, n}$ for cells $C_0$ and $C_{N+1}$ and use the extrapolated values in the central difference seen in Eq (3.16); that is, we use

$D_x{\left[ {\frac{{\mathcal{R_1}}^{\rho u, n}}{\Psi_1^n}} \right]} = \frac{1}{{{\Delta x}}}{\left[ {\frac{{\mathcal{R_2}}^{\rho u, n}}{\Psi_2^n} - \frac{{\mathcal{R_1}}^{\rho u, n}}{\Psi_0^n}} \right]} \qquad \text{and} \qquad D_x{\left[ {\frac{{\mathcal{R_N}}^{\rho u, n}}{\Psi_N^n}} \right]} = \frac{1}{{{\Delta x}}}{\left[ {\frac{{\mathcal{R_N}}^{\rho u, n}}{\Psi_{N+1}^n} - \frac{{\mathcal{R}}_{N-1}^{\rho u, n}}{\Psi_{N-1}^n}} \right]}.$

Both this and the first-order reconstructions at pipe intersections imply a first-order method at the junctions and, thus, everywhere in the domain. In the near future, we intend to expand upon this, the first-order time discretization, and the piecewise-constant initial data within the half-Riemann problem (3.23) to enhance the proposed scheme to be fully second-order accurate.

3.6. The discrete low mach/high friction limit

3.6.1. The asymptotic-preserving property

In this section, we show that in the limit as ${\varepsilon} \rightarrow 0$ , the method converges to the asymptotic state in Eq (2.8). Consider the following asymptotic expansions for the discretization of $\rho$ , $u$ , and $p$ , which follow from §2.2:

$\begin{equation} \begin{aligned} \bar\rho_j^n & = \bar\rho^{(0), n}_j + {\varepsilon}^2 \bar\rho^{(2), n}_j + \cdots, \\ u_j^n & = u^{(0), n}_j + {\varepsilon}^2 u^{(2), n}_j + \cdots, \\ p(\bar\rho_j^n) & = {\left( {\rho^{(0), n}_j} \right)}^\gamma + {\varepsilon}^2 \gamma {\left( {\rho^{(0), n}_j } \right)}^{\gamma-1}\rho^{(2), n}_j+ \cdots, \end{aligned} \end{equation}$

(3.28)

and analogously at time $t^{n+1}$ . We want to show that as the proposed method advances in time, it keeps the asymptotic state up to the predicted order of accuracy. To this end, we wish to find a bound on the asymptotic expansion for the time-advanced solution

${\left[ {{\left( {\bar\rho_j^{(0), n+1}} \right)}^\gamma} \right]}_{x}+\frac{C_\delta \kappa}{2}\bar{\rho}^{(0), n+1}_j u^{(0), n+1}_j {{\left| {u_j^{(0), n+1}} \right|}}.$

Since the pressure to friction balance in the asymptotic state arises from the momentum equation in system (2.1), we start with a slightly manipulated form of the discrete approximation for ${\left({\bar{\rho u}} \right)}^{n+1}_j$ in Eq (3.17):

$\frac{{\varepsilon}^2}{{{\Delta t}}}{\left( {{\left( {\bar{\rho u}} \right)}^n_j + {{\Delta t}} {\mathcal{{R}}}^{\rho u, n}_j - \frac{a^n {{\Delta t}}}{{\varepsilon}^2}D_x\bar\rho^{n+1}_j-\Psi^n_j{\left( {\bar{\rho u}} \right)}^{n+1}_j } \right)} = 0.$

Note that the left-hand side is effectively the local truncation error multiplied by ${\varepsilon}^2$ . Here, we first substitute in for ${\mathcal{{R}}}^{\rho u, n}_j$ using Eq (3.22), resulting in the form

$\frac{{\varepsilon}^2}{{{\Delta t}}}{\left( {{\left( {\bar{\rho u}} \right)}^n_j - {{\Delta t}} D_x\widetilde{F}^{\rho u}({\mathit{\boldsymbol{U}}}_j^n) - \frac{a^n {{\Delta t}}}{{\varepsilon}^2}D_x\bar\rho^{n+1}_j-\Psi^n_j{\left( {\bar{\rho u}} \right)}^{n+1}_j } \right)} = \frac{{\varepsilon}^2}{{{\Delta t}}}\cdot C_1\frac{{{\Delta t}} ({{\Delta x}})^2}{{\varepsilon}^2},$

where $C_1$ comes from the leading order of diffusion ${\mathcal{D}}_{{j+\frac{1}{2}}}^{\rho u, n}$ that appears within the ${\mathcal{O}}({\varepsilon}^{-2}({{\Delta x}})^2)$ term of Eq (3.22). Substituting in $\widetilde{F}^{\rho u}({\mathit{\boldsymbol{U}}}_j^n)$ from Eq (3.3) and $\Psi_j^n$ from Eq (3.10), we obtain

${\varepsilon}^2{\left[ {\frac{{\left( {\bar{\rho u}} \right)}^n_j-{\left( {\bar{\rho u}} \right)}^{n+1}_j}{{{\Delta t}}}-D_x{\left( {\rho u^2} \right)}_j^n} \right]} - D_x{\left( {p(\bar\rho_j^n)-a^n \bar\rho_j^n} \right)} - a^n D_x\bar\rho^{n+1}_j-\frac{C_\delta \kappa}{2}{\left( {\bar{\rho u}} \right)}^{n+1}_j {{\left| {u_j^n} \right|}} = C_1\cdot({{\Delta x}})^2.$

We then apply the asymptotic expansions from Eq (3.28) and take ${\varepsilon} \rightarrow 0$ , resulting in

$D_x {\left[ {{\left( {\bar\rho_j^{(0), n}} \right)}^\gamma} \right]}+\frac{C_\delta \kappa}{2}\bar{\rho}^{(0)n+1}_j u^{(0), n+1}_j {{\left| {u_j^{(0), n}} \right|}} = a^n D_x{\left( {\bar\rho^{(0), n}_j-\bar\rho^{(0), n+1}_j} \right)} - C_1\cdot({{\Delta x}})^2.$

Since $D_x$ represents the second-order central difference operator, we equivalently have

${\left[ {{\left( {\bar\rho_j^{(0), n}} \right)}^\gamma} \right]}_x+\frac{C_\delta \kappa}{2}\bar{\rho}^{(0), n+1}_j u^{(0), n+1}_j {{\left| {u_j^{(0), n}} \right|}} = a^n {\left( {\bar\rho^{(0), n}_j-\bar\rho^{(0), n+1}_j} \right)}_x +(C_2 - C_1)({{\Delta x}})^2,$

where $C_2$ represents the the leading error coefficient from using central difference. Lastly, applying Taylor series to the terms at time $t^n$ about time $t^{n+1}$ , we obtain

${\left[ {{\left( {\bar\rho_j^{(0), n+1}} \right)}^\gamma} \right]}_x+\frac{C_\delta \kappa}{2}\bar{\rho}^{(0), n+1}_j u^{(0), n+1}_j{{\left| {u_j^{(0), n+1}} \right|}} = C_3{{\Delta t}}+(C_2 - C_1)({{\Delta x}})^2.$

Here, $C_3$ denotes the leading error term from the Taylor series. Therefore, we can bound the asymptotic expansion at time $t^{n+1}$ :

$\begin{equation*} {\left\| {{\left[ {{\left( {\bar\rho_j^{(0), n+1}} \right)}^\gamma} \right]}_x+\frac{C_\delta \kappa}{2}\bar{\rho}^{(0), n+1}_j u^{(0), n+1}_j{{\left| {u_j^{(0), n+1}} \right|}}} \right\|} \leq C({{\Delta t}}+({{\Delta x}})^2), \end{equation*}$

with $C = \max{\left\{ {{\left| {C_3} \right|}, {\left| {C_2-C_1} \right|}} \right\}}$ , hence, implying that the proposed method preserves the asymptotic state up to the order of the scheme.

The limiting scheme is accompanied by the (limiting) coupling conditions (3.24) at the boundary implemented again using ghost cells. Note again that by the discussion made in §3.5.2, the pipe entries and exits occurring at intersections would have a bound that is instead first-order in both space and time.

3.6.2. Full discrete scheme at the low mach/high friction limit

We now present the reduced form of the proposed scheme shown in Eqs (3.16) and (3.17). As we will show, this reduced form properly approximates the system (2.1) at the limit ${\varepsilon} \rightarrow 0$ , which reads

$\begin{equation} \begin{aligned} & \rho_t + {\left( {\rho u} \right)}_x = 0, \\ & p_x = -\frac{C_\delta \kappa}{2}\rho u {\left| {u} \right|}. \end{aligned} \end{equation}$

(3.29)

Prior to the limiting discretization, we would like to note that in the actual numerical implementation, any division of ${\varepsilon}$ is replaced with $\max{\left\{ {{\varepsilon}, \epsilon^m} \right\}}$ , with $\epsilon^m$ denoting machine precision, to avoid any division-by-zero errors. Furthermore, all numerical divisions by ${\varepsilon} = 0$ are safe, including that within $\Psi_j^n$ in Eq (3.10). We bring this up since the limit

$\lim\limits_{{\varepsilon} \rightarrow 0}\ \frac{1}{{\varepsilon}^2\Psi_j^n} = \frac{2}{{{\Delta t}} C_\delta \kappa |u_j^n|},$

arises often in the limiting scheme—and although the above appears problematic, the true numerical implementation of the above limit reads

$\lim\limits_{{\varepsilon} \rightarrow 0}\ \frac{1}{{\varepsilon}^2\Psi_j^n} = {\left( {\epsilon^m + {{\Delta t}}\frac{C_\delta \kappa}{2}|u_j^n|} \right)}^{-1},$

and is free of issues when ${\varepsilon} = 0$ or ${\left| {u_j^n} \right|} = 0$ . For the sake of computing the limiting scheme, we will proceed as if ${\varepsilon} = 0$ is not replaced by the machine precision parameter $\epsilon^m$ , and continue knowing any division of ${\left| {u_j^n} \right|}$ does indeed avoid division-by-zero errors.

We start by looking at the limiting scheme for $(\rho u)_j^{n+1}$ written in Eq (3.17). After applying the above limit and simplifying the numerical fluxes ${\mathcal{{R}}}_j^{\rho u, n}$ when ${\varepsilon} = 0$ , we obtain

$\begin{equation} {\left( {\bar{\rho u}} \right)}_j^{n+1} = \frac{2}{C_\delta \kappa |u_j^n|}{\left[ {-D_x p_j^n + a^n D_x{\left( {\bar\rho_j^{n+1} - \bar\rho_j^n} \right)} + \frac{1}{{{\Delta x}}}{\left( {{\left( {{\mathcal{{D}}}_0^{\rho u}} \right)}_{{j-\frac{1}{2}}}^n - {\left( {{\mathcal{{D}}}_0^{\rho u}} \right)}_{{j+\frac{1}{2}}}^n} \right)} } \right]}, \end{equation}$

(3.30)

where $D_x$ again denotes the central difference operator, and ${\left({{\mathcal{{D}}}_0^{\rho u}} \right)}_{{j+\frac{1}{2}}}^n$ denotes the diffusion term from Eq (3.19) for the $\rho u$ numerical flux at the limit ${\varepsilon} \rightarrow 0$ , and reads

$\begin{equation*} \begin{aligned} {\left( {{\mathcal{{D}}}_0^{\rho u}} \right)}_{{j+\frac{1}{2}}}^n = & \frac{s^+_{{j+\frac{1}{2}}}}{2{\left( {s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} \right)}}{\left[ {2 {\left( {p^-_{{j+\frac{1}{2}}} - a^n\rho^-_{{j+\frac{1}{2}}}} \right)}-{\left( {p^n_j - a^n\bar\rho^n_j} \right)}-{\left( {p^n_{j+1} - a^n\bar\rho^n_{j+1}} \right)}} \right]}\\[4pt] &-\frac{s^-_{{j+\frac{1}{2}}}}{2{\left( {s^+_{{j+\frac{1}{2}}}-s^-_{{j+\frac{1}{2}}}} \right)}}{\left[ {2 {\left( {p^+_{{j+\frac{1}{2}}} - a^n\rho^+_{{j+\frac{1}{2}}}} \right)}-{\left( {p^n_j - a^n\bar\rho^n_j} \right)}-{\left( {p^n_{j+1} - a^n\bar\rho^n_{j+1}} \right)}} \right]}. \end{aligned} \end{equation*}$

We again can see that the limiting scheme (3.30) does indeed satisfy, with first-order accuracy, the limiting equation for $p_x$ seen in Eq (3.29) and, thus, the AP property. The limiting scheme for $\rho_j^{n+1}$ looks fairly similar to the original scheme written in Eq (3.16), as the majority of terms do not contain ${\varepsilon}$ . The resulting limiting scheme for $\rho_j^{n+1}$ reads

$\begin{equation} \begin{aligned} &\bar{\rho}^{n+1}_j - \frac{{{\Delta t}}}{({{\Delta x}})^2}a^n{\left[ {(\phi_0)^n_{{j+\frac{1}{2}}} \bar{\rho}^{n+1}_{j+1} - {\left( {(\phi_0)^n_{{j+\frac{1}{2}}} + (\phi_0)^n_{{j-\frac{1}{2}}}} \right)}\bar{\rho}^{n+1}_j + (\phi_0)^n_{{j-\frac{1}{2}}} \bar{\rho}^{n+1}_{j-1}} \right]}\\ & \quad = \bar{\rho}^n_j + {{\Delta t}} D_x {\left[ { \frac{2}{C_\delta \kappa |u_j^n|}{\left( {D_x p_j^n + a^n D_x\bar\rho_j^n - \frac{1}{{{\Delta x}}}{\left( {{\left( {{\mathcal{{D}}}_0^{\rho u}} \right)}_{{j-\frac{1}{2}}}^n - {\left( {{\mathcal{{D}}}_0^{\rho u}} \right)}_{{j+\frac{1}{2}}}^n} \right)} } \right)} } \right]}, \end{aligned} \end{equation}$

(3.31)

where

$(\phi_0)^n_{{j+\frac{1}{2}}} = \frac{1}{C_\delta \kappa}{\left( {\frac{1}{|u_j^n|} + \frac{1}{|u_{j+1}^n|}} \right)}.$

While it is unclear in the limiting scheme, one can show Eq (3.31) satisfies the following equation up to the first-order consistency error:

$\rho_t = {\left( {\frac{2}{C_\delta \kappa |u|}p_x} \right)}_x,$

which is equivalent to Eq (3.29) after a substitution.

4. Numerical results

In this section, we demonstrate the accuracy and performance of the proposed AP scheme on pipe networks in several numerical experiments conducted on the isentropic Euler equations in (2.1). For all examples, we follow the CFL condition in (3.18) and take the CFL number $\nu = 0.45$ ; in the source term of (2.1), choose $C_\delta = 1$ ; take the minmod parameter $\theta = 1.3$ ; and set the tolerance for the Newton solve at the junctions to be $10^{-8}$ .

We start with a test on a single pipe with zero friction in Example 1 to compare to previous work. In Examples 2–4, we will conduct verifications on the two types of T-junctions—the so-called 1-to-2 T-junction, which has one pipe entering and two pipes leaving the junction, and the so-called 2-to-1 T-junction, which has two pipes entering and one pipe leaving the junction. Illustrations of these are presented in Figure 2. While we only consider the example T-junction networks with all pipes having the same length, the method is of course not restricted to this and can be extended to large and more complex pipe networks, such as those in which there are different speed regimes in different pipe segments.

Figure 2. Illustrations of the 1-to-2 T-junction (left) and 2-to-1 T-junction (right).

DownLoad: Full-Size Img PowerPoint

Note that according to ^[26], $\alpha$ from Eq (3.5) must be $0 < \alpha < 1$ . Thus, in any ${\varepsilon} = 1$ case considered in the following examples, we follow ^[26] and take $\alpha = 0.5$ instead of using $\alpha = {\varepsilon}^{b}$ as in Eq (3.5).

4.1. Example 1—Shocktube problem

In this first example, we look at a single pipe with no junctions and turn off friction; i.e., we set $\kappa = 0$ to confirm the proposed method against existing AP methods. Note that in the presence of no friction and no junctions, the proposed method simplifies to that of ^[26]. We follow the initial conditions of the first example of ^[25], which considers the Riemann problem

$(\rho u)(x, 0) = 1, \qquad \rho(x, 0) = \left\{ \begin{aligned} &1+ {\varepsilon}^2 && x < 0.5 \\ & 1 && \text{otherwise}, \end{aligned}\right.$

on the domain [0, 1] with $\gamma = 1.4$ and Neumann boundary conditions. Note that since ^[25] defined their system slightly different than that in Eq (2.1), ${\varepsilon}$ in this paper corresponds to ${\varepsilon}^2$ in ^[25].

We compare results to the exact solution of the Riemann problem for ${\varepsilon} = 1, 0.1,$ and 0.01 at the end times of $t = 0.125, 0.02,$ and 2.5 $\times10^{-3}$ , respectively, in . The coarser mesh sizes are selected, i.e., taking 50 discretization cells for ${\varepsilon} = 1$ , 125 cells for ${\varepsilon} = 0.1$ , and 500 cells for ${\varepsilon} = 0.01$ , are taken to match those in ^[25], and the refined meshes are included for further comparison. As expected, the results shown in are more diffusive as ${\varepsilon}$ approaches zero. This is due to the fact that the stiff portion of the system (3.1), computed with an implicit method, significantly outweighs the non-stiff portion of the system as ${\varepsilon} \rightarrow 0$ .

Figure 3. Example 1: Solutions to

$\rho$ (top row) and

$\rho u$ (bottom row) for the Riemann problem with

${\varepsilon} = 1$ at

$t = 0.125$ (left column),

${\varepsilon} = 0.1$ at

$t = 0.02$ (middle column), and

${\varepsilon} = 0.01$ at

$t = 2.5\times 10^{-3}$ (right column).

DownLoad: Full-Size Img PowerPoint

4.2. Example 2—Convergence tests

In this first example including junctions, we consider a problem with a smooth density and velocity profile to confirm the convergence rate of the proposed AP method. We consider both the 1-to-2 and 2-to-1 T-junction setups with $\gamma = 5/3$ and $\kappa = 10^{-3}$ . For ingoing pipelines, we take the initial conditions

$u(x, 0) = 0, \qquad \rho(x, 0) = \left\{ \begin{aligned} &1.1 && x \leq \frac{0.4}{{\varepsilon}} \\ &1+0.1\sin{\left( {\frac{\pi{\varepsilon}}{0.8}x} \right)} && \frac{0.4}{{\varepsilon}} < x < \frac{0.8}{{\varepsilon}}, \\ & 1 && x > \frac{0.8}{{\varepsilon}}, \end{aligned}\right.$

and for outgoing pipes we take the initial conditions $\rho(x, 0) = 1$ , $u(x, 0) = 0$ . We will take each pipe to have a length of ${\varepsilon}^{-1}$ , as we are more interested in the investigating behavior near the junction than the inlet/outlet boundaries. These pipe lengths allow for the smooth profile to pass through the junction before the final time of $t = 0.2$ , but the smooth profile will not reach the inlets/outlets by final time. At the inlets, we take a Dirichlet boundary condition of $\rho = 1.1$ , and at the outlets we take Neumann boundary conditions.

Since the exact solution is unknown, we obtain the experimental $L^1$ convergence rates by computing the solution on a number of different meshes, then use the Runge formula

$\begin{equation*} {\rm Rate}_{{\Delta x}}(\rho) = \log_2{\left( {\frac{\|\rho(x, 1)_{{{\Delta x}}}-\rho(x, 1)_{2{{\Delta x}}}\|_1}{\|\rho(x, 1)_{{{\Delta x}}/2}-\rho(x, 1)_{{\Delta x}}\|_1}} \right)}, \end{equation*}$

where $\rho_{{\Delta x}}$ denotes the density solution computed on a uniform mesh with all cells having width ${{\Delta x}}$ . The convergence rates for $u$ can be computed analogously. We present the $L^1$ errors and rates for $\rho$ and $u$ for the 1-to-2 junction in and for the 2-to-1 junction in , both of which confirm the expected first order of accuracy when the fluid passes through the junction. Note that for the ${\varepsilon} = 1$ case, we use a more refined grid to indeed confirm the first-order convergence for $u$ . In addition, note that as ${\varepsilon} \rightarrow 0$ , there seems to be a partially higher order convergence rate in and . However, this appearance is only due to the evolution of $u$ being near zero everywhere (in this case, starting at $u = 0$ ). In the case when $u = 0$ , one can show that the resulting asymptotic scheme in Eq (3.31) approaches second-order accuracy, but this is only specific to small $u$ and is not expected otherwise.

Table 1. Example 2 (1-to-2 T-junction):

$L^1$ errors and corresponding experimental convergence rates.

	${{\Delta x}}$	$\\|\rho(x, 1)_{{{\Delta x}}/2}-\rho(x, 1)_{{\Delta x}}\\|_1$	Rate $_{{\Delta x}}(\rho)$	$\\|u(x, 1)_{{{\Delta x}}/2}-u(x, 1)_{{\Delta x}} \\|_1$	Rate $_{{\Delta x}}(u)$
${\varepsilon} = 1$	1/80	3.54e-04	–	2.44e-04	–
	1/160	1.82e-04	0.96	1.46e-04	0.75
	1/320	9.29e-05	0.97	8.25e-05	0.82
	1/640	4.65e-05	1.00	4.36e-05	0.92
	1/1280	2.32e-05	1.00	2.27e-05	0.94
${\varepsilon} = 0.1$	1/10	1.43e-02	–	1.39e-01	–
	1/20	7.72e-03	0.89	7.57e-02	0.88
	1/40	3.89e-03	0.99	3.93e-02	0.94
	1/80	1.97e-03	0.98	2.05e-02	0.94
	1/160	9.85e-04	1.00	1.05e-02	0.96
${\varepsilon} = 0.01$	1/10	8.59e-02	–	1.20e+01	–
	1/20	3.08e-02	1.48	3.43e+00	1.81
	1/40	9.53e-03	1.69	1.08e+00	1.67
	1/80	2.82e-03	1.76	3.08e-01	1.80
	1/160	9.65e-04	1.55	9.94e-02	1.63
${\varepsilon} = 0.001$	1/10	2.89e-02	–	8.21e+00	–
	1/20	9.10e-03	1.67	1.38e+00	2.58
	1/40	3.06e-03	1.57	5.28e-01	1.38
	1/80	1.01e-03	1.60	2.22e-01	1.25
	1/160	3.64e-04	1.47	1.04e-01	1.09

| Show Table

DownLoad: CSV

Table 2. Example 2 (2-to-1 T-junction):

$L^1$ errors and corresponding experimental convergence rates.

	${{\Delta x}}$	$\\|\rho(x, 1)_{{{\Delta x}}/2}-\rho(x, 1)_{{\Delta x}}\\|_1$	Rate $_{{\Delta x}}(\rho)$	$\\|u(x, 1)_{{{\Delta x}}/2}-u(x, 1)_{{\Delta x}} \\|_1$	Rate $_{{\Delta x}}(u)$
${\varepsilon} = 1$	1/80	3.44e-04	–	2.65e-04	–
	1/160	1.76e-04	0.97	1.54e-04	0.78
	1/320	9.02e-05	0.96	8.84e-05	0.80
	1/640	4.53e-05	0.99	4.70e-05	0.91
	1/1280	2.27e-05	0.99	2.47e-05	0.93
${\varepsilon} = 0.1$	1/10	1.49e-02	–	1.41e-01	–
	1/20	6.98e-03	1.10	8.33e-02	0.76
	1/40	3.35e-03	1.06	4.42e-02	0.91
	1/80	1.67e-03	1.01	2.27e-02	0.96
	1/160	8.42e-04	0.99	1.14e-02	0.99
${\varepsilon} = 0.01$	1/10	9.23e-02	–	1.19e+01	–
	1/20	3.84e-02	1.26	3.01e+00	1.98
	1/40	1.19e-02	1.69	9.56e-01	1.66
	1/80	2.99e-03	1.99	3.12e-01	1.61
	1/160	9.44e-04	1.66	1.21e-01	1.37
${\varepsilon} = 0.001$	1/10	2.93e-02	–	8.16e+00	–
	1/20	9.33e-03	1.65	1.35e+00	2.60
	1/40	3.18e-03	1.55	5.12e-01	1.40
	1/80	1.07e-03	1.57	2.11e-01	1.27
	1/160	4.02e-04	1.42	9.78e-02	1.11

| Show Table

DownLoad: CSV

4.3. Example 3—Inlet discontinuity

For the next example, we look at the evolution of an initial jump discontinuity that starts at the inlets, resulting in a wave that travels through the junction. For all pipelines, we consider the initial conditions

$\rho(x, 0) = 1, \qquad u(x, 0) = 0,$

with $\gamma = 5/3$ , $\kappa = 10^{-3}$ , and each pipe having a length of 100. At the inlets, we take a Dirichlet boundary condition of $\rho = 1.3$ , and prescribe Neumann boundary conditions at the outlets. We look at this initial set up for both the 1-to-2 and 2-to-1 T-junction networks.

We use the proposed AP scheme to compute the solution for ${\varepsilon}$ values of $0.1$ , $0.01$ , and $0.001$ to the final times of $t = 10$ , $t = 1$ , and $t = 0.1$ , respectively, to capture the behavior of the discontinuity through the junction. The resulting solutions for $\rho$ and $u$ when taking a spatial mesh of 4000 cells in each pipeline are presented in for the 1-to-2 T-junction and in for the 2-to-1 T-junction. In these figures, the influence of the friction source term is very clear, as it appears as if the flow did not yet reach the junction in the ${\varepsilon} = 0.001$ results. This is, however, not the case and is only a result of the strong friction; see , where we plot the outgoing pipeline results for ${\varepsilon} = 0.001$ at $t = 0.1$ . In addition, we see that in the cases that produce discontinuities in the solution, there are no spurious oscillations present.

Figure 4. Example 3 (1-to-2 T-junction): Solutions for

$\rho$ (left column) and

$u$ (right column) for

${\varepsilon} = 0.1$ at

$t = 10$ (top row),

${\varepsilon} = 0.01$ at

$t = 1$ (middle row), and

${\varepsilon} = 0.001$ at

$t = 0.1$ (bottom row).

DownLoad: Full-Size Img PowerPoint

Figure 5. Example 3 (2-to-1 T-junction): Solutions for

$\rho$ (left column) and

$u$ (right column) for

${\varepsilon} = 0.1$ at

$t = 10$ (top row),

${\varepsilon} = 0.01$ at

$t = 1$ (middle row), and

${\varepsilon} = 0.001$ at

$t = 0.1$ (bottom row).

DownLoad: Full-Size Img PowerPoint

Figure 6. Example 3: Results for

$\rho$ (left) and

$u$ (right) of the outgoing pipeline(s) for the 1-to-2 T-junction (blue) and the 2-to-1 T-junction (red dotted).

DownLoad: Full-Size Img PowerPoint

Lastly, for this example, we run the simulation to long times to ensure we obtain the expected behavior of density approaching a constant state. Due to friction, the steady state is $p_x = -\frac{C_\delta \kappa}{2 {\varepsilon}^2} \rho u {\left| {u} \right|}$ and takes quite long to converge. Therefore, we show this by plotting the densities for the 1-to-2 setup at multiple times ( $t = 20,200,500$ ) for cases ${\varepsilon} = 0.01$ , ${\varepsilon} = 0.001$ , which are presented in , and clearly present a convergence to the constant state $\rho = 1.3$ as time approaches infinity.

Figure 7. Example 3 (1-to-2 T-junction): Solutions for

$\rho$ for the

${\varepsilon} = 0.01$ case (left column) and the

${\varepsilon} = 0.001$ case (right column) at times

$t = 20$ (top row),

$t = 200$ (middle row), and

$t = 500$ (bottom row).

DownLoad: Full-Size Img PowerPoint

4.4. Example 4—Comparison with the CU explicit scheme

In this example, we test the importance of the ${\varepsilon}$ -independent time-step restriction and compare the proposed AP scheme against the CU finite volume discretization with a forward Euler time-step. The forward Euler time-step was selected so that the explicit scheme matches the first-order in time, second-order in space accuracy of the proposed AP scheme.

We consider the 1-to-2 T-junction with the same $\gamma$ , $\kappa$ , initial conditions, and boundary conditions as that of Example 3, and run both the proposed AP scheme and the related explicit scheme to a final time $t = 10$ for the ${\varepsilon}$ values of 0.1, 0.01, and 0.001 on various spatial discretizations. We present the run times of both schemes in . In the run times, we see that for ${\varepsilon} = 0.1$ , the run times for the two methods are on the same order of magnitude. However, the difference becomes drastic as ${\varepsilon} \rightarrow 0$ . It is clear in that the AP scheme keeps all run times roughly the same for varying ${\varepsilon}$ . Conversely, the explicit scheme has the expected ${\varepsilon}$ -dependence within the experimental run times—which, for ${\varepsilon} = 0.001$ , results in simulations that are a daunting 140 $\times$ slower than that of the AP scheme.

Table 3. Example 4: Run time comparisons of the proposed AP scheme and CU discretization with forward Euler time-stepping for various

${\varepsilon}$ and mesh sizes.

	${{\Delta x}}$	AP Scheme Run Time	Explicit Scheme Run Time
${\varepsilon} = 0.1$	1/20	2.44 s	4.05 s
	1/40	9.14 s	16.1 s
	1/80	37.4 s	65.5 s
${\varepsilon} = 0.01$	1/20	2.89 s	46.8 s
	1/40	11.0 s	192 s
	1/80	44.0 s	782 s
${\varepsilon} = 0.001$	1/20	2.91 s	406 s
	1/40	11.0 s	1610 s
	1/80	43.7 s	6450 s

| Show Table

DownLoad: CSV

5. Conclusion

In this study, we developed a novel asymptotic-preserving method for the isentropic Euler equations with a friction source on pipe networks. As a direct result, this would allow the simulation of full network—possibly containing flow in different regimes in different segments—to be computed in its entirety using solely the proposed method.··· To address the difficulties and stiffness of the system in low Mach or high friction regimes, we split the flux into pieces that represent the slow and fast dynamics. This allows for an explicit time-step on the non-stiff (slow) dynamics portion, and, by using a method based on Rosenbrock-type Runge-Kutta schemes, allows for an implicit time-step for the stiff (fast) dynamics that do not require a nonlinear solver. In turn, we were able to confirm both experimentally and theoretically that the proposed scheme is AP in the sense that it provides a consistent and stable solution in low Mach and high friction regimes. Most importantly, the CFL stability restriction is only dependent on mesh size, rather than requiring dependence on both the mesh and the small limiting parameter ${\varepsilon}$ within the system. This method within each pipeline is then extended to entire pipe networks, in which coupling conditions must be used at pipe-to-pipe intersections to ensure a mathematically well-posed problem. We show that, even in the limiting regime, the coupling conditions remain well-posed. These coupling conditions are used to set the ghost cell values on each pipeline, thus allowing a seamless coupling of the AP method across pipe junctions within the network.

Since the AP method is combined with a friction source term and pipe junctions, both the asymptotic state in the low Mach/high friction regime and the boundary conditions are nontrivial. Because of this, in addition to a first-order approximation in time, the scheme is currently limited to first-order. The second-order extension to the proposed method is, to our knowledge, nontrivial, as it would require, e.g., one-sided reconstructions, adjustments to the half-Riemann problems, and approximating fluxes in the ghost cells. In future work, these difficulties are something we plan to address in hopes of making a fully second-order AP scheme for isentropic flow in pipe networks. In addition, there is a question in the consideration of making the scheme well-balanced; i.e., preserving the discrete steady states of the system exactly, in addition to asymptotic-preserving property. However, due to the flux splitting, and different spatial discretizations within the stiff and non-stiff portions, we believe this to be nontrivial. Thus, ensuring the combination of the well-balanced and asymptotic-preserving property for pipe-networks is also planned to be addressed in future works.

Author contributions

Michael Redle: Project administration, Conceptualization, Investigation, Methodology, Formal analysis, Software, Validation, Visualization, Writing—original draft. Michael Herty: Funding Acquisition, Conceptualization, Investigation, Methodology, Writing—original draft.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for the financial support through 463312734/FOR5409, 320021702/GRK2326, and within SPP 2410 Hyperbolic Balance Laws in Fluid Mechanics: Complexity, Scales, Randomness (CoScaRa) within the Project(s) HE5386/26-1 (Numerische Verfahren für gekoppelte Mehrskalenprobleme, 525842915) and (Zufällige kompressible Euler Gleichungen: Numerik und ihre Analysis, 525853336) HE5386/27-1. Support received funding from the European Union's Horizon Europe research and innovation programme under the Marie Sklodowska-Curie Doctoral Network Datahyking (Grant No. 101072546) is acknowledged.

Conflict of interest

The authors declare there is no conflict of interest.

References

[1]	A. Osiadacz, Simulation and analysis of gas networks, 1987. Available from: https://www.osti.gov/biblio/5141539.
[2]	A. Bressan, S. Čanić, M. Garavello, M. Herty, B. Piccoli, Flows on networks: Recent results and perspectives, EMS Surv. Math. Sci., 1 (2014), 47–111. https://doi.org/10.4171/emss/2 doi: 10.4171/emss/2
[3]	J. Brouwer, I. Gasser, M. Herty, Gas pipeline models revisited: Model hierarchies, nonisothermal models, and simulations of networks, Multiscale Model. Simul., 9 (2011), 601–623. https://doi.org/10.1137/100813580 doi: 10.1137/100813580
[4]	H. Guillard, C. Viozat, On the behaviour of upwind schemes in the low Mach number limit, Comput. Fluids, 28 (1999), 63–86. https://doi.org/10.1016/S0045-7930(98)00017-6 doi: 10.1016/S0045-7930(98)00017-6
[5]	H. Guillard, A. Murrone, On the behavior of upwind schemes in the low Mach number limit: Ⅱ. Godunov type schemes, Comput. Fluids, 33 (2004), 655–675. https://doi.org/10.1016/j.compfluid.2003.07.001 doi: 10.1016/j.compfluid.2003.07.001
[6]	F. Rieper, On the dissipation mechanism of upwind-schemes in the low Mach number regime: A comparison between roe and hll, J. Comput. Phys., 229 (2010), 221–232. https://doi.org/10.1016/j.jcp.2009.09.043 doi: 10.1016/j.jcp.2009.09.043
[7]	S. Jin, Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms, J. Comput. Phys., 122 (1995), 51–67. https://doi.org/10.1006/jcph.1995.1196 doi: 10.1006/jcph.1995.1196
[8]	S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations, SIAM J. Sci. Comput., 21 (1999), 441–454. https://doi.org/10.1137/S1064827598334599 doi: 10.1137/S1064827598334599
[9]	J. Hu, S. Jin, Q. Li, Asymptotic-preserving schemes for multiscale hyperbolic and kinetic equations, In: Handbook of Numerical Analysis, Elsevier, 18 (2007), 103–129. https://doi.org/10.1016/bs.hna.2016.09.001
[10]	J. Hu, R. Shu, X. Zhang, Asymptotic-preserving and positivity-preserving implicit-explicit schemes for the stiff BGK equation, SIAM J. Numer. Anal., 56 (2018), 942–973. https://doi.org/10.1137/17M1144362 doi: 10.1137/17M1144362
[11]	S. Jin, L. Pareschi, G. Toscani, Uniformly accurate diffusive relaxation schemes for multiscale transport equations, SIAM J. Numer. Anal., 38 (2000), 913–936. https://doi.org/10.1137/S0036142998347978 doi: 10.1137/S0036142998347978
[12]	S. Jin, Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: A review, Riv. Math. Univ. Parma (N.S.), 3 (2012), 177–216.
[13]	S. Jin, Asymptotic-preserving schemes for multiscale physical problems, Acta Numer., 31 (2022), 415–489. https://doi.org/10.1017/S0962492922000010 doi: 10.1017/S0962492922000010
[14]	W. Ren, H. Liu, S. Jin, An asymptotic-preserving Monte Carlo method for the Boltzmann equation, J. Comput. Phys., 276 (2014), 380–404. https://doi.org/10.1016/j.jcp.2014.07.029 doi: 10.1016/j.jcp.2014.07.029
[15]	B. Zhang, H. Liu, S. Jin, An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation, J. Comput. Phys., 305 (2016), 575–588. https://doi.org/10.1016/j.jcp.2015.11.006 doi: 10.1016/j.jcp.2015.11.006
[16]	F. Fei, A time-relaxed Monte Carlo method preserving the Navier-Stokes asymptotics, J. Comput. Phys., 486 (2023), 112128. https://doi.org/10.1016/j.jcp.2023.112128 doi: 10.1016/j.jcp.2023.112128
[17]	F. Fei, A Navier-Stokes asymptotic preserving Direct Simulation Monte Carlo method for multi-species gas flows, arXiv preprint arXiv: 2410.20322, 2024.
[18]	S. Jin, Z. Ma, K. Wu, Asymptotic-preserving neural networks for multiscale time-dependent linear transport equations, J. Sci. Comput., 94 (2023), 57. https://doi.org/10.1007/s10915-023-02100-0 doi: 10.1007/s10915-023-02100-0
[19]	S. Jin, Z. Ma, K. Wu, Asymptotic-preserving neural networks for multiscale kinetic equations, Commun. Comput. Phys., 35 (2024), 693–723. https://doi.org/10.4208/cicp.oa-2023-0211 doi: 10.4208/cicp.oa-2023-0211
[20]	K. R. Arun, S. Samantaray, An asymptotic preserving time integrator for low Mach number limits of the Euler equations with gravity, arXiv preprint arXiv: 1902.00221, 2019.
[21]	K. R. Arun, S. Samantaray, Asymptotic preserving low Mach number accurate IMEX finite volume schemes for the isentropic Euler equations, J. Sci. Comput., 82 (2020), 1–32. https://doi.org/10.1007/s10915-020-01138-8 doi: 10.1007/s10915-020-01138-8
[22]	S. Boscarino, J. M. Qiu, G. Russo, T. Xiong, A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system, J. Comput. Phys., 392 (2019), 594–618. https://doi.org/10.1016/j.jcp.2019.04.057 doi: 10.1016/j.jcp.2019.04.057
[23]	P. Degond, M. Tang, All speed scheme for the low Mach number limit of the isentropic Euler equations, Commun. Comput. Phys., 10 (2011), 1–31. https://doi.org/10.4208/cicp.210709.210610a doi: 10.4208/cicp.210709.210610a
[24]	G. Dimarco, R. Loubère, M. H. Vignal, Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit, SIAM J. Sci. Comput., 39 (2017), A2099–A2128. https://doi.org/10.1137/16M1069274 doi: 10.1137/16M1069274
[25]	G. Dimarco, R. Loubère, V. Michel-Dansac, M. H. Vignal, Second-order implicit-explicit total variation diminishing schemes for the Euler system in the low Mach regime, J. Comput. Phys., 372 (2018), 178–201. https://doi.org/10.1016/j.jcp.2018.06.022 doi: 10.1016/j.jcp.2018.06.022
[26]	J. Haack, S. Jin, J. G. Liu, An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations, Commun. Comput. Phys., 12 (2012), 955–980. https://doi.org/10.4208/cicp.250910.131011a doi: 10.4208/cicp.250910.131011a
[27]	S. Samantaray, Asymptotic preserving linearly implicit additive IMEX-RK finite volume schemes for low Mach number isentropic Euler equations, arXiv preprint arXiv: 2409.05854, 2024.
[28]	S. Avgerinos, F. Bernard, A. Iollo, G. Russo, Linearly implicit all Mach number shock capturing schemes for the Euler equations, J. Comput. Phys., 393 (2019), 278–312. https://doi.org/10.1016/j.jcp.2019.04.020 doi: 10.1016/j.jcp.2019.04.020
[29]	G. Bispen, M. Lukáčová-Medvi $\check{{\text{d}}}$ ová, L. Yelash, Asymptotic preserving IMEX finite volume schemes for low Mach number Euler equations with gravitation, J. Comput. Phys., 335 (2017), 222–248. https://doi.org/10.1016/j.jcp.2017.01.020 doi: 10.1016/j.jcp.2017.01.020
[30]	R. Klein, Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics Ⅰ: One-dimensional flow, J. Comput. Phys., 121 (1995), 213–237. https://doi.org/10.1016/S0021-9991(95)90034-9 doi: 10.1016/S0021-9991(95)90034-9
[31]	V. Kučera, M. Lukáčová-Medvi $\check{{\text{d}}}$ ová, S. Noelle, J. Schütz, Asymptotic properties of a class of linearly implicit schemes for weakly compressible Euler equations, Numer. Math., 150 (2022), 79–103. https://doi.org/10.1007/s00211-021-01240-5 doi: 10.1007/s00211-021-01240-5
[32]	S. Noelle, G. Bispen, K. R. Arun, M. Lukáčová-Medvi $\check{{\text{d}}}$ ová, C. D. Munz, A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics, SIAM J. Sci. Comput., 36 (2014), B989–B1024. https://doi.org/10.1137/120895627 doi: 10.1137/120895627
[33]	J. Zeifang, J. Schütz, K. Kaiser, A. Beck, M. Lukáčová-Medvi $\check{{\text{d}}}$ ová, S. Noelle, A novel full-Euler low Mach number IMEX splitting, Commun. Comput. Phys., 27 (2020), 292–320. https://doi.org/10.4208/cicp.oa-2018-0270 doi: 10.4208/cicp.oa-2018-0270
[34]	W. Boscheri, M. Dumbser, M. Ioriatti, I. Peshkov, E. Romenski, A structure-preserving staggered semi-implicit finite volume scheme for continuum mechanics, J. Comput. Phys., 424 (2021), 109866. https://doi.org/10.1016/j.jcp.2020.109866 doi: 10.1016/j.jcp.2020.109866
[35]	F. Cordier, P. Degond, A. Kumbaro, An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations, J. Comput. Phys., 231 (2012), 5685–5704. https://doi.org/10.1016/j.jcp.2012.04.025 doi: 10.1016/j.jcp.2012.04.025
[36]	M. Dumbser, V. Casulli, A conservative, weakly nonlinear semi-implicit finite volume scheme for the compressible Navier-Stokes equations with general equation of state, Appl. Math. Comput., 272 (2016), 479–497. https://doi.org/10.1016/j.amc.2015.08.042 doi: 10.1016/j.amc.2015.08.042
[37]	I. Peshkov, M. Dumbser, W. Boscheri, E. Romenski, S. Chiocchetti, M. Ioriatti, Simulation of non-Newtonian viscoplastic flows with a unified first order hyperbolic model and a structure-preserving semi-implicit scheme, Comput. Fluids, 224 (2021), 104963. https://doi.org/10.1016/j.compfluid.2021.104963 doi: 10.1016/j.compfluid.2021.104963
[38]	G. Bispen, K. R. Arun, M. Lukáčová-Medvi $\check{{\text{d}}}$ ová, S. Noelle, IMEX large time step finite volume methods for low Froude number shallow water flows, Commun. Comput. Phys., 16 (2014), 307–347. http://doi.org/10.4208/cicp.040413.160114a doi: 10.4208/cicp.040413.160114a
[39]	W. Boscheri, M. Tavelli, C. E. Castro, An all Froude high order IMEX scheme for the shallow water equations on unstructured Voronoi meshes, Appl. Numer. Math., 185 (2023), 311–335. https://doi.org/10.1016/j.apnum.2022.11.022 doi: 10.1016/j.apnum.2022.11.022
[40]	G. Huang, Y. Xing, T. Xiong, High order well-balanced asymptotic preserving finite difference WENO schemes for the shallow water equations in all Froude numbers, J. Comput. Phys., 463 (2022), 111255. https://doi.org/10.1016/j.jcp.2022.111255 doi: 10.1016/j.jcp.2022.111255
[41]	S. Vater, R. Klein, A semi-implicit multiscale scheme for shallow water flows at low Froude number, Commun. Appl. Math. Comput. Sci., 13 (2018), 303–336. https://doi.org/10.2140/camcos.2018.13.303 doi: 10.2140/camcos.2018.13.303
[42]	X. Xie, H. Dong, M. Li, High order well-balanced asymptotic preserving IMEX RKDG schemes for the two-dimensional nonlinear shallow water equations, J. Comput. Phys., 510 (2024), 113092. https://doi.org/10.1016/j.jcp.2024.113092 doi: 10.1016/j.jcp.2024.113092
[43]	H. Zakerzadeh, The RS-IMEX scheme for the rotating shallow water equations with the Coriolis force, In: Finite volumes for complex applications Ⅷ—Hyperbolic, Elliptic and Parabolic Problems, Cham: Springer, 200 (2017), 199–207. https://doi.org/10.1007/978-3-319-57394-6-22
[44]	S. Busto, M. Dumbser, A staggered semi-implicit hybrid finite volume / finite element scheme for the shallow water equations at all Froude numbers, Appl. Numer. Math., 175 (2022), 108–132. https://doi.org/10.1016/j.apnum.2022.02.005 doi: 10.1016/j.apnum.2022.02.005
[45]	A. Duran, F. Marche, R. Turpault, C. Berthon, Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes, J. Comput. Phys., 287 (2015), 184–206. https://doi.org/10.1016/j.jcp.2015.02.007 doi: 10.1016/j.jcp.2015.02.007
[46]	A. Kurganov, Y. Liu, M. Lukáčová-Medvi $\check{{\text{d}}}$ ová, A well-balanced asymptotic preserving scheme for the two-dimensional rotating shallow water equations with nonflat bottom topography, SIAM J. Sci. Comput., 44 (2022), A1655–A1680. https://doi.org/10.1137/21M141573X doi: 10.1137/21M141573X
[47]	X. Liu, A. Chertock, A. Kurganov, An asymptotic preserving scheme for the two-dimensional shallow water equations with Coriolis forces, J. Comput. Phys., 391 (2019), 259–279. https://doi.org/10.1016/j.jcp.2019.04.035 doi: 10.1016/j.jcp.2019.04.035
[48]	H. Egger, J. Giesselmann, T. Kunkel, N. Philippi, An asymptotic-preserving discretization scheme for gas transport in pipe networks, IMA J. Numer. Anal., 43 (2023), 2137–2168. https://doi.org/10.1093/imanum/drac032 doi: 10.1093/imanum/drac032
[49]	H. Egger, N. Philippi, An asymptotic preserving hybrid-dG method for convection-diffusion equations on pipe networks, arXiv preprint arXiv: 2209.04238, 2022. https://doi.org/10.48550/arXiv.2209.04238
[50]	A. Kurganov, S. Noelle, G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), 707–740. https://doi.org/10.1137/S1064827500373413 doi: 10.1137/S1064827500373413
[51]	A. Kurganov, E. Tadmor, Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers, Numer. Methods Partial Differ. Equations, 18 (2002), 584–608. https://doi.org/10.1002/num.10025 doi: 10.1002/num.10025
[52]	G. Wanner, E. Hairer, Solving Ordinary Differential Equations II, New York: Springer Berlin Heidelberg, 1996.
[53]	X. Zhong, Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows, J. Comput. Phys., 128 (1996), 19–31. https://doi.org/10.1006/jcph.1996.0193 doi: 10.1006/jcph.1996.0193
[54]	M. K. Banda, M. Herty, A. Klar, Gas flow in pipeline networks, Netw. Heterog. Media, 1 (2006), 41–56. https://doi.org/10.3934/nhm.2006.1.41 doi: 10.3934/nhm.2006.1.41
[55]	H. Egger, J. Giesselmann, Stability and asymptotic analysis for instationary gas transport via relative energy estimates, Numer. Math., 153 (2023), 701–728. https://doi.org/10.1007/s00211-023-01349-9 doi: 10.1007/s00211-023-01349-9
[56]	K. Ehrhardt, M. C. Steinbach, Nonlinear optimization in gas networks, In: Modeling, Simulation and Optimization of Complex Processes, Berlin Heidelberg: Springer, 2005,139–148. https://doi.org/10.1007/3-540-27170-8-11
[57]	M. Herty, N. Izem, M. Seaid, Fast and accurate simulations of shallow water equations in large networks, Comput. Math. Appl., 78 (2019), 2107–2126. https://doi.org/10.1016/j.camwa.2019.03.049 doi: 10.1016/j.camwa.2019.03.049
[58]	J. Caputo, D. Dutykh, B. Gleyse, Coupling conditions for water waves at forks, Symmetry, 11 (2019), 434. https://www.mdpi.com/2073-8994/11/3/434
[59]	R. M. Colombo, M. Garavello, A well posed Riemann problem for the $p$ -system at a junction, Netw. Heterog. Media, 1 (2006), 495–511. https://doi.org/10.3934/nhm.2006.1.495 doi: 10.3934/nhm.2006.1.495
[60]	F. M. White, H. Xue, Fluid Mechanics, New York: McGraw-hill, 2003.
[61]	A. Ambroso, C. Chalons, F. Coquel, E. Godlewski, F. Lagoutière, P. A. Raviart, et al., Coupling of general Lagrangian systems, Math. Comp., 77 (2008), 909–941. https://doi.org/10.1090/S0025-5718-07-02064-9 doi: 10.1090/S0025-5718-07-02064-9
[62]	R. Borsche, Numerical schemes for networks of hyperbolic conservation laws, Appl. Numer. Math., 108 (2016), 157–170. https://doi.org/10.1016/j.apnum.2016.01.006 doi: 10.1016/j.apnum.2016.01.006
[63]	G. Bretti, R. Natalini, B. Piccoli, Fast algorithms for the approximation of a traffic flow model on networks, Discrete Contin. Dyn. Syst. Ser. B, 6 (2006), 427–448. https://doi.org/10.3934/dcdsb.2006.6.427 doi: 10.3934/dcdsb.2006.6.427
[64]	C. Chalons, P. A. Raviart, N. Seguin, The interface coupling of the gas dynamics equations, Quart. Appl. Math., 66 (2008), 659–705. https://doi.org/10.1090/S0033-569X-08-01087-X doi: 10.1090/S0033-569X-08-01087-X
[65]	F. Dubois, P. LeFloch, Boundary conditions for nonlinear hyperbolic systems of conservation laws, J. Differ. Equations, 71 (1988), 93–122. https://doi.org/10.1016/0022-0396(88)90040-X doi: 10.1016/0022-0396(88)90040-X
[66]	E. Godlewski, P. A. Raviart, The numerical interface coupling of nonlinear hyperbolic systems of conservation laws: Ⅰ. The scalar case, Numer. Math., 97 (2004), 81–130. http://link.springer.com/10.1007/s00211-002-0438-5 doi: 10.1007/s00211-002-0438-5
[67]	E. Godlewski, P. A. Raviart, A method of coupling non-linear hyperbolic systems: Examples in CFD and plasma physics, Int. J. Numer. Methods Fluids, 47 (2005), 1035–1041. https://doi.org/10.1002/fld.856 doi: 10.1002/fld.856
[68]	H. Holden, N. H. Risebro, A mathematical model of traffic flow on a network of unidirectional roads, SIAM J. Math. Anal., 26 (1995), 999–1017. https://doi.org/10.1137/S0036141093243289 doi: 10.1137/S0036141093243289
[69]	L. O. Müller, P. J. Blanco, A high order approximation of hyperbolic conservation laws in networks: Application to one-dimensional blood flow, J. Comput. Phys., 300 (2015), 423–437. https://doi.org/10.1016/j.jcp.2015.07.056 doi: 10.1016/j.jcp.2015.07.056
[70]	A. Harten, P. D. Lax, B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev., 25 (1983), 35–61. https://doi.org/10.1137/1025002 doi: 10.1137/1025002
[71]	K. A. Lie, S. Noelle, On the artificial compression method for second-order nonoscillatory central difference schemes for systems of conservation laws, SIAM J. Sci. Comput., 24 (2003), 1157–1174. https://doi.org/10.1137/S1064827501392880 doi: 10.1137/S1064827501392880
[72]	H. Nessyahu, E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), 408–463. https://doi.org/10.1016/0021-9991(90)90260-8 doi: 10.1016/0021-9991(90)90260-8
[73]	P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21 (1984), 995–1011. https://doi.org/10.1137/0721062 doi: 10.1137/0721062
[74]	M. Garavello, K. Han, B. Piccoli, Models for Vehicular Traffic on Networks, American Institute of Mathematical Sciences, 2016.
[75]	O. Kolb, J. Lang, P. Bales, An implicit box scheme for subsonic compressible flow with dissipative source term, Numer. Algor., 53 (2010), 293–307. https://doi.org/10.1007/s11075-009-9287-y doi: 10.1007/s11075-009-9287-y
[76]	H. Egger, A robust conservative mixed finite element method for isentropic compressible flow on pipe networks, SIAM J. Sci. Comput., 40 (2018), A108–A129. https://doi.org/10.1137/16M1094373 doi: 10.1137/16M1094373
[77]	J. G. Caputo, D. Dutykh, Nonlinear waves in networks: Model reduction for the sine-Gordon equation, Phys. Rev. E, 90 (2014), 022912. https://link.aps.org/doi/10.1103/PhysRevE.90.022912 doi: 10.1103/PhysRevE.90.022912
[78]	D. Dutykh, J. G. Caputo, Wave dynamics on networks: Method and application to the sine-Gordon equation, Appl. Numer. Math., 131 (2018), 54–71. https://doi.org/10.1016/j.apnum.2018.03.010 doi: 10.1016/j.apnum.2018.03.010
[79]	R. Borsche, J. Kall, ADER schemes and high order coupling on networks of hyperbolic conservation laws, J. Comput. Phys., 273 (2014), 658–670. https://doi.org/10.1016/j.jcp.2014.05.042 doi: 10.1016/j.jcp.2014.05.042
[80]	M. Herty, N. Kolbe, S. Müller, Central schemes for networked scalar conservation laws, Netw. Heterog. Media, 18 (2023), 310–340. https://doi.org/10.3934/nhm.2023012 doi: 10.3934/nhm.2023012

Reader Comments

Your name:*

Email:*
© 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Networks and Heterogeneous Media

1.2 1.8

Metrics

Article views(101) PDF downloads(19) Cited by(0)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(7) / Tables(3)

Networks and Heterogeneous Media

An asymptotic-preserving scheme for isentropic flow in pipe networks

Related Papers:

Abstract

1. Introduction

2. Dimensional analysis

2.1. Continuous extension to pipe networks

2.2. The low mach/high friction limit

2.3. Numerical Difficulties in Low Mach/High Friction regimes

3. Asymptotic-preserving scheme

3.1. Time discretization

3.2. Spatial discretization of non-stiff flux terms

3.3. Fully discrete method

3.4. Numerical diffusion and stability

3.5. Boundary conditions on pipe networks

3.5.1. Coupling conditions at pipe intersections

3.5.2. Additional boundary condition requirements

3.6. The discrete low mach/high friction limit

3.6.1. The asymptotic-preserving property

3.6.2. Full discrete scheme at the low mach/high friction limit

4. Numerical results

4.1. Example 1—Shocktube problem

4.2. Example 2—Convergence tests

4.3. Example 3—Inlet discontinuity

4.4. Example 4—Comparison with the CU explicit scheme

5. Conclusion

Author contributions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog

Abstract

1. Introduction

2. Dimensional analysis

2.1. Continuous extension to pipe networks

2.2. The low mach/high friction limit

2.3. Numerical Difficulties in Low Mach/High Friction regimes

3. Asymptotic-preserving scheme

3.1. Time discretization

3.2. Spatial discretization of non-stiff flux terms

3.3. Fully discrete method

3.4. Numerical diffusion and stability

3.5. Boundary conditions on pipe networks

3.5.1. Coupling conditions at pipe intersections

3.5.2. Additional boundary condition requirements

3.6. The discrete low mach/high friction limit

3.6.1. The asymptotic-preserving property

3.6.2. Full discrete scheme at the low mach/high friction limit

4. Numerical results

4.1. Example 1—Shocktube problem

4.2. Example 2—Convergence tests

4.3. Example 3—Inlet discontinuity

4.4. Example 4—Comparison with the CU explicit scheme

5. Conclusion

Author contributions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References