Coupling conditions for the transition from supersonic to subsonic fluid states

1.   Introduction

Natural gas pipeline transportation systems are an important part of our infrastructure. Besides transport of gas, also $CO_2$ pipeline transport plays an important role in carbon capture and storage, see for example [18]. A well-established mathematical model for such flow is the p-system (1). The p-system in conservative variables density $\rho > 0$ and flux $q$ is given by

where $t\geq0$ is time, $x$ is space, and $\rho \to p(\rho)$ is the pressure law. The system is stated in Eulerian coordinates and it is a simplified model for gas flow in pipes [21] or water flow in open canals [11]. In the latter case $\rho=h$ is the water height and $p(h)=\frac{g}2h^2$ with $g$ denoting the gravitational acceleration. We are interested in the situation of a single junction connected to pipes depicted in Figure 1. The pipes are modeled as directed arcs and the junction as node. We suppose to have $n$ incoming pipes parameterized by $x \in (-\infty {,0}]$ and labeled by $k \in \delta^- := \{ -n, \dots, -1\}.$ Similarly, we suppose to have $m$ outgoing pipes $k \in \delta^+:=\{ 1,\dots, m\}$ and parameterized by $x\in[0,\infty).$ Along each pipe the evolution of the gas is supposed to fulfill (1) and the state of the system in pipe $k$ is denoted by $U_k:=(\rho_k,q_k)^t.$ The dynamics of the pipes is then coupled at $x=0$ through coupling conditions [2]. General results on the p-system at a junction are given in [9]. For a recent overview of theoretical and numerical results for flows on networks governed by hyperbolic equations we refer to [4]. Except for [14] all current coupling conditions rely on a transversality condition [10] to establish well-posedness. This condition is replaced here by a solvability condition prior to Theorem 2.3 having a similar meaning, namely, that Lax curves are not parallel in the point of intersection in phase space.

It has been shown that the transversality condition is fulfilled for subsonic data, i.e. initial conditions $U^0_i(x)$ with $i \in \delta^+ \cup \delta^-$ fulfilling $\lambda_1(U^0_i) < 0 < \lambda_2(U^0_i)$ where $\lambda_{{j}}, {j}=1,2$ are the characteristic velocities given by equation (5). To the best of our knowledge, no results in the supersonic case exist. However, in the scalar case on networks this case has been well-understood and discussed intensively, see for example [17,6] or [12]. We are interested in the coupling of supersonic and subsonic states for the p-system. We focus on the most common coupling condition for gas networks as well as coupled open canals, namely, we prescribe

for some unknown nodal pressure $p^*(t)$ and where $0\pm$ denote the trace values coming from the interior of the pipe. We refer to [2,15,9] for a more detailed discussion of coupling conditions. Our main result is Theorem 2.3 showing that there exists a solution to the coupling conditions (2) under suitable assumptions on the initial data. The key assumption is on the sign eigenvalues of the initial datum, i.e., equation (17). Briefly summarized the assumption states that the data on incoming arcs is supersonic and the data on outgoing arcs is subsonic.

The study of the coupling of subsonic and supersonic flow is also important for algorithmic purposes. For example in an iterative method for the solution of optimal control problems (see for example [7]), it can happen that during the iteration in an intermediate step, supersonic states appear even if the optimal state is subsonic. The transition from supersonic to subsonic flow is important for high-speed aircraft, where a supersonic external freestream flow must be decelerated efficiently to subsonic speed at the turbine inlet. This transition from supercritical flow to subcritical flow also occurs in Fanno flows, see [22], (p. 199). The transition from supercritical to subcritical regime in free surface flow of yield stress fluids has been studied in [23] where a hydraulic jump during laminar flow of an Herschel-Bukley fluid in a rectangular channel is studied. An experimental assessment of scale effects affecting two-phase flow properties in hydraulic jumps is given in [20].

In this work we will introduce a notion of supply and demand for the state of the system on the adjacent pipes similar to the notion in traffic flow [19,14]. This will allow to solve for constant initial data $U^0_k=cst$ the problem (1) and (2) also in the supersonic case. Necessary conditions on the constants $U^0_k$ are presented ensuring solvability. In the case of $n=m=1$ the solution coincides with a classical solution to a Riemann problem.

2.   Theoretical results

We recall the Lax-curves for the p-system for a pressure law satisfying the following assumption on the equation of state, i.e. $p \in C^3(\mathbb{R}^+;\mathbb{R}^+)$ with

