An Asymptotic-Preserving Scheme for Isentropic Flow in 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., [7, 29]. 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, [42, 43, 39] 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 $\bm{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:

The isentropic Euler system (2.1) paired with (2.3) and one of (2.4)–(2.5) has been proven a well-posed problem under the condition that the initial data (i) is not in a vacuum state ($\rho^{(k)}(x,t=0)>0\ \forall k$); (ii) has a flow direction that does not change ($u^{(k)}(x,t=0)\geq 0\ \forall k$); and (iii) is under subsonic conditions ($u^{(k)}(x,t=0)\leq c\ \forall k$) [6, 7, 19]. 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., [33, 34, 1, 18, 23, 44, 64, 15, 10, 5].

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

	$\displaystyle\rho$	$\displaystyle=\rho^{(0)}+\varepsilon^{2}\rho^{(2)}+\cdots,$
	$\displaystyle u$	$\displaystyle=u^{(0)}+\varepsilon^{2}u^{(2)}+\cdots.$

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

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

(2.7)

Note that the $\varepsilon^{1}$ terms are skipped since there are no $\mathcal{O}(\varepsilon^{-1})$ 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:

	$\displaystyle\mathcal{O}(\varepsilon^{-2}):$	$\displaystyle{\left[{\left(\rho^{(0)}\right)}^{\gamma}\right]}_{x}=-\frac{C_{\delta}\kappa}{2}\rho^{(0)}u^{(0)}{\left|u^{(0)}\right|},$		(2.8)
	$\displaystyle\mathcal{O}(1):$	$\displaystyle\rho^{(0)}_{t}+{\left(\rho^{(0)}u^{(0)}\right)}_{x}=0,$
		$\displaystyle{\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}$
		$\displaystyle\hskip 56.9055pt=-\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)}.$

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 (i) the condition related to conservation of momentum will remain as it is $\varepsilon$-independent; and (ii) regardless of your selection of coupling conditions (2.4), (2.5), or (2.6), all simplify to requiring equal pressures at the junction in the low Mach/high friction limit. Note that the different limiting equations have already been presented e.g. in [16]. The main aspect here is on their numerical treatment.

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

$${\left\{u\pm\frac{1}{\varepsilon}\sqrt{p^{\prime}(\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

$${\Delta t}_{\rm{ex}}\leq\nu\frac{{\Delta x}}{\displaystyle{\max_{x}}{\left\{{\left|u\right|}+\frac{1}{\varepsilon}\sqrt{p^{\prime}(\rho)}\right\}}}={\mathcal{O}(\varepsilon{\Delta x})},$$

(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)}}$ [37], 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}_{\rm{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 (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 (2.8), and maintains the correct coupling conditions in the $\varepsilon\rightarrow 0$ limit. The scheme proposed in the following section aims to meet these desired properties on the discrete level.

To discretize the split system (3.2)–(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{\bm{F}}(\bm{U})_{x}$ explicitly, and use an implicit approximation for the stiff flux terms $\widehat{\bm{F}}(\bm{U})_{x}$ and the friction source $\bm{S}(\bm{U}).$ Since the source term $\bm{S}(\bm{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., [73, 79] and references therein. To this end, we use the following first-order IMEX time discretization:

	$\displaystyle\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,$		(3.6)
	$\displaystyle\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}.$		(3.7)

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

$$\bm{R}^{n}:=(R^{\rho,n},R^{\rho u,n})^{\top}=-\widetilde{\bm{F}}(\bm{U})_{x}.$$

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

	$\displaystyle\rho^{n+1}$	$\displaystyle=\rho^{n}+{\Delta t}R^{\rho,n}-{\Delta t}(1-\alpha){\left(\rho u\right)}^{n+1}_{x},$		(3.8)
	$\displaystyle{\left(\rho u\right)}^{n+1}$	$\displaystyle=\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]},$		(3.9)

where

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

(3.10)

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

$$\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}.$$

(3.11)

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

We use the Central-Upwind (CU) finite volume scheme to discretize the non-stiff flux terms $\widetilde{\bm{F}}(\bm{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

$$\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}^{n}=\frac{1}{{\Delta x}}\int_{C_{j}}\bm{U}(x,t^{n})\ dx.$$

Then we approximate the contribution of the non-stiff flux terms $\bm{R}:=-\widetilde{\bm{F}}(\bm{U})_{x}$ in each cell using

$${\bm{{\mathcal{R}}}}^{n}_{j}=-\frac{\widetilde{{\bm{{\mathcal{F}}}}}^{n}_{j+\frac{1}{2}}-\widetilde{{\bm{{\mathcal{F}}}}}^{n}_{j-\frac{1}{2}}}{{\Delta x}},$$

(3.12)

where $\widetilde{{\bm{{\mathcal{F}}}}}^{n}_{j\pm\frac{1}{2}}$ are computed using the CU fluxes; see, e.g., [38, 60, 61]:

$$\widetilde{{\bm{{\mathcal{F}}}}}^{n}_{j+\frac{1}{2}}=\frac{s^{+}_{j+\frac{1}{2}}\widetilde{\bm{F}}{\left(\bm{U}^{-}_{j+\frac{1}{2}}\right)}-s^{-}_{j+\frac{1}{2}}\widetilde{\bm{F}}{\left(\bm{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(\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}}\right)}.$$

(3.13)

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

$$\widetilde{\bm{U}}^{n}(x):=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j}+{\left(\bm{U}_{x}\right)}^{n}_{j}{\left(x-x_{j}\right)},\quad x\in C_{j}.$$

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

$$\bm{U}^{-}_{j+\frac{1}{2}}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j}+\frac{{\Delta x}}{2}{\left(\bm{U}_{x}\right)}^{n}_{j},\qquad\bm{U}^{+}_{j+\frac{1}{2}}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j+1}-\frac{{\Delta x}}{2}{\left(\bm{U}_{x}\right)}^{n}_{j+1},$$

(3.14)

where the slopes ${\left(\bm{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., [62, 65, 71]):

$${\left(\bm{U}_{x}\right)}^{n}_{j}={\textrm{minmod}}{\left(\theta\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j+1}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j}}{{\Delta x}},\ \frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j+1}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j-1}}{2{\Delta x}},\ \theta\frac{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{n}_{j-1}}{{\Delta x}}\right)},$$

