Modeling and simulation of gas networks coupled to power grids

TL;DR

This paper develops a mathematical framework for coupling gas networks with electric grids, modeling the transition from gas to power, and demonstrates how this coupling influences network dynamics through simulations.

Contribution

It introduces a novel coupling model using pressure laws that ensure well-posedness and provides numerical methods for simulating gas-power network interactions.

Findings

Gas-to-power transition affects network dynamics

Coupling model ensures mathematical well-posedness

Simulations illustrate impact on gas and power flows

Abstract

In this paper, a mathematical framework for the coupling of gas networks to electric grids is presented to describe in particular the transition from gas to power. The dynamics of the gas flow are given by the isentropic Euler equations, while the power flow equations are used to model the power grid. We derive pressure laws for the gas flow that allow for the well-posedness of the coupling and a rigorous treatment of solutions. For simulation purposes, we apply appropriate numerical methods and show in a experimental study how gas-to-power might influence the dynamics of the gas and power network, respectively.

Tables

Table 1. Table 1 : (r)arefaction waves and (s)hocks for different values of the density ρ 𝜌 \rho .

value of $\rho$					left wave	right wave
		$\rho$	$\le$	${\rho }_{\text{min}}$	invalid
${\rho }_{\text{min}}$		$\rho$	$\le$	${\rho }_l$	r	r
${\rho }_l$	$\le$	$\rho$	$\le$	${\rho }_r$	s	r
${\rho }_r$	$\le$	$\rho$			s	s

Table 2. (a) wave types for ρ l ≤ ρ r subscript 𝜌 𝑙 subscript 𝜌 𝑟 \rho_{l}\leq\rho_{r} .

value of $\rho$					left wave	right wave
		$\rho$	$\le$	${\rho }_{\text{min}}$	invalid
${\rho }_{\text{min}}$		$\rho$	$\le$	${\rho }_l$	r	r
${\rho }_l$	$\le$	$\rho$	$\le$	${\rho }_r$	s	r
${\rho }_r$	$\le$	$\rho$			s	s

Table 3. (b) wave types for ρ r ≤ ρ l subscript 𝜌 𝑟 subscript 𝜌 𝑙 \rho_{r}\leq\rho_{l} .

value of $\rho$					left wave	right wave
		$\rho$	$\le$	${\rho }_{\text{min}}$	invalid
${\rho }_{\text{min}}$		$\rho$	$\le$	${\rho }_r$	r	r
${\rho }_r$	$\le$	$\rho$	$\le$	${\rho }_l$	r	s
${\rho }_l$	$\le$	$\rho$			s	s

Table 4. Table 2 : Solution structure for ρ l ≤ ρ r subscript 𝜌 𝑙 subscript 𝜌 𝑟 \rho_{l}\leq\rho_{r} .

value of $\epsilon$					left wave	right wave
		$\epsilon$	$\ge$	$\left(L_l-L_r\right)\left({\rho }_{\text{min}}\right)$	invalid
$\left(L_l-L_r\right)\left({\rho }_{\text{min}}\right)$	$>$	$\epsilon$	$\ge$	$\left(L_l-L_r\right)\left({\rho }_l\right)$	r	r
$\left(L_l-L_r\right)\left({\rho }_l\right)$	$\ge$	$\epsilon$	$\ge$	$\left(L_l-L_r\right)\left({\rho }_r\right)$	s	r
$\left(L_l-L_r\right)\left({\rho }_r\right)$	$\ge$	$\epsilon$			s	s

Table 5. Table 3 : Parameters of the gas network.

Pipe	From	To	Length [km]	Diameter [mm]	Roughness [mm]
P10	S4	S20	20.322	600	0.05
P20	S5	S17	20.635	600	0.05
P21	S17	S4	10.586	600	0.05
P22	S17	S8	10.452	600	0.05
P24	S8	S20	19.303	600	0.05
P25	S20	S25	66.037	600	0.05
P99	S4	S8	5.000	600	0.05

Table 6. Table 4 : Parameters of the power grid.

Node	$G$	$B$
N1	0.0000	-17.3611
N2	0.0000	-16.0000
N3	0.0000	-17.0648
N4	3.3074	-39.3089
N5	3.2242	-15.8409
N6	2.4371	-32.1539
N7	2.7722	-23.3032
N8	2.8047	-35.4456
N9	2.5528	-17.3382

Table 7. (a) Busses

Node	$G$	$B$
N1	0.0000	-17.3611
N2	0.0000	-16.0000
N3	0.0000	-17.0648
N4	3.3074	-39.3089
N5	3.2242	-15.8409
N6	2.4371	-32.1539
N7	2.7722	-23.3032
N8	2.8047	-35.4456
N9	2.5528	-17.3382

Table 8. (b) Transmission lines

Edge	From	To	$G$	$B$
TL14	N1	N4	0.0000	17.3611
TL45	N4	N5	-1.9422	10.5107
TL56	N5	N6	-1.2820	5.5882
TL36	N3	N6	0.0000	17.0648
TL67	N6	N7	-1.1551	9.7843
TL78	N7	N8	-1.6171	13.6980
TL82	N8	N2	0.0000	16.0000
TL89	N8	N9	-1.1876	5.9751
TL94	N9	N4	-1.3652	11.6041

Table 9. Table 5 : Initial boundary conditions of the power grid.

Node	$P$	$Q$	$\left|V\right|$	$\varphi$
N1	-	-	1	0
N2	163	-	1	-
N3	85	-	1	-
N4	0	0	-	-
N5	-90	-30	-	-
N6	0	0	-	-
N7	-100	-35	-	-
N8	0	0	-	-
N9	-125	-50	-	-

Equations

$${\left(\begin{array}{c}\rho \\ q\end{array}\right)}_t+{\left(\begin{array}{c}q\\ p(\rho )+\frac{q^2}{\rho }\end{array}\right)}_x=\left(\begin{array}{c}0\\ S(\rho ,q)\end{array}\right),$$

$$p\in C^1({\mathrm{IR}}^{+},{\mathrm{IR}}^{+}),p^{\prime }(\rho )>0\text{ for all }\rho \in {\mathrm{IR}}^{+},$$

$$p(\rho )=c^2\rho ,$$

$$p(\rho )=\kappa {\rho }^{\gamma },$$

$${\lambda }_1(\rho ,q)=\frac{q}{\rho }-\sqrt{p^{\prime }(\rho )},{\lambda }_2(\rho ,q)=\frac{q}{\rho }+\sqrt{p^{\prime }(\rho )},$$