and such that $p(\rho), p'(\rho) \to 0$ for $\rho\to 0.$ Note that the first condition (3a) implies the usual assumption on the equation of state

see e.g. [9,5]. This also implies that both characteristic fields are genuinely nonlinear. In addition, the second condition (3b) is required later in order to obtain strict convex (resp. concave) wave curves. The assumptions (3) are fulfilled by the pressure law of the isentropic Euler equations $p(\rho)=\rho^\gamma$ for ${\gamma > 1}$ as well as for the shallow-water equations where the hydrostatic pressure is $p(\rho)=\frac{g}2 \rho^2.$ In particular, the isothermal equation $p(\rho)=\rho$ is frequently applied in network simulations. We note that a sufficient condition for (3b) is $p'''(\rho)\ge 0$ that for isentropic Euler equations only holds for $\gamma\ge 2$.

Under the assumption (3a) the p-system (1) is strictly hyperbolic in the non-vacuum states $\{ U: \rho>0 \}.$ The speed of sound is given by $c(\rho)=\sqrt{ p'(\rho) }$ and the eigenvalues of the $i$th characteristic field are for $i=1,2$

Both characteristic fields are genuinely nonlinear if the equation of state fulfills (4) and (3a). The 1-Lax-curve through ${\bar{U}=}U_\ell$ is given by

The reversed 2-Lax-curve through ${\bar{U}=}U_r$ is given by

Here, for $i=1,2$ we have the positive and negative branch of

The precise form and derivation of the curves can be found e.g. in [9,13]. Therein, also the following properties are established: The shock curves $S_i, i=1,2,$ are concave and convex with respect to $U$ in the $\rho$-$q$-plane. The rarefaction curves $R_i,i=1,2$, are concave and convex, respectively, in the $\rho$-$q$-plane. The Riemann problem has a global unique solution, see [5,Section 6.3]. The remaining Lax curves are given for example in [9,Equation 4.1,Equation 4.2].

From the explicit form of $R_i$ we obtain

Since $\sigma> \bar{\rho }>0$, the discriminant $s(\sigma) = \frac{\sigma}{\bar{\rho }} ( \sigma-\bar{\rho }) ( p(\sigma)-p(\bar{\rho }))$ in $S_{i}(\sigma; \bar{U})$ is positive. Furthermore, due to assumption (3), $s$ is strict convex, i.e.

Hence,

The function $\tilde{s}(\sigma)= 2 s''(\sigma) s(\sigma) - s'(\sigma)^2$ vanishes at $\sigma=\bar{\rho }$ and $\tilde{s}'(\sigma) = 2 s'''(\sigma) s(\sigma).$

We also note that

provided that condition (3b) holds true. Then, $\tilde{s}$ is a monotone increasing function for all $\sigma>\bar{\rho }$ starting at zero. This yields the strict concavity and convexity of $S_i,$ respectively.

Under assumption (3b) we therefore obtain

Hence, $L_1$ and $L_2^-$ have a single maximum and minimum, respectively. We denote by $q^*=q^*(\bar{U})$ the maximum of the functions $\rho \to L_1(\rho;{\bar{U}})$ and $\rho \to - L_2^{-}(\rho;\bar{U})$. The corresponding density is denoted by $\rho^*= \rho^{*}(\bar{U})$. The {Heaviside} function is denoted by $H(x).$

2.1.   The notion of demand and supply for the p-system

We now define the demand and supply for the p-system extending the notion in the scalar case [19]. In the $\rho$-$q$-plane and for $\rho >0$ the demand $\rho {\,\mapsto\,} d(\rho; \bar{U})$ is an extension of the non-decreasing part of the curve $\rho {\,\mapsto\,} L_1(\rho;U).$ In the $\rho$-$q$-plane and for $\rho >0$ the supply function $\rho {\,\mapsto\,} s(\rho;\bar{U})$ is an extension of the non-increasing part of the curve $\rho {\,\mapsto\,} - L_2^-(\rho;U)$. The maximum of $L_1$ and $-L_2^-$ depend on $\bar{U}.$ For readability we skip this dependence in the following formulas for $d$ and $s.$

A sketch of the supply and demand for different states $\bar{U}=(\bar{\rho },\bar{q})$ is given in Figure 2 and Figure 3, respectively.