where the minmod function is

$$\mbox{minmod}(z_{1},z_{2},\ldots)=\left\{\begin{aligned} &\min(z_{1},z_{2},\cdots)&&\mbox{if}\leavevmode\nobreak\ z_{i}>0,\leavevmode\nobreak\ \forall i,\\ &\max(z_{1},z_{2},\cdots)&&\mbox{if}\leavevmode\nobreak\ z_{i}<0,\leavevmode\nobreak\ \forall i,\\ &0&&\mbox{otherwise},\end{aligned}\right.$$

and is applied component-wise. Lastly, $s^{\pm}_{j+\frac{1}{2}}$ of (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{\bm{F}}/\partial\bm{U}$ seen in (3.4):

	$\displaystyle s^{+}_{j+\frac{1}{2}}=$	$\displaystyle\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^{\prime}{\left(\rho^{-}_{j+\frac{1}{2}}\right)}-a^{n}\right]}},\right.$		(3.15)
		$\displaystyle\hskip 56.9055pt\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^{\prime}{\left(\rho^{+}_{j+\frac{1}{2}}\right)}-a^{n}\right]}},0\right\},$
	$\displaystyle s^{-}_{j+\frac{1}{2}}=$	$\displaystyle\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^{\prime}{\left(\rho^{-}_{j+\frac{1}{2}}\right)}-a^{n}\right]}},\right.$
		$\displaystyle\hskip 56.9055pt\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^{\prime}{\left(\rho^{+}_{j+\frac{1}{2}}\right)}-a^{n}\right]}},0\right\}.$

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 (i) compute non-stiff flux terms $\bm{R}^{n}=-\widetilde{\bm{F}}(\bm{U})_{x}$ using the CU numerical fluxes described in §3.2; (ii) discretize the first spatial derivative terms using standard second-order central difference; and (iii) 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 $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j}^{n+1}$ reads

		$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n+1}_{j}-\frac{{\Delta t}^{2}}{{\Delta x}^{2}\varepsilon^{2}}a^{n}(1-\alpha){\left[\phi^{n}_{j+\frac{1}{2}}\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n+1}_{j+1}-{\left(\phi^{n}_{j+\frac{1}{2}}+\phi^{n}_{j-\frac{1}{2}}\right)}\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n+1}_{j}+\phi^{n}_{j-\frac{1}{2}}\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n+1}_{j-1}\right]}$		(3.16)
		$\displaystyle\hskip 113.81102pt=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n}_{j}+{\Delta t}{\mathcal{R}}^{\rho,n}_{j}-{\Delta t}(1-\alpha)D_{x}{\left[\frac{1}{\Psi^{n}_{j}}{\left({\left(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}}\right)}^{n}_{j}+{\Delta t}{\mathcal{R}}^{\rho u,n}_{j}\right)}\right]},$

where ${\bm{{\mathcal{R}}}}^{n}_{j}=({\mathcal{R}}^{\rho,n}_{j},{\mathcal{R}}^{\rho u,n}_{j})^{\top}$ is defined in (3.12)–(3.13), $\Psi_{j}^{n}=\Psi^{n}(x_{j})$ from (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 (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 right hand side of (3.16) is non-singular (by Gershgorin Theorem).

After solving the linear system in (3.16) for ${\left\{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j}^{n+1}\right\}}$, it can be directly substituted into the spatially discretized version of (3.9) to obtain the fully discrete approximation for ${\left(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}}\right)}^{n+1}_{j}$:

$${\left(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}}\right)}^{n+1}_{j}=\frac{1}{\Psi^{n}_{j}}{\left[{\left(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}}\right)}^{n}_{j}+{\Delta t}{\mathcal{R}}^{\rho u,n}_{j}-\frac{a^{n}{\Delta t}}{\varepsilon^{2}}D_{x}\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n+1}_{j}\right]}.$$

(3.17)

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 [37], 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 (3.6)–(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 (3.4), resulting in the following time-step restriction:

$${\Delta t}_{\rm{AP}}\leq\nu\frac{{\Delta x}}{\displaystyle{\max_{x}}{\left\{{\left|u\right|}+\sqrt{{\left(1-\alpha\right)}u^{2}+\frac{\alpha}{\varepsilon^{2}}{\left[p^{\prime}(\rho)-a(t)\right]}}\right\}}},$$

(3.18)

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

To ensure ${\Delta t}_{\rm{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 equations (3.16)–(3.17). We start by rewriting the CU flux in (3.13) in the form

$$\widetilde{{\bm{{\mathcal{F}}}}}^{n}_{j+\frac{1}{2}}=\frac{\widetilde{{\bm{F}}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}^{n})+\widetilde{{\bm{F}}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1}^{n})}{2}+{\bm{{\mathcal{D}}}}^{n}_{j+\frac{1}{2}},$$

where

	$\displaystyle{\bm{{\mathcal{D}}}}^{n}_{j+\frac{1}{2}}=$	$\displaystyle\frac{s^{+}_{j+\frac{1}{2}}}{2{\left(s^{+}_{j+\frac{1}{2}}-s^{-}_{j+\frac{1}{2}}\right)}}{\left[2\widetilde{{\bm{F}}}(\bm{U}^{-}_{j+\frac{1}{2}})-\widetilde{{\bm{F}}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}^{n})-\widetilde{{\bm{F}}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1}^{n})\right]}$		(3.19)
		$\displaystyle-\frac{s^{-}_{j+\frac{1}{2}}}{2{\left(s^{+}_{j+\frac{1}{2}}-s^{-}_{j+\frac{1}{2}}\right)}}{\left[2\widetilde{{\bm{F}}}(\bm{U}^{+}_{j+\frac{1}{2}})-\widetilde{{\bm{F}}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}^{n})-\widetilde{{\bm{F}}}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{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(\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}}\right)},$

is the numerical diffusion term of the CU discretization, with ${\bm{{\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 (3.15) and $\bm{U}^{\pm}_{j+\frac{1}{2}}$ are defined in (3.14). To compute the leading order of the numerical diffusion in (3.19), we introduce the following asymptotic expansions:

	$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{n}_{j}$	$\displaystyle=\rho^{(0),n}_{j}+\varepsilon^{2}\rho^{(1),n}_{j}+\cdots,$
	$\displaystyle\rho^{\pm}_{j+\frac{1}{2}}$	$\displaystyle=\rho^{(0),\pm}_{j+\frac{1}{2}}+\varepsilon^{2}\rho^{(1),\pm}_{j+\frac{1}{2}}+\cdots.$

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 equation (2.7). In addition, we expand $a^{n}$ in a way that follows its definition in (3.5):

$$a^{n}=\min_{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$:

	$\displaystyle{\mathcal{D}}^{\rho,n}_{j+\frac{1}{2}}=$	$\displaystyle\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]}$		(3.20)
		$\displaystyle-\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]}$
		$\displaystyle+\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)}.$

Since $\alpha=\varepsilon^{b},\ b\geq 2$ from (3.5), this is ${\mathcal{O}}(1)$ by the last term of (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{\bm{F}}(\bm{U})$ seen in equation (3.2). However, since $a^{n}$ and its expansion are related to $p^{\prime}(\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 (3.12), we can obtain the following approximations to ${\bm{{\mathcal{R}}}}_{j}^{n}$:

	$\displaystyle{\mathcal{R}}_{j}^{\rho,n}$	$\displaystyle=-D_{x}\widetilde{F}^{\rho}(\bm{U}_{j}^{n})+{\mathcal{O}}({\Delta x}^{2}),$		(3.21)
	$\displaystyle{\mathcal{R}}_{j}^{\rho u,n}$	$\displaystyle=-D_{x}\widetilde{F}^{\rho u}(\bm{U}_{j}^{n})+{\mathcal{O}}{\left(\frac{{\Delta x}^{2}}{\varepsilon^{2}}\right)},$		(3.22)

where $\widetilde{\bm{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 (3.21) and the ${\mathcal{O}}(\varepsilon^{-2}{\Delta x}^{2})$ term in (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 (3.22) into the fully discrete AP scheme seen (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 (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 (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.