$$r_1(\rho ,q)=\left(\begin{array}{c}-1\\ -{\lambda }_1(\rho ,q)\end{array}\right),r_2(\rho ,q)=\left(\begin{array}{c}1\\ {\lambda }_2(\rho ,q)\end{array}\right).$$

$$p_e=p_{\text{coupling}}$$

$${\rho }_e={\rho }_{\text{coupling}},$$

$$\sum _{\text{incoming pipes}}{q}_{\text{pipe}}=\sum _{\text{outgoing pipes}}{q}_{\text{pipe}}$$

$$\left|\frac{q}{\rho }\right|<c(\rho )=\sqrt{p^{\prime }(\rho )}.$$

$$L_l(\rho ;{\rho }_l,q_l)=\{\begin{array}{ll}\rho \left(\frac{q_l}{{\rho }_l}+{\int }_{\rho }^{{\rho }_l}\frac{c(s)}{s}ds\right)& \text{ for }\rho \le {\rho }_l\text{ (rarefaction)}\\ \rho \frac{q_l}{{\rho }_l}-\sqrt{f(\rho ,{\rho }_l)}& \text{ for }{\rho }_l\le \rho \text{ (shock)},\end{array}$$

$$f(\rho ,{\rho }_l)=a(\rho ,{\rho }_l)\Delta p(\rho ,{\rho }_l)\text{mbox}witha(\rho ,{\rho }_l)=\frac{\rho }{{\rho }_l}(\rho -{\rho }_l),\Delta p(\rho ,{\rho }_l)=p(\rho )-p({\rho }_l)$$

$$L_r(\rho ;{\rho }_r,q_r)=\{\begin{array}{ll}\rho \left(\frac{q_r}{{\rho }_r}-{\int }_{\rho }^{{\rho }_r}\frac{c(s)}{s}ds\right)& \text{ for }\rho \le {\rho }_r\text{ (rarefaction)}\\ \rho \frac{q_r}{{\rho }_r}+\sqrt{f(\rho ,{\rho }_r)}& \text{ for }{\rho }_r\le \rho \text{ (shock)}.\end{array}$$

$$L_l(\rho )-L_r(\rho )=0.$$

$$-L_r(\rho ;{\rho }_r,q_r)=L_l(\rho ;{\rho }_r,-q_r).$$

$$p^{\prime }(\rho )=\alpha {\rho }^{\delta }.$$

$$p(\rho )=\frac{\alpha }{\gamma }{\rho }^{\gamma }+\text{const}$$

$$p(\rho )=\alpha \log (\rho )+\text{const}.$$

$$L_l^{\prime }(\rho )$$

$$L_l^{\prime \prime }(\rho )$$

$$L_r^{\prime }(\rho )$$

$$L_r^{\prime \prime }(\rho )$$

$$2p^{\prime }(\rho )+\rho p^{\prime \prime }(\rho )$$

$$\Delta p^2+a^2\left(2\Delta p{p}^{\prime \prime }-(p^{\prime }{)}^2\right)+\frac{1}{2}(a^2{)}^{\prime }(\Delta p^2{)}^{\prime }$$

$$p_{\infty }-\rho p^{\prime }(\rho )-p(\rho )\ge 0.$$

$$2f{f}^{\prime \prime }-(f^{\prime }{)}^2=\Delta p^2+a^2\left(2\Delta p{p}^{\prime \prime }-(p^{\prime }{)}^2\right)+\frac{1}{2}(a^2{)}^{\prime }(\Delta p^2{)}^{\prime }.$$

$$0<l=\underset{\rho \to 0}{\lim }{L}_l(\rho )=\underset{\rho \to 0}{\lim }\rho {\int }_{\rho }^{{\rho }_l}\frac{c(s)}{s}ds=\underset{\rho \to 0}{\lim }\rho F(\rho ).$$

$${\left({\rho }^2\frac{c(\rho )}{\rho }\right)}^{\prime }=\frac{1}{2\sqrt{p^{\prime }(\rho )}}(2p^{\prime }(\rho )+\rho p^{\prime \prime }(\rho ))\ge 0$$

$$\frac{c(\rho )}{\rho }=-F^{\prime }(\rho )=\frac{l}{{\rho }^2}+\underset{\rho \to 0}{o}\left(\frac{1}{{\rho }^2}\right)$$

$$p^{\prime }(\rho )=c(\rho {)}^2=\frac{l^2}{{\rho }^2}+\underset{\rho \to 0}{o}\left(\frac{1}{{\rho }^2}\right),$$

$$p(\rho )=-\frac{l^2}{\rho }+\underset{\rho \to 0}{o}\left(\frac{1}{\rho }\right).$$

$$\underset{\rho \to 0}{\lim }{L}_l(\rho )=0\iff p\in \underset{\rho \to 0}{o}\left(\frac{1}{\rho }\right).$$

Full text

AMS Classification: 35L65, 65M08

Keywords: Gas networks, pressure laws, power flow equations, simulation

1 Introduction

Power generation is changing drastically with more and more available energy coming from renewable sources. Most of these, like wind and solar power, are much more volatile than traditional sources. The problems of incorporating them into the power system have sparked much research. Due to their flexibility in comparison to other power plants, gas turbines and gas engine power plants are often proposed as a means to balance power demands that cannot be met with renewable power sources at a given time. The repercussions on the gas pipeline networks by this balancing have been studied [6] and also joint optimal control of gas and power networks are investigated [24, 27]. But so far only steady-state flow of the gas dynamics is examined. Although this is often sufficient, the comparably low acoustic speed allows for short term effects that cannot be seen with steady state methods.

We propose a fully time-dependent gas model on a pipeline network coupled to a power network. Gas networks have been studied during the last two decades and questions of node conditions have been analyzed from a theoretical point of view, see for example [1, 2, 4, 5, 7]. Therein, a system of conservation or balance laws based on the Euler equations is used and the concept of Lax curves [22] allows for a discussion on well-posedness of solutions. Finite volume or finite difference methods are typically applied to solve the network problem numerically [15, 19, 21]. Power grids are mainly described by the power flow equations consisting of a nonlinear system of equations to determine power and voltage, see for example [3, 17].

Although the power network could also be modeled with full time-dependency [9], the high signal speed of electricity renders this unnecessary. Coupling of the networks happens at nodes common to both, where a gas-fired generator is placed such that it can convert gas to electricity, see [16, 27]. The coupling is modeled by an equality constraint relating generated power to consumed gas flow. In section 2 and 3, we derive conditions on the gas pressure and power demand, respectively, that guarantee that the coupling conditions lead to well-posedness of the gas-power coupling. For the numerical study in section 4, we first give some validation results for the proposed discretizations, then we present the influence of different pressure laws and lastly we showcase the solution of a coupled gas-power system.

2 Gas networks

Gas networks have been investigated very intensively during the last ten years, see for example [2, 5, 21]. Coupling conditions at nodes have been established to ensure well-posedness of the network solution [7] and a rigorous numerical treatment [8, 15, 21, 25]. As underlying structure for gas networks we consider a directed graph with nodes $V_{\text{Gas}}$ and arcs $E_{\text{Gas}}$. Nodes as well as arcs may have state variables, which are in general time-dependent. We briefly recall the common setting.

In each pipe $e\in E_{\text{Gas}}$, the gas flow is described by the isentropic Euler equations in one space dimension:

[TABLE]

where $t\in\mathrm{IR}^{+}$ is the time, $x\in[0,l_{e}]$ is the position along the pipeline of length $l_{e}$, $\rho$ is the density of the gas, $q$ is the momentum of the gas, $p$ is the pressure and $S$ includes source terms that influence the momentum of the gas. The pressure law must be a function of the form

[TABLE]

to ensure strict hyperbolicity. Usual examples are the isothermal pressure law

[TABLE]

where $c$ is the acoustic speed of the gas, or, more general, the $\gamma$-law

[TABLE]

for suitable constants $\kappa$ and $\gamma$, which we examine in Proposition 2. In fact, we show that the class of possible pressure laws can be enlarged, leading to non-standard pressure functions. From literature we know that the eigenvalues $\lambda$ and the eigenvectors $r$ for the system (1) are given by

[TABLE]

At the (gas) network inflow boundaries, typically the pressure $p_{\text{in}}$ is prescribed and at the outflow boundaries the flow $q_{\text{out}}$ is given. Further, the following coupling conditions, an introduction to which can be found in [2], are used at all inner nodes of $V_{\text{Gas}}$:

• •

equality of pressure: the pressure values at the ends of all arcs connected to the same node must be equal, that is, there is a coupling pressure $p_{\text{coupling}}$ such that

[TABLE]

at the end of all arcs $e$ connected to the junction. This condition is equivalent to a condition

[TABLE]

because $p$ as a strictly increasing function is one-to-one.

• •

conservation of mass: the sum of all incoming fluxes must equal the sum of all outgoing fluxes (including eventual source/sink terms), i.e.,

[TABLE]

2.1 Well-posedness of the Riemann problem for the isentropic Euler equation

In this section, we investigate conditions for the pressure function (3) and the flows that guarantee well-posedness of the Riemann problem for the isentropic Euler equations without a source term, see [4] for an introduction. By the front-tracking technique this is sufficient for the well-posedness of the isentropic Euler equations, see [12]. Assuming strict hyperbolicity, we start with the assumption $p^{\prime}>0$. We also assume that the initial data $U_{l}=(\rho_{l},q_{l})$ and $U_{r}=(\rho_{r},q_{r})$ of the Riemann problem is sub-sonic, i.e.,

[TABLE]

A solution to the Riemann problem with left state $U_{l}$ and right state $U_{r}$ is given by an intersection point of the Lax curves through these two points, again see [4]. The Lax curve of all states $U=(\rho,q)$ reachable via $1$-rarefaction or $1$-shock (that is, waves corresponding to $\lambda_{1}$ above) from a left state $U_{l}=(\rho_{l},q_{l})$ is given by (taken from [7])

[TABLE]

where $f$ is defined by

[TABLE]

such that $U=(\rho,L_{l}(\rho))$ holds.

Analogously, the Lax curve of all states $U=(\rho,q)$ reachable via $2$-rarefaction or $2$-shock from a right state $U_{r}=(\rho_{r},q_{r})$ is given by

[TABLE]

Whenever we write $L_{l}(\rho)$ and $L_{r}(\rho)$, it will mean $L_{l}(\rho,\rho_{l},q_{l})$ and $L_{r}(\rho,\rho_{r},q_{r})$, respectively. Beware that later in section 2.2 different second and third arguments will appear and we will write them out again.

A solution is then at a point $\rho$, where

[TABLE]

An example of this is shown in figure 1.

Note that

[TABLE]

So $L_{l}$ and $-L_{r}$ share most properties. Therefore equation (10) is best understood not as a difference of two functions but as a sum of two very similar functions. Similar to [16], we show certain properties of the Lax curves that are used to define the new pressure laws. There will appear three propositions, which are numbered by A, B and C. These propositions present lists of conditions that are closely related. To make things tractable, each list bears the same numbering but has the letter of the corresponding proposition in front. For example, the conditions A1, B1 and C1 all govern concavity of the Lax curves.

Proposition A** (Proposition 1).**

Let $L\in C^{2}(\mathbb{R}^{+},\mathbb{R}^{+})$ and let it fulfill conditions A1 through A3 (where only one of A3A$3$(a) and A3A$3$(b) must hold).

1. (A1)

$L^{\prime\prime}\leq 0$, 2. (A2)

$\exists\rho>0\mathrel{\mathop{\ordinarycolon}}L(\rho)<0$, 3. (A3)

1. (a)

$0<\lim_{\rho\to 0}L(\rho)<\infty$* or* 2. (b)

$\lim_{\rho\to 0}L(\rho)=0$* and $\lim_{\rho\to 0}L^{\prime}(\rho)>0$.*

Then there is a unique $\rho>0$ with $L(\rho)=0$ such that for all $\hat{\rho}>\rho$ holds $L(\hat{\rho})<0$ and $L(\rho)\to-\infty$ for $\rho\to\infty$.

Proof.

Condition A3 ensures positive values near [math], condition A2 ensures negative values for some $\rho>0$, together yielding a zero in between and condition A1 makes $L$ concave, guaranteeing uniqueness of the zero and implying then that $L(\rho)\to-\infty$ for $\rho\to\infty$. ∎

Note if two functions $L_{1}$, $L_{2}$ satisfy the prerequisites of Proposition A, so does $L_{1}+L_{2}$. Therefore we search for conditions on pressure functions that make $L_{l}$ and $-L_{r}$ satisfy the conditions of Proposition A, as then their sum in equation (10) will, too, and so it will have a unique zero.

Our main result is the following. It generalizes the usual notion of $\gamma$-laws to allow for $\gamma<1$. Further, as we will show below (Proposition C), positive linear combinations of such valid pressure laws again lead to valid pressure laws.

Proposition 2** (generalized $\gamma$-laws).**

Let the derivative of the pressure be given by

[TABLE]

This translates to pressure functions given by

[TABLE]

for $\gamma\neq 0$ and

[TABLE]

For these pressure functions there holds

• •

The Lax curves through any sub-sonic initial states $U_{l}$, $U_{r}$ have a unique intersection point at some $\rho>0$ if and only if $\mathinner{\!\left\lvert\delta\right\rvert}\leq 2$ and $\alpha>0$.

• •

For every $\delta$ with $\mathinner{\!\left\lvert\delta\right\rvert}>2$ there are sub-sonic states $U_{l}$, $U_{r}$ such that the Lax curves have no intersection at all.

Expressed in the usual form of $\gamma$-laws, this result means $p(\rho)=\kappa\rho^{\gamma}$ is a valid pressure function if and only if $0<\gamma<3$ and $\kappa>0$ or $-1<\gamma<0$ and $\kappa<0$. The proof of Proposition 2 will be given at the end of this section. To prove Proposition 2 we will reformulate conditions A1 through A3 in Proposition A in terms of the pressure function.

Because of equation (11), both Lax curves behave essentially identically and therefore we prove our findings for $L_{l}$ only as the proofs for $-L_{r}$ are the same. Along the way the derivatives of the Lax curves will be important and so we provide them here,

[TABLE]

and

[TABLE]

Next, we intend to give conditions under which Proposition A is applicable to $L_{l}$. Before we do so, we provide the following lemma:

Lemma 3**.**

Let $g\in C^{1}(\mathbb{R}^{+},\mathbb{R}^{+})$ be a non-negative function, $g\geq 0$, such that also $(\rho^{2}g(\rho))^{\prime}\geq 0$ for all $\rho>0$ and let $G$ be given by $G(\rho)=\int_{\rho}^{\rho_{l}}g(s)\operatorname{d\!}s$. Then, there holds $\rho^{2}g(\rho)\overset{\rho\to 0}{\longrightarrow}0$ if and only if $\rho G(\rho)\overset{\rho\to 0}{\longrightarrow}0$.

Proof.

The proof can be found in Appendix A. ∎

The role of $g$ in Lemma 3 will be played by $p^{\prime}(\rho)$ and $\frac{c(\rho)}{\rho}$ in the following. Let us now put the focus on the Lax curve $L_{l}$ again. The following Proposition B will turn conditions A1 through A3 in Proposition A into conditions on the pressure function.

Proposition B** (Proposition 4).**

Let $p\in C^{2}(\mathbb{R}^{+})$ with $p^{\prime}>0$. Then conditions A1, A2, A3 hold for $L_{l}$ for all $\rho_{l}$, $q_{l}$ with $\mathinner{\!\left\lvert\frac{q_{l}}{\rho_{l}}\right\rvert}<\sqrt{p^{\prime}(\rho_{l})}$ if and only if conditions B1 through B3 hold (again with only one of B3B$3$(a) and B3B$3$(b) fulfilled).

1. (B1)

These inequalities hold:

[TABLE]

where the arguments $\rho$ and $\rho_{l}$ have been omitted for readability, see equation (8). 2. (B2)

Let $p_{\infty}=\lim_{\rho\to\infty}p(\rho)\in\mathbb{R}\cup\{\infty\}$. For all $\rho>0$ there holds

[TABLE] 3. (B3)

1. (a)

There is $p_{0}>0$ such that $p(\rho)=-\frac{p_{0}}{\rho}+\underset{\rho\to 0}{o}\left(\frac{1}{\rho}\right)$ or 2. (b)

$p(\rho)\in\underset{\rho\to 0}{o}\left(\frac{1}{\rho}\right)$* and $\lim_{\rho\to 0}\int_{\rho}^{\rho_{l}}\frac{c(s)}{s}\operatorname{d\!}s-c(\rho)-c(\rho_{l})\geq 0$ for all $\rho_{l}>0$.*

Proof.

• B1 $\Leftrightarrow$ A1:

The equivalence B1 $\Leftrightarrow$ A1 is immediate from the definitions, inequality (14) is for the rarefaction part, inequality (15) for the shock part. Note that

[TABLE]

• A1, A3A$3$(a)$\Rightarrow$ B1,B3B$3$(a):

Assume conditions A1 and A3A$3$(a). Let $F(\rho)=\int_{\rho}^{\rho_{l}}\frac{c(s)}{s}\operatorname{d\!}s$.

[TABLE]

Therefore $F(\rho)=\frac{l}{\rho}+\underset{\rho\to 0}{o}\left(\frac{1}{\rho}\right)$. Note now that

[TABLE]

due to condition A1. Lemma 3 therefore shows that

[TABLE]

and hence

[TABLE]

which again with lemma 3 yields

[TABLE]

• A1, A3A$3$(a)$\Leftarrow$ B1,B3B$3$(a):

Assume now conditions B1 and B3B$3$(a). We now note that $\left(\rho^{2}p^{\prime}(\rho)\right)^{\prime}=\rho\left(2p^{\prime}(\rho)+\rho p^{\prime\prime}(\rho)\right)\geq 0$, because of B1 and use the lemmas to arrive at A3A$3$(a).

• A1, A3A$3$(b)$\Leftarrow$ B1,B3B$3$(b):

We assume condition B3B$3$(b). The last proof (here we need again conditions A1 and B1 respectively) also shows that

[TABLE]

For the derivative we find

[TABLE]

due to the sub-sonic condition. This is non-negative due to condition B3B$3$(b).

• A1, A3A$3$(b)$\Rightarrow$ B1,B3B$3$(b):

Assume that B3B$3$(b) does not hold, that is

[TABLE]

Choosing $\rho_{l},q_{l}$ such that $\frac{q_{l}}{\rho_{l}}=-c(\rho_{l})+\frac{1}{2}\delta$ yields

[TABLE]

• A2$\Leftarrow$B2:

First of all we note that the limit exists, as $p^{\prime}>0$. For $\rho>\rho_{l}$ we find

[TABLE]

This is less than zero if and only if the bracket is less then zero. If now $p_{\infty}=\infty$, this is fulfilled for some $\rho>\rho_{l}$. If $p_{\infty}<\infty$, we have the limit

[TABLE]

Let now $\delta$ and $\epsilon_{\rho}$ be defined by

[TABLE]

Then $\delta>0$ as the state is sub-sonic and $\epsilon_{\rho}\to 0$. Therefore let $\rho$ be so great that $\delta-\epsilon_{\rho}>0$. Then, we find

[TABLE]

which implies A2.

• A2$\Rightarrow$B2:

In case $p_{\infty}=\infty$ condition B2 is fulfilled. Let now $p_{\infty}<\infty$ and let there be for all sub-sonic $(\rho_{l},q_{l})$ a $\rho_{-}>0$ such that

[TABLE]

that is, let A2 be true. Note that the right-hand-side is increasing in $\rho_{-}$ because its derivative is positive. Therefore taking suprema of (18) yields

[TABLE]

With Proposition B we have a list of conditions B1, B2, B3 which is equivalent to conditions A1, A2, A3 from Proposition A but now only involves the pressure function. This is fortunate as we want to prove Proposition 2, which only contains pressure functions. In principle, we could just plug the pressure functions into our conditions B and check, what generalized $\gamma$-laws are allowed. But there is more insight to be gained by simplifying the conditions further. We will do so in Propositions 5, 6 and 7. Yet this comes at a price, the next list of conditions in Proposition C below will only be sufficient conditions for Propositions A and B to hold.

Proposition 5**.**

Let $p\in C^{3}(\mathbb{R}^{+})$. Condition B1 holds if

[TABLE]

holds for all $\rho>0$.

Proof.

We note that $\left(2ff^{\prime\prime}-(f^{\prime})^{2}\right)(\rho_{l},\rho_{l})=0$ and that

[TABLE]

As $f\geq 0$, we see that $2ff^{\prime\prime}-(f^{\prime})^{2}\geq 0$ for all $\rho\geq\rho_{l}$ if $f^{\prime\prime\prime}(\rho,\rho_{l})\geq 0$ for all $\rho\geq\rho_{l}$. This was also proved in [13]. This is the case if and only if $\rho_{l}f^{\prime\prime\prime}(\rho,\rho_{l})\geq 0$. For this we find

[TABLE]

which is an affine function in $\rho_{l}$. Hence for given $\rho>0$ this function takes its minimum in one of the edges of the simplex $\{\rho_{l}\ |\ 0\leq\rho_{l}\leq\rho\}$. The values on these are given by

[TABLE]

The next condition we will replace is condition B3B$3$(b).

Proposition 6**.**

Condition B3B$3$(b) is fulfilled if either of these conditions is satisfied:

1. C3(b)

2. (i)

$\lim_{\rho\to 0}c(\rho)=0$* and $2p^{\prime}(\rho)-\rho p^{\prime\prime}(\rho)\geq 0$ for all $\rho>0$.* 2. (ii)

$0<\lim_{\rho\to 0}c(\rho)<\infty$. 3. (iii)

There is a $\eta\in(0,1)$ such that $\lim_{\rho\to 0}\rho^{\eta}c(\rho)$ exists and $0<\lim_{\rho\to 0}\rho^{\eta}c(\rho)<\infty$.

Proof.

1. (i)

Here we have

[TABLE] 2. (ii)

In this case the integral is unbounded near zero, but $c(\rho)+c(\rho_{l})$ is finite. 3. (iii)

In this case there is $a>0$ such that $\rho^{\eta}c(\rho)=a\left(1+r(\rho)\right)$ with $r(\rho)\to 0$. Then, for every $m\in\mathbb{N}$ there is $\rho_{m}>0$ such that for all $\rho<\rho_{m}$ there holds $c(\rho)\geq a\left(1-\frac{1}{m}\right)\rho^{-\eta}$. Choose now $m$ so great that

[TABLE]

and then $\rho_{m,0}$ so small that for all $\rho<\rho_{m,0}$ there holds

[TABLE]

with $\theta>0$. Define also

[TABLE]

Then we find for $\rho<\rho_{m,0}$

[TABLE]

Proposition 7**.**

All conditions only depend on differences in pressures, that is, if $c\in\mathbb{R}$ is a constant then if a pressure function $p$ satisfies our conditions, so does $p+c$. This means that if $p_{\infty}<\infty$ in condition B2, we can instead choose $p_{\infty}=0$, yielding the clearer condition

[TABLE]

Having all we need, we now sum up our findings in the following Proposition C:

Proposition C** (Proposition 8).**

The Riemann problem for arbitrary sub-sonic left- and right-hand states is well-posed if the pressure function $p$ satisfies the following conditions:

1. (C1)

both inequalities of Proposition 5:

[TABLE] 2. (C2)

one of the conditions according to Proposition 7:

1. (a)

$p\to\infty$* for $\rho\to\infty$ or* 2. (b)

$p^{\prime}(\rho)\leq-\frac{p(\rho)}{\rho}$* for all $\rho>0$,* 3. (C3)

according to Proposition 6:

1. (a)

There is $p_{0}>0$ such that $p(\rho)=-\frac{p_{0}}{\rho}+o\left(\frac{1}{\rho}\right)$ for $\rho\to 0$ or 2. (b)

one of the following conditions is fulfilled:

1. (i)

$\lim_{\rho\to 0}c(\rho)=0$* and $2p^{\prime}(\rho)-\rho p^{\prime\prime}(\rho)\geq 0$ for all $\rho>0$,* 2. (ii)

$0<\lim_{\rho\to 0}c(\rho)<\infty$, 3. (iii)

there is a $\eta\in(0,1)$ such that $\lim_{\rho\to 0}\rho^{\eta}c(\rho)$ exists and $0<\lim_{\rho\to 0}\rho^{\eta}c(\rho)<\infty$.

These conditions are positively linear in $p$ and hence all pressure functions satisfying them define a convex cone. In addition, because the integral is monotone and linear, we can also integrate over pressure functions to find new pressure functions.

With Proposition C we have a useful list of conditions which can be checked easily for any given candidate for a pressure function. The proposition now makes it easy to prove our central result Proposition 2 and we will do so now.

Proof of Proposition 2.

• •

For the first part of the claim we note that for $\mathinner{\!\left\lvert\delta\right\rvert}\leq 2$ and $\alpha>0$ the generalized $\gamma$-law fulfills the conditions in C. Note that at $\delta=0$ the conditions in C3 switch and at $\delta=-1$ the conditions in C2 switch.

• •

For the second part of the claim we observe: For $\delta>2$, choose $-c(\rho_{l})<\frac{q_{l}}{\rho_{l}}<-\frac{2}{\delta}c(\rho_{l})$, which is obviously sub-sonic. Then $L_{l}$ is strictly negative for $\rho>0$, as is easily computed. A similar range for $\frac{q_{r}}{\rho_{r}}$ shows the same for $-L_{r}$. For $\delta<-2$, choose $\sqrt{\frac{-1}{\delta+1}}c(\rho_{l})<\frac{q_{l}}{\rho_{l}}<c(\rho_{l})$ and $L_{l}$ is strictly positive for $\rho>0$ and similarly for $-L_{r}$.∎

Example 9**.**

As just proved, the functions

[TABLE]

are all valid pressure functions. Even a and more exotic

[TABLE]

might be possible.

2.2 Extension to junctions

Our findings guarantee well-posedness of Riemann problems with sub-sonic initial conditions under certain conditions on the pressure function. We now examine the coupling conditions (4) and (5) and do so by considering the generalized Riemann problem at the junction in accordance with [1, 2, 7]. Consider a junction with incoming pipes indexed by $e\in{I}$ and outgoing pipes indexed by $f\in{O}$. At the junction-facing end of each pipe there are initial states $U_{i}=(\rho_{i},q_{i})$ for each $i\in I$ and $i\in O$, respectively. To make things tractable, we restrict the solution of the junction Riemann problem to be of this form: In each pipe $i\in I\cup O$, there appears exactly one new state $V_{i}$ next to the junction such that the $V_{i}$ satisfy the junction conditions and are connected to their respective $U_{i}$ by a shock or a rarefaction wave. On ingoing pipes these must be $1$-waves, on outgoing pipes these must be $2$-waves. A sketch of this is shown in figure 2.

This is not a great restriction as this is the only solution structure found in pipelines with low Mach number, yet for different solutions, see [13]. For these new states to appear at all, the wave speed between $V_{i}$ and $U_{i}$ must be negative on ingoing pipes and positive on outgoing ones. To keep things simple, we examine a single ingoing pipe $i$ with initial condition $U_{i}$ and new state $V_{i}=(\rho,L_{l}(\rho;U_{i}))$.

Proposition 10**.**

Let Proposition B and hence A be fulfilled for the pressure. Let $\rho_{i,\text{min}}$ be such that $L_{l}^{\prime}(\rho_{i,\text{min}};U_{i})=0$, if it exists, otherwise let $\rho_{i,\text{min}}=0$. The wave speed of the wave between $U_{i}$ and $V_{i}$ is negative if and only if $\rho>\rho_{i,\text{min}}$.

Proof.

One easily sees that $\lambda_{1}(\rho,L_{l}(\rho;U_{i}))=L_{l}^{\prime}(\rho;U_{i})$. Also using the definitions (7) and (12) and the concavity of $L_{l}$, one can see that $\rho_{i,\text{min}}$ is always unique. Because $L_{l}$ is concave, $L_{l}^{\prime}$ is decreasing, so $\lambda_{1}(\rho,L_{l}(\rho;U_{i}))<0$ if and only if $\rho>\rho_{i,\text{min}}$.

• $\Leftarrow$:

Let $\rho>\rho_{i,\text{min}}$. Then $\lambda_{1}(V_{i})<0$. If $V_{i}$ is in the rarefaction part of $L_{l}(\rho;U_{i})$ the wave speed is negative. If $V_{i}$ is in the shock part of $L_{l}(\rho;U_{i})$, we compute:

[TABLE]

since $L_{l}$ is concave. So

[TABLE]

because $U_{i}$ is sub-sonic.

• $\Rightarrow$:

Let $\rho\leq\rho_{i,\text{min}}$. Because $U_{i}$ is sub-sonic, we have $\rho\leq\rho_{i,\text{min}}<\rho_{i}$ and are dealing with a rarefaction wave. But then the wave speed is just given by $\lambda_{1}(\rho,L_{l}(\rho;U_{i})$ which is non-negative. ∎

For outgoing pipes one defines $\rho_{i,\text{min}}\geq 0$ such that $L_{r}^{\prime}(\rho_{i,\text{min}};U_{i})=0$, if possible, or $\rho_{i,\text{min}}=0$ otherwise and obtains in the same way

Proposition 11**.**

Let Proposition B be fulfilled for the pressure. On an outgoing pipe the wave speed between $U_{i}$ and $V_{i}=(\rho,L_{r}(\rho;U_{i}))$ is positive if and only if $\rho>\rho_{i,\text{min}}$.

Coming back to junctions we define the junction minimal density as

[TABLE]

A solution $(V_{i})_{i\in I\cup O}$ to the junction Riemann problem is admissible if and only if the density $\rho$ at the junction fulfills $\rho>\rho_{\text{min}}$. If it is not admissible, then in at least one pipe there is a super-sonic gas flow.

Note that a usual Riemann problem with sub-sonic initial conditions can be treated as a junction with one ingoing and one outgoing pipe. In this case only one new state $V$ is created and the admissibility criterion guarantees that $\lambda_{1}(V)<0<\lambda_{2}(V)$. For the sake of completeness, we classify the solutions to the usual Riemann problem by wave types in Table 1.

2.3 Additional constraints for consistency

We have found a complete list of conditions for a usual Riemann problem and a junction Riemann problem to be well-posed. For the usual Riemann problem this is enough, as just noted. Yet our conditions do not guarantee that the solution states at the junction are sub-sonic, as only one of each eigenvalue is restricted. Since sub-sonic states are necessary at the junction, we introduce an additional condition by hand. We again focus on one incoming pipe $i$ with initial condition $U_{i}$ and new state $V_{i}=(\rho,L_{l}(\rho;U_{i}))$. We define

[TABLE]

and further

[TABLE]

Although there is no knowledge about $\lambda_{2}(V_{i})$ as it is not conveniently given by the derivative of a Lax curve, we still have $\lambda_{2}(V_{i})=\lambda_{1}(V_{i})+2c(\rho)>\lambda_{1}(V_{i})$, so there holds at least $\rho_{i,\text{max}}>\rho_{i,\text{min}}$ in every pipe. But there is no guarantee that $\rho_{\text{max}}>\rho_{\text{min}}$.

3 Coupling to power grids

In this section, we focus on the coupling of the gas network to power grids. Following the ideas in [16, 27], the coupling is done extending the already existing coupling conservation of fluxes (5) by a gas power plant. To model the power grid, we use the well-known powerflow equations [11], which describe real and reactive power nodewise. We have to investigate again under which assumptions this kind of coupling is well-posed.

3.1 Power flow model

Power grids are usually modeled as a graph, whose nodes (buses) are indexed by $k$ and whose edges (transmission lines) are indexed by the indices of the nodes connected by the edge. For each node $k\in V_{{\text{PG}}}$ of the power grid (PG) there are four (time-dependent) state variables, namely real power $P_{k}(t)$, reactive power $Q_{k}(t)$, voltage magnitude $\mathinner{\!\left\lvert V_{k}(t)\right\rvert}$ and voltage angle $\phi_{k}(t)$. We consider three different types of nodes: load/PQ buses ($P$ und $Q$ are given), generater/PV buses ($P$ and $\mathinner{\!\left\lvert V\right\rvert}$ are given) and a single slack bus ($\mathinner{\!\left\lvert V\right\rvert}$ and $\phi$ are given). Thus we always have two unknowns per node $k\in V_{\text{PG}}$. Accordingly, we apply the powerflow equations for real and reactive power for each node:

[TABLE]

where $G_{kj}$ is the real part of the entry $y_{kj}$ in the bus admittance matrix and $B_{kj}$ is the imaginary part. For $k\neq j$ we consider $G_{kj}$ and $B_{kj}$ as properties of the arc/transmission line connecting nodes $k$ and $j$. $G_{kk}$ and $B_{kk}$ are considered as node properties.

3.2 Gas-to-power coupling

We focus on a junction with one incoming pipeline, one outgoing pipeline and an outlet that draws a set amount of flow $\varepsilon$, which is converted to electric power. The coupling condition is taken from [16] and reads

[TABLE]

where the outflow $\varepsilon$ is non-negative. The fuel is converted to gas via the following heat rate formula, taken from [27],

[TABLE]

where $P$ is the real power and $a_{0}$, $a_{1}$, $a_{2}$ are constants.

To find what $\varepsilon$ are allowed, we examine another generalized Riemann problem at the junction. We start with two constant states on the ingoing and outgoing pipes and solve the Riemann problem where the flow $q$ jumps by an amount $\varepsilon$ at the junction. We again demand the waves to leave the junction, giving rise to a picure like that in Figure 4,

where next to the junction in red two new states $V_{l}$, $V_{r}$ appear that fulfill the coupling conditions (25). The analysis of this setting is similar to a three-way junction but with the flow on one outgoing pipe fixed to $\varepsilon$. We must now solve the equation

[TABLE]

which is similar to equation (10) of a usual (two-way) Riemann problem but with a non-zero right-hand side. We even can forego the additional consistency constraints of section 2.3: One half of the eigenvalues is on the right side of zero due to the $\rho_{\text{min}}$-criterion, which we need here again to make the waves have the right direction.

[TABLE]

Instead of invoking $\rho_{\text{max}}$, we note that $\varepsilon\geq 0$ and compute for the remaining two eigenvalues

[TABLE]

and

[TABLE]

Remark 12**.**

In our setting this is sufficient. If we also considered power-to-gas plants, we would need to also treat negative $e$. In this case, we would have to adhere to the $\rho_{\text{max}}$-criterion (23) again.

Solution structure for different outflows $\varepsilon$

For the usual Riemann problem we had a unique non-zero solution due to our findings in Section 2.1. For $\varepsilon>0$ we now have two solutions to $L_{l}(\rho)-L_{r}(\rho)=\varepsilon$, one of which is not admissible as it lies to the left of the maximum and hence has non-negative derivative (which means it would be super-sonic). As $L_{l}-L_{r}$ is decreasing in the admissible regime, greater $\varepsilon$ result in smaller $\rho$. One of the two solution structures is given in table 2.

4 Numerical results

Within the following numerical examples, we consider two different discretization schemes. The first is a third-order CWENO scheme (CWENO3) with suitable boundary treatment [20, 23], which relies on a local Lax-Friedrichs flux function for the inner discretization points of each pipe and handles coupling points by explicitly solving equation (27).

The second scheme is an implicit box scheme (IBOX) [21], suitable for sub-sonic flows. For a general system of balance laws

[TABLE]

the considered scheme reads

[TABLE]

Here, $\Delta t$ and $\Delta x$ are the temporal and spatial mesh size, respectively, and the numerical approximation is thought in the following sense:

[TABLE]

The implicit box scheme has to obey an inverse CFL condition [21], which is beneficial for problems with large characteristic speeds whereas the solution is merely quasi-stationary. This is usually the case for daily operation tasks in gas networks and therefore motivates the choice of this scheme for the real-world scenario below.

The first test example in section 4.1 is supposed to demonstrate the different cases revealed in the analysis above (section 3.2). Further, since the applied implicit box scheme does not explicitly make use of any Riemann solver, this scenario is also considered as a numerical validation of its applicability, where the CWENO3 scheme with Riemann solver at the junction serves as reference.

Within the second example (section 4.2) we briefly demonstrate the differences resulting from various pressure functions, which are all covered by our theoretical results above.

The third example (section 4.3) considers a more complex scenario: An increasing power demand within the power grid leads to an increasing fuel demand of a gas-to-power generator and further to a significant pressure drop in the gas network.

4.1 Validation

We consider the isentropic Euler equations with pressure law $p(\rho)=\kappa\rho^{\gamma}$ and parameters $\kappa=1.0$, $\gamma=1.4$, and a Riemann problem with left state $U_{l}=\begin{pmatrix}4.0\\ 1.0\end{pmatrix}$ and right state $U_{r}=\begin{pmatrix}3.0\\ -1.0\end{pmatrix}$. Further we assume a gas demand $\varepsilon$ at the coupling point of the two states (here $x=0$). Then, from table 2, we get the following solution structure:

• •

s-s solution for $\varepsilon\leq 0.57877$,

• •

r-s solution for $0.57877\leq\varepsilon\leq 3.0594$,

• •

r-r solution for $3.0594\leq\varepsilon$.

Further, one can easily compute $\rho_{1,\min}\approx 1.8819$, $\rho_{2,\min}\approx 1.5041$, and therewith $\rho_{\min}\approx 1.8819$ and maximum gas demand $\varepsilon\approx 4.3892$. We will consider the numerical simulation of the described setting until time $t=0.1$ for $\varepsilon\in\{0.25,1.75,3.25\}$ and the following discretization parameters:

• •

CWENO3: $\Delta t=5\cdot 10^{-5}$, $\Delta x=5\cdot 10^{-4}$,

• •

IBOX: $\Delta t=5\cdot 10^{-4}$, $\Delta x=5\cdot 10^{-5}$.

The different choices result from the (usual) CFL condition the explicit CWENO3 scheme has to obey, in contrast to the inverse CFL condition of the IBOX scheme. Figures 5 to 7 show the computed densities at the final time. Both schemes show the correct solution structure (shock/rarefaction waves), where CWENO3 expectably achieves the sharper resolution.

4.2 Different pressure laws

In our second test case, we apply various pressure laws, which are all covered by our theoretical study, and are interested in the different dynamics one may observe even on a single pipeline. Therefore, we consider a single pipe with length $l=0.1$ and the following pressure laws:

• •

$p(\rho)=\frac{1}{\gamma}\rho^{\gamma}$ with $\gamma=1.4$ (“gamma law”),

• •

$p(\rho)=\frac{-1}{\rho}$ (“inverse”, corresponding to $\gamma=-1$),

• •

$p(\rho)=\ln(\rho)$ (“logarithmic”),

• •

$p(\rho)=\frac{1}{10}\sum\limits_{i=1}^{10}\frac{\rho^{1+i/5}}{1+i/5}$ (“sum of gamma laws”).

Note that all considered pressure functions are scaled in such a way that $p^{\prime}(\rho=1)=1$. Initially, we have $\rho=1$ and $q=0$ in the whole pipe. Further, we fix $\rho=1$ on the left-hand boundary, whereas $q$ at the right-hand boundary linearly increases from [math] to $0.2$ until time $t=0.1$ and stays constant afterwards until the final time $t=0.5$. We approximate the solution to this problem by CWENO3 with discretization parameters $\Delta t=5\cdot 10^{-4}$ and $\Delta x=10^{-3}$. The variety of the resulting dynamics is demonstrated in figures 8 and 9, which show the density in the pipeline at times $t=0.25$ and $t=0.5$, respectively.

4.3 Coupled gas and power grid

We consider the network(s) depicted in figure 10, containing a power grid from the example “case9” of the MATPOWER Matlab programming suite [26] and a small part of the GasLib-40 network [18], extended by a compressor station in front of node S17 and a gas-to-power generator between S4 and N1, providing the necessary power at the latter node. The parameters of the gas network are given in table 3, the parameters of the power grid in table 4.

Within each pipe $e\in E_{\text{Gas}}$ of the gas network (with length $l_{e}$, diameter $d_{e}$, cross section $A_{e}$, roughness $k_{e}$) we consider the isothermal Euler equations with acoustic speed $c=340\frac{\text{m}}{\text{s}}$. The source term in the momentum equation is given by

[TABLE]

with friction factor $\lambda(q)$, which is determined by the Prandtl-Colebrook formula:

[TABLE]

with Reynolds number

[TABLE]

and dynamic viscosity $\eta=10^{-5}\frac{\text{kg}}{\text{ms}}$.

Initially, the gas network is in a stationary state: The pressure at S5 is fixed at 60bar, the ouflow at S25 is $q=100\frac{\text{m}^{3}}{\text{s}}\cdot\frac{\rho_{0}}{A_{e}}$ with $\rho_{0}=0.785\frac{\text{kg}}{\text{m}^{3}}$, and there is an additional gas consumption at S4 resulting from the gas to power transformation ($a_{0}=2$, $a_{1}=5$, $a_{2}=10$) due to the power demand at the slack bus N1. The initial (stationary) state of the power grid is determined by boundary conditions given in table 5.

In the course of the simulation, the power (and reactive power) demand at N5 is(/are) linearly increased between $t=1$ hour and $t=1.5$ hours from $0.9$ p.u. to $1.8$ p.u. (reactive power from $0.3$ p.u. to $0.6$ p.u.) - see also figure 11 (left). Accordingly, the power demand at the slack bus N1 increases (see figure 11 (right)) and therewith the gas consumption at S4, which also results in an increase of the inflow at S5 (see figure 12). Due to the increased flow values, the pressure in the gas network decreases (see figure 13 for the pressure at the nodes S20 and S25).

5 Conclusion and future work

We have presented a coupled model for gas and power allowing for a mathematically well-defined transition from gas to power. The framework involves also the consideration of non-standard pressure functions. Various simulation results show the properties of the presented approach.

Future work includes the investigation of optimal control problems for the coupled model as for example the control of compressor stations [14] or the inclusion of uncertain customer demands [10]. This will immediately lead to the class of nonlinear (stochastic) optimization problems and tailored solution techniques.

Acknowledgments

The authors gratefully thank the BMBF project ENets (05M18VMA) for the financial support.

Appendix A Proof of Lemma 3

The following two technical results are the key ingredients to prove Lemma 3.

Lemma 13**.**

Let $g\in C^{1}(\mathbb{R}^{+},\mathbb{R}^{+})$be a non-negative function, $g\geq 0$, and let $G$ be given by $G(\rho)=\int_{\rho}^{\rho_{l}}g(s)\operatorname{d\!}s$. Then There holds

If $\rho^{2}g(\rho)\overset{\rho\to 0}{\to}0$, then $\rho G(\rho)\overset{\rho\to 0}{\to}0$. 2. 2.

If $\rho G(\rho)\overset{\rho\to 0}{\to}0$, then $\liminf_{\rho\to 0}\rho^{2}g(\rho)=0$.

Proof.

By assumption $\rho^{2}g(\rho)\overset{\rho\to 0}{\to}0$. For $m\in\mathbb{N}$ choose $\rho_{m}>0$ such that $\rho^{2}g(\rho)\leq\frac{1}{m\rho^{2}}$ for $\rho<\rho_{m}$. Now choose $\rho_{m,0}<\rho_{m}$ so small that

[TABLE]

Then, for $\rho<\rho_{m,0}$, there holds

[TABLE]

for small enough $\rho$. Therefore $\lim_{\rho\to 0}\rho\int_{\rho}^{\rho_{l}}g(s)\operatorname{d\!}s\leq\frac{2}{m}$ for all $m\in\mathbb{N}$. As $g\geq 0$, we also have $0\leq\lim_{\rho\to 0}\rho\int_{\rho}^{\rho_{l}}g(s)\operatorname{d\!}s$. Summarizing, we get $\lim_{\rho\to 0}\rho G(\rho)=\lim_{\rho\to 0}\rho\int_{\rho}^{\rho_{l}}g(s)\operatorname{d\!}s=0$. 2. 2.

We prove by contradiction: Assume there are $\rho_{0}>0$ and $a>0$ such that $\rho^{2}g(\rho)\geq a$ for $\rho<\rho_{0}$. Then $g(\rho)\geq\frac{a}{2}\frac{1}{\rho^{2}}$ for such $\rho$. Therefore

[TABLE]

which contradicts $\rho G(\rho)\overset{\rho\to 0}{\longrightarrow}0$.∎

Lemma 14**.**

Let $g\in C^{1}(\mathbb{R}^{+},\mathbb{R}^{+})$, $g\geq 0$ and $\liminf_{\rho\to 0}\rho^{2}g(\rho)=0$. Let also $\left(\rho^{2}g(\rho)\right)^{\prime}\geq 0$. Then, $\limsup_{\rho\to 0}\rho^{2}g(\rho)=\liminf_{\rho\to 0}\rho^{2}g(\rho)=0$.

Proof.

We prove by contradiction. Let $\limsup_{\rho\to 0}\rho^{2}g(\rho)>\liminf\rho^{2}g(\rho)$. Then it is easily seen that $\liminf_{\rho\to 0}\left(\rho^{2}g(\rho)\right)^{\prime}=-\infty<0$ resulting in a contradiction. ∎