The demand and supply are used to determine admissible fluxes for waves of non-positive and non-negative speed, respectively. A Riemann problem is a Cauchy problem for equation (1) and initial data $U_0(x) = \begin{pmatrix} U_\ell&x<0 \\ U_r&x>0\end{pmatrix}.$ A weak solution to the Riemann problem is defined e.g. in [3] and consists of a superposition of self-similar waves. In the case of the $2\times 2$ system we obtain at most two waves separated by an intermediate state. The general structure in the present case therefore is a superposition of either a shock or rarefaction wave of the first family with a shock or rarefaction wave of the second family.

Theorem 2.1. Let $\bar{U}$ with $\bar{\rho }>0$ be given and assume (3). Then, for all $q \leq d(\bar{\rho };\bar{U}),$ there exists a state $\tilde{U}=(\tilde{\rho},q)$ such that the Riemann problem for (1) with $U_\ell=\bar{U}$ and $U_r = \tilde{U}$ admits a weak self-similar solution composed of waves of non-positive speed.

Proof. We have to distinguish two cases. Suppose $\bar{\rho }<\rho^*(\bar{U})$, see Figure 3 (right). Then, $d(\bar{\rho };\bar{U})=\bar{q}.$ Denote by $\tau \not = \bar{\rho }$ the value of $\rho,$ such that $\bar{q} = L_1(\tau;\bar{U}).$ Due to the strict concavity of $L_1$ and for $q< d(\bar{\rho };\bar{U})$ there exists a unique density $\tilde{\rho}>\tau,$ such that $L_1(\tilde{\rho};\bar{U})=q.$ The solution to the Riemann problem consists of a single shock propagating with non-positive velocity. Suppose $\bar{\rho }<\rho^*(\bar{U})$ and assume $q=\bar{q}=d(\bar{\rho };\bar{U}).$ Then, $\tilde{U}=\bar{U}$ as well as $\tilde{U}=( \tau, q)$ where $\tau\not=\bar{\rho }$ fulfills $L_1(\tau;\bar{U})=q.$ Note that in the latter case the solution is a shock of zero speed. Since this is a stationary discontinuity located at a single point we can exclude this solution from our construction as in [16].

In the case $\bar{\rho }>\rho^*(\bar{U})$, see Figure 3 (left), we obtain $d(\bar{\rho };\bar{U})=q^*(\bar{U}).$ For $\bar{q} \leq q \leq q^*(\bar{U})$ we find a density $\tilde{\rho}$ such that $\rho^* \leq \tilde{\rho}{(\bar{U})} \leq \bar{\rho }$ and $q = L_1(\tilde{\rho}; \bar{U}).$ The solution to the corresponding Riemann problem consists of a rarefaction wave with non-positive speed. For $q<\bar{q},$ we find $\tilde{\rho}>\bar{\rho }$ and $q=L_1(\tilde{\rho}; \bar{U}),$ such that the solution to the corresponding Riemann problem consists of a shock of negative speed.

Corollary 1. Let $\bar{U}$ with $\bar{\rho }>0$ be given and assume (3). Then, for all$q \geq s(\bar{\rho };\bar{U}),$ there exists a state $\tilde{U}=(\tilde{\rho},q)$ such that the Riemann problem for (1) with$U_\ell=\tilde{U}$ and $U_r = \bar{U}$ admits a weak self-similar solution composed of waves of non-negative speed.

The proof is similar to the proof of Theorem 2.1 and omitted. It should be noted that in the formulation of Corollary 1 the right state $\bar{U}$ is given instead of the left state in Theorem 2.1.

2.2.   The coupling at a single junction

We introduce the following subsets of $\mathbb{R}^+\times \mathbb{R}$ depending on the sign of the eigenvalues $\lambda_i$ given by (5).

The set $A$ comprises of the states in the subsonic region, $A^+ \cup A^-$ is the supersonic region. For the junction located at $x=0$ we consider the Riemann problem

Depending on the pipe only one of the Riemann data is defined at $t=0,$ i.e. $U_k^- = U_{0,k}$ if $k\in\delta^-$ and $U_k^+ = U_{0,k}$ if $k\in\delta^+$. For each pipe the remaining unknown state $U_k^\pm$ for $k \in \delta^\pm$ is determined in such a way that equation (2) holds true. We construct a Riemann solver at the node $(U_k(t,x))_{k \in \delta^\pm }$ to (13) and (2). An explicit construction is provided for constant initial data self-similar solutions $U_k(t,x)$ consisting of wave of non-positive speed and non-negative speed for $k \in \delta^-$ and $k\in\delta^+,$ respectively. This solution is restricted to $U_k(t,x)$ for $x<0$ and $x>0,$ when $k\in \delta^-$ and $k\in\delta^+,$ respectively.

The demand and supply allows to construct solutions of certain wave speeds. Consider $k \in \delta^-.$ Let $U_k^-:=U_{0,k}$ and $q\leq d(\rho_{0,k}; U_{0,k})$ be given. Then, according to Theorem 2.1, there exist a state $U_k^+$ with $q_k^+=q$ and such that the Riemann problem (13) admits a solution $U_k(t,x)$. The solution consists of a wave of non-positive speed. Similarly, for $k \in\delta^+$, data $U_k^+:=U_{0,k}$ and any $q \geq s(\rho_{0,k}; U_{0,k})$ a self-similar solution $U_k(t,x)$ consisting of waves of non-negative speed exist according to Corollary 1.

In order to obtain a unique density at the node we define the function $\bar{U} {\,\mapsto\,} \tau_k(\bar{U})$ by

Similarly,

Clearly, $\tau_k>0$ if $\bar{\rho }>0$. $\tau_k$ is well-defined due to equation (8). For $k\in \delta^-$ and $\bar{\rho }<\rho^*(\bar{U})$ we obtain $L_1(\tau_k(\bar{U}); \bar{U})= L_1(\bar{\rho };\bar{U}){=\bar{q}}.$ Hence, in this case $\tau_k(\bar{U})$ is the density such that the wave connecting $U_\ell=(\bar{\rho },\bar{q})$ and $U_r=(\tau_k(\bar{U}), \bar{q})$ has zero speed and is omitted as in the proof of Theorem 2.1.

Finally, we turn to the description of the Riemann solver for problem (13) and (2). Due to assumption (3) we have that the condition (2) is equivalent to

and

Note that in the case $n=m=1,$ i.e. $\delta^\pm=\{\pm 1\},$ the conditions (15) and (2) imply by means of the Rankine-Hugoniot jump conditions in particular

For consistency a solution to (13) should therefore be also a solution to a (classical) Riemann problem for (1) with initial data $U(0,x)= \begin{pmatrix} U_{-1}^-&x<0 \\ U_1^+&x>0 \end{pmatrix}.$

We are interested in supersonic and subsonic coupling at the junction. Therefore, we assume in the following that

Remark 1. Due to the symmetry of the Lax curves for the p-system, the discussion also includes the case where $U_k^- \in A$ for $k\in \delta^-$ and $U_k^+ \in A^-$ for $k\in\delta^+.$ The case $U_k^\pm \in A$ for $k\in\delta^\pm$ has been discussed for example in [8].

We now plan to construct a solution fulfilling the coupling conditions. As preliminary result we prove the monotonicity of the Lax curves.

Lemma 2.2. Assume $U_{0,k}$ for $k \in \delta^{\pm}$ with $\rho_{0,k}>0$ are given. Suppose further that (3) holds true.

Then, for all $\sigma \geq \sigma_{lb}$ we have that $\sigma {\,\mapsto\,} L_1(\sigma;U_{0,k})$ for $k\in \delta^{-}$ is strictly monotone decreasing and $\sigma {\,\mapsto\,} L_2^{-}(\sigma;U_{0,k})$ for $k\in\delta^{+}$ is strictly monotone increasing.

Proof. In the case $\rho_{0,k} < \rho^{*}(U_{0,k})$ we have that $\tau_k(U_{0,k})>\rho^{*}(U_{0,k})$ by definition (14). Therefore, the proof follows from the monotonicity properties of $L_1(\cdot)$ and $L_2^{-}(\cdot)$, respectively. The monotonicity implies that $\sigma {\,\mapsto\,} \sum\nolimits_{k\in\delta^{-} } L_1(\sigma; U_{0,k})$ and $\sigma {\,\mapsto\,} \sum\nolimits_{k\in\delta^{+} } L_2^{-}(\sigma; U_{0,k})$ are monotone decreasing and increasing in $[\sigma_{lb},\infty)$, respectively. Furthermore, $L_1(\sigma_{lb};U_{0,k}) \leq d(\sigma_{lb};U_{0,k}), k\in \delta^{-}$ and $L_2(\sigma_{lb};U_{0,k}) \geq s(\sigma_{lb};U_{0,k}), k\in\delta^{+}.$ This finishes the proof.

In general, under assumption (17), the solution to a Riemann problem for (1) consists of a superposition of waves of the first family and waves of the second family. In order to solve the problem (13) we mimic this behavior by constructing of a solution as follows. Denote by

Then, equation (16) is fulfilled provided that there exists $\sigma$ with $\sigma\geq \sigma_{lb}$ and (16)

The following Theorem guarantees the solvability of the previous equation and presents the main result.

Theorem 2.3. Assume $U_{0,k}$ for $k \in \delta^{\pm}$ with $\rho_{0,k}>0$ are given. Suppose further that (3) holds true. Let $\sigma_{lb}$ be given by equation (18).

Provided that

there exists a unique $\sigma\geq \sigma_{lb},$ such that equation (19) is fulfilled.

The proof follows immediately from the monotonicity properties established before.

Now, we continue our construction: Suppose (19) has a unique solution $\overline{\sigma},$ see Theorem 2.3. Then, define

Then, we have $L_1(\overline{\sigma}; U_{0,k}) \leq d( \rho_{0,k}; U_{0,k})$ for $k\in \delta^{-}$ and $L_2^{-}(\overline{\sigma}; U_{0,k}) \geq s( \rho_{0,k}; U_{0,k})$ for $k\in\delta^{+}$ due to Lemma 2.2. According to Theorem 2.1 and Corollary 1 this ensures that the solution to problem (13) has non-positive and non-negative speed, respectively. Therefore, a solution to the problem (13) fulfills equation (16) and (15).

Some remarks are in order.

Remark 2. Note that this construction is also valid in the case $U_k^- \in A$ and coincides with the solution to the equal pressure condition proposed in [8,2].

In the case $n=m=1$ we have $U_{-1}^+ = U_1^-=:\bar{U}.$ Since $\lambda_1(\bar{U})< \lambda_2(\bar{U})$ the solution (13) consists by construction of at most a wave of the first family on $k=-1$. On $k=1$ we may have at most a wave of the second family. This is consistent with the solution to a Riemann problem (13) with initial data $U(0,x)= \begin{pmatrix} U_{-1}^-&x<0 \\ U_1^+&x>0 \end{pmatrix}.$

Remark 3. Note that linearization techniques in order to solve (2) and (16) will fail to provide a solution to (13) due to the fact that at $U_{0,k}, k\in \delta^{-}$ both eigenvalues are non-negative. Therefore, a linear system obtained through linearization at $U_{0,k}, k\in\delta^{\pm}$ would have only $m$ degrees of freedom to prescribe boundary conditions. However, the given set of conditions are in fact $n+m.$ This is different compared with the case $U_{0,k} \in A$ for all $k\in\delta^{\pm}$.

In the case $U_{0,k}\in A^{-}$ for $k\in\delta^{+}$ the problem (13), (15) and (16) is more complicated in the following way. The characteristic families can no longer be uniquely related to whether they are incoming $(k\in\delta^{-})$ or outgoing $(k\in\delta^{+})$. One can easily find examples such that we observe two waves of different families on the outgoing and a wave on the incoming pipe. Therefore, new additional conditions have to be imposed to treat the additional degree of freedom.

2.3.   Example

We consider $p(\rho)=\rho^2$ as an example. Further, we assume $\delta^-=\{-3,\dots, -1\}$ and $\delta^+=\{ 1,\dots, 4\}.$ Initial data $U_{0,k}$ is chosen randomly such that $U_{0,k} \in A^+$ for $k\in\delta^-$ and $U_{0,k}\in A$ for $k\in\delta^+.$ The initial data $U_{0,k}$ does not fulfill the coupling condition (16) and (15) as seen in Table 1. We then compute the states at the junction states at junction (21) according to the discussion of the previous section. Numerically, we solve equation (19) using Newton's method. The results are given in Table 2. As expected the density and the pressure fulfill (2) and (15). Further, we report on the difference between incoming and outgoing fluxes. The equation (19) has been solved up to machine precision.

3.   Summary

While coupling conditions for the subsonic flow in gas pipelines are well-established, in this paper the coupling conditions for the transonic transition from supersonic flow to subsonic flow are studied that have not been thoroughly analyzed before. In this way we fill a gap in the theory, which is important in order to have a complete theory that covers all cases. This can also be important for algorithmic purposes, since during an iterative method for the solution of an optimal control problem, it may happen that supersonic states occur, even if the optimal state is subsonic. The purely subsonic case has been treated explicitly in [1] and in general e.g. in [8]. For a complete theory of all possible cases the fully supersonic case is still at large. Future work may then include a complete statement provided the challenges in the supersonic regime can be dealt with.