An asymptotic preserving mixed finite element method for wave propagation in pipelines

TL;DR

This paper develops an asymptotic preserving mixed finite element method for wave propagation in pipelines, effectively capturing the transition to steady states and parabolic limits with proven stability and accuracy.

Contribution

It introduces a novel mixed finite element approach that preserves asymptotic behavior in the parabolic limit and demonstrates uniform exponential stability regardless of discretization.

Findings

Method is asymptotic preserving and uniformly exponentially stable.

Numerical tests confirm sharpness of theoretical estimates.

Approach effectively handles parameter-dependent damped hyperbolic equations.

Abstract

We consider a parameter dependent family of damped hyperbolic equations with interesting limit behavior: the system approaches steady states exponentially fast and for parameter to zero the solutions converge to that of a parabolic limit problem. We establish sharp estimates and elaborate their dependence on the model parameters. For the numerical approximation we then consider a mixed finite element method in space together with a Runge-Kutta method in time. Due to the variational and dissipative nature of this approximation, the limit behavior of the infinite dimensional level is inherited almost automatically by the discrete problems. The resulting numerical method thus is asymptotic preserving in the parabolic limit and uniformly exponentially stable. These results are further shown to be independent of the discretization parameters. Numerical tests are presented for a simple model problem which illustrate that the derived estimates are sharp in general.

Tables

Table 1. Table 1: Distance ‖ p h , τ ϵ ​ ( t ) − p ¯ h ‖ 2 + ϵ 2 ​ ‖ m h , τ ϵ − m ¯ h ‖ 2 ≤ C ​ e − γ ​ t superscript norm superscript subscript 𝑝 ℎ 𝜏 italic-ϵ 𝑡 subscript ¯ 𝑝 ℎ 2 superscript italic-ϵ 2 superscript norm superscript subscript 𝑚 ℎ 𝜏 italic-ϵ subscript ¯ 𝑚 ℎ 2 𝐶 superscript 𝑒 𝛾 𝑡 \|p_{h,\tau}^{\epsilon}(t)-\bar{p}_{h}\|^{2}+\epsilon^{2}\|m_{h,\tau}^{\epsilon}-\bar{m}_{h}\|^{2}\leq Ce^{-\gamma t} of the numerical solution to the discrete steady state for different values of ϵ italic-ϵ \epsilon and times t = 0 , 0.1 , 0.5 , 1.0 𝑡 0 0.1 0.5 1.0 t=0,0.1,0.5,1.0 and estimated exponential convergence rate γ 𝛾 \gamma . Discretization parameters were set to h = 0.01 ℎ 0.01 h=0.01 and τ = 10 − 5 𝜏 superscript 10 5 \tau=10^{-5} .

$t\backslash ϵ$	1/4	1/8	1/16	1/32	1/64	1/128
0.0	5.00e-01	5.00e-01	5.00e-01	5.00e-01	5.00e-01	5.00e-01
0.1	2.72e-01	9.09e-02	7.26e-02	7.02e-02	6.96e-02	6.95e-02
0.5	3.56e-04	5.35e-06	1.94e-05	2.42e-05	2.54e-05	2.57e-05
1.0	8.51e-08	2.71e-11	6.64e-10	1.13e-09	1.28e-09	1.32e-09
$\gamma$	15.59	23.64	20.44	19.90	19.78	19.75

Table 2. Table 2: Error ‖ p h , τ ϵ ​ ( t k ) − p h , τ 0 ​ ( t k ) ‖ 2 + ∑ j = 1 k a ¯ ​ ‖ m h , τ ϵ ​ ( t j ) − m h , τ 0 ​ ( t j ) ‖ 2 = O ​ ( ϵ α ) superscript norm superscript subscript 𝑝 ℎ 𝜏 italic-ϵ superscript 𝑡 𝑘 superscript subscript 𝑝 ℎ 𝜏 0 superscript 𝑡 𝑘 2 superscript subscript 𝑗 1 𝑘 ¯ 𝑎 superscript norm superscript subscript 𝑚 ℎ 𝜏 italic-ϵ superscript 𝑡 𝑗 superscript subscript 𝑚 ℎ 𝜏 0 superscript 𝑡 𝑗 2 𝑂 superscript italic-ϵ 𝛼 \|p_{h,\tau}^{\epsilon}(t^{k})-p_{h,\tau}^{0}(t^{k})\|^{2}+\sum_{j=1}^{k}\underline{a}\|m_{h,\tau}^{\epsilon}(t^{j})-m_{h,\tau}^{0}(t^{j})\|^{2}=O(\epsilon^{\alpha}) between the discrete approximations for the hyperbolic problem and the parabolic limit problem for different values of ϵ italic-ϵ \epsilon and time steps t k superscript 𝑡 𝑘 t^{k} and observed convergence rate α 𝛼 \alpha . Discretization with h = 0.01 ℎ 0.01 h=0.01 and τ = 10 − 5 𝜏 superscript 10 5 \tau=10^{-5} .

$t^k\backslash ϵ$	1/4	1/8	1/16	1/32	1/64	1/128	$\alpha$
0.1	9.81e-02	3.47e-02	9.41e-03	2.38e-03	5.89e-04	1.39e-04	1.87
0.5	1.18e-01	3.58e-02	9.44e-03	2.39e-03	5.89e-04	1.39e-04	1.93
1.0	1.18e-01	3.58e-02	9.44e-03	2.39e-03	5.89e-04	1.39e-04	1.93

Table 3. Table 3: Error ‖ p h , τ ϵ ​ ( t k ) − p h , τ 0 ​ ( t k ) ‖ 2 + ∑ j = 1 k a ¯ ​ ‖ m h , τ ϵ ​ ( t j ) − m h , τ 0 ​ ( t j ) ‖ 2 = O ​ ( ϵ α ) superscript norm superscript subscript 𝑝 ℎ 𝜏 italic-ϵ superscript 𝑡 𝑘 superscript subscript 𝑝 ℎ 𝜏 0 superscript 𝑡 𝑘 2 superscript subscript 𝑗 1 𝑘 ¯ 𝑎 superscript norm superscript subscript 𝑚 ℎ 𝜏 italic-ϵ superscript 𝑡 𝑗 superscript subscript 𝑚 ℎ 𝜏 0 superscript 𝑡 𝑗 2 𝑂 superscript italic-ϵ 𝛼 \|p_{h,\tau}^{\epsilon}(t^{k})-p_{h,\tau}^{0}(t^{k})\|^{2}+\sum_{j=1}^{k}\underline{a}\|m_{h,\tau}^{\epsilon}(t^{j})-m_{h,\tau}^{0}(t^{j})\|^{2}=O(\epsilon^{\alpha}) between the discrete approximations for the hyperbolic problem and the parabolic limit problem for time t k = 1 superscript 𝑡 𝑘 1 t^{k}=1 and different values of ϵ italic-ϵ \epsilon and the discretization parameters h ℎ h and τ 𝜏 \tau .

$ϵ$	1/4	1/8	1/16	1/32	1/64	1/128	$\alpha$
$h=0.010,\tau ={10}^{-5}$	1.18e-01	3.58e-02	9.44e-03	2.39e-03	5.89e-04	1.39e-04	1.93
$h=0.002,\tau ={10}^{-5}$	1.18e-01	3.58e-02	9.44e-03	2.39e-03	5.89e-04	1.39e-04	1.99
$h=0.010,\tau ={10}^{-6}$	1.18e-01	3.58e-02	9.46e-03	2.40e-03	6.00e-04	1.49e-04	1.90
$h=0.002,\tau ={10}^{-6}$	1.18e-01	3.58e-02	9.46e-03	2.40e-03	6.00e-04	1.49e-04	1.90

Equations

$${\partial }_t{p}^{ϵ}+{\partial }_x{m}^{ϵ}$$

$${ϵ}^2{\partial }_t{m}^{ϵ}+{\partial }_x{p}^{ϵ}+a{m}^{ϵ}$$

$$\Vert p^{ϵ}(t)-p^0(t){\Vert }^2+{\int }_0^t\Vert m^{ϵ}(s)-m^0(s){\Vert }^{2}ds\le C{ϵ}^2$$

$$\Vert p^{ϵ}(t)-\bar{p}{\Vert }^2+{ϵ}^2\Vert m^{ϵ}(t)-\bar{m}{\Vert }^2\le C{e}^{-\gamma t}$$

$$\Vert p_{h,\tau }^{ϵ}(t)-p_{h,\tau }^0(t){\Vert }^2+{\int }_0^t\Vert m_{h,\tau }^{ϵ}(s)-m_{h,\tau }^0(s){\Vert }^{2}ds\le C{ϵ}^2$$

$$\Vert p_{h,\tau }^{ϵ}(t)-{\bar{p}}_h{\Vert }^2+{ϵ}^2\Vert m_{h,\tau }^2(t)-{\bar{m}}_h{\Vert }^2\le C{e}^{-\gamma t}$$

$${\partial }_t{p}^{ϵ}(x,t)+{\partial }_x{m}^{ϵ}(x,t)$$

$${ϵ}^2{\partial }_t{m}^{ϵ}(x,t)+{\partial }_x{p}^{ϵ}(x,t)+a(x)m^{ϵ}(x,t)$$

$$p^{ϵ}(0,t)=g_0,p^{ϵ}(x,t)$$

$$(p,m)\in C^1({\mathbb{R}}^{+};L^2(0,1)\times L^2(0,1))\times C({\mathbb{R}}^{+};H^1(0,1)\times H^1(0,1))$$

$$p^0\in C^1({\mathbb{R}}^{+};L^2(0,1))\times C({\mathbb{R}}^{+};H^1(0,1)),m^0\in C({\mathbb{R}}^{+};L^2(0,1))$$

$${\partial }_x\bar{m}(x)$$

$${\partial }_x\bar{p}(x)+a(x)\bar{m}(x)$$

$$\bar{p}(0)=g_0,\bar{p}(1)$$

$$\Vert p^{ϵ}(t)-\bar{p}{\Vert }^2+{ϵ}^2\Vert m^{ϵ}(t)-\bar{m}{\Vert }^2+2{\int }_0^t\underline{a}\Vert m^{ϵ}(s)-\bar{m}{\Vert }^{2}ds$$

$$\le \Vert p_0-\bar{p}{\Vert }^2+{ϵ}^2\Vert m_0-\bar{m}{\Vert }^2.$$

$$\Vert {\partial }_t{p}^{ϵ}(t){\Vert }^2+{ϵ}^2\Vert {\partial }_t{m}^{ϵ}(t){\Vert }^2+2{\int }_0^t\underline{a}\Vert {\partial }_t{m}^{ϵ}(s){\Vert }^{2}ds$$

$$\le \Vert {\partial }_x{m}_0{\Vert }^2+\frac{1}{{ϵ}^2}\Vert {\partial }_x{p}_0+a{m}_0{\Vert }^2.$$

$$\frac{1}{2}\frac{d}{dt}\Vert p^{ϵ}{\Vert }^2$$

$$\frac{d}{dt}\Vert p^{ϵ}{\Vert }^2+{ϵ}^2\frac{d}{dt}\Vert m^{ϵ}{\Vert }^2\le -2\underline{a}\Vert m^{ϵ}{\Vert }^2.$$

$$\Vert {\partial }_t{p}^{ϵ}(t){\Vert }^2+{ϵ}^2\Vert {\partial }_t{m}^{ϵ}(t){\Vert }^2+2{\int }_0^t\underline{a}\Vert {\partial }_t{m}^{ϵ}(t){\Vert }^2$$

$$\le \Vert {\partial }_t{p}^{ϵ}(0){\Vert }^2+{ϵ}^2\Vert {\partial }_t{m}^{ϵ}(0){\Vert }^2.$$

$$\Vert p^{ϵ}(t)-p^0(t){\Vert }^2+{\int }_0^t\underline{a}\Vert m^{ϵ}(s)-m^0(s){\Vert }^{2}ds$$

$$\le \frac{{ϵ}^4}{2{\underline{a}}^2}(\Vert {\partial }_x{m}_0{\Vert }^2+\frac{1}{{ϵ}^2}\Vert {\partial }_x{p}_0+a{m}_0{\Vert }^2).$$

$${\partial }_t{r}^{ϵ}+{\partial }_x{w}^{ϵ}$$

$${\partial }_x{r}^{ϵ}+a{w}^{ϵ}$$

$$\frac{1}{2}\frac{d}{dt}\Vert r^{ϵ}(t){\Vert }^2+\underline{a}\Vert w^{ϵ}(t){\Vert }^2$$

$$\Vert r^{ϵ}(t){\Vert }^2+{\int }_0^t\underline{a}\Vert w^{ϵ}(s){\Vert }^{2}ds\le \Vert r^{ϵ}(0){\Vert }^2+\frac{{ϵ}^4}{\underline{a}}{\int }_0^t\Vert {\partial }_t{m}^{ϵ}(s){\Vert }^{2}ds.$$

$$\Vert p^{ϵ}(t)-\bar{p}{\Vert }^2+{ϵ}^2\Vert m^{ϵ}(t)-\bar{m}{\Vert }^2\le C{e}^{-\gamma (t-s)}(\Vert p^{ϵ}(s)-\bar{p}{\Vert }^2+{ϵ}^2\Vert m^{ϵ}(s)-\bar{m}{\Vert }^2)$$

$${\partial }_t{\pi }^{ϵ}+{\partial }_x{\mu }^{ϵ}$$

$${\partial }_t{\mu }^{ϵ}+{\partial }_x{\pi }^{ϵ}+\frac{a}{ϵ}{\mu }^{ϵ}$$

$$\Vert {\pi }^{ϵ}(\tau )-\bar{\pi }{\Vert }^2+\Vert {\mu }^{ϵ}(\tau )-\bar{\mu }{\Vert }^2\le C{e}^{-cϵ(\tau -\sigma )}(\Vert {\pi }^{ϵ}(\sigma )-\bar{\pi }{\Vert }^2+\Vert {\mu }^{ϵ}(\sigma )-\bar{\mu }{\Vert }^2).$$

$$Q_h=\{q\in L^2(0,1):q{|}_e\in P_0(e)\}\text{mbox}and{V}_h=\{v\in C[0,1]:v{|}_e\in P_1(e)\}$$

$$({\bar{\partial }}_{\tau }{p}_{h,\tau }^{ϵ}(t^k),q_h)+({\partial }_x{m}_{h,\tau }^{ϵ}(t^k),q_h)$$

$${ϵ}^2({\bar{\partial }}_{\tau }{m}_{h,\tau }^{ϵ}(t^k),v_h)-(p_{h,\tau }^{ϵ}(t^k),{\partial }_x{v}_h)+(a{m}_{h,\tau }^{ϵ}(t^k),v_h)$$

$$\Vert p_{h,\tau }^{ϵ}(t^k)-{\bar{p}}_h{\Vert }^2+{ϵ}^2\Vert m_{h,\tau }^{ϵ}(t^k)-{\bar{m}}_h{\Vert }^2+2\underline{a}\sum _{j=1}^k\tau \Vert m_{h,\tau }^{ϵ}(t^j)-{\bar{m}}_h{\Vert }^2$$

Full text

1 Introduction

Pipeline networks in gas or water supply systems are usually made up of rather long pipes and the time scales of interest are typically large as well. The propagation of pressure waves in such long pipes may then be described by a hyperbolic system

[TABLE]

together with appropriate initial and boundary conditions. Here $p^{\epsilon}$ corresponds to the pressure, $m^{\epsilon}$ to the momentum or mass flux, and $a$ is a generalized friction coefficient which encodes information about the pipe diameter and roughness. This system can be derived by a parabolic rescaling $t=\tilde{t}\epsilon^{2}$, $x=\tilde{x}\epsilon$ of the physical space and time variables $\tilde{x},\tilde{t}$ from the Euler equations or the shallow water equations under some simplifying assumptions BrouwerGasserHerty11 ; Osi96 and $\epsilon$ can be assumed to be small.

The parameter dependent hyperbolic problem (1)–(2) has an interesting limit behavior for long time $t\to\infty$ and in the parabolic limit $\epsilon\to 0$ which has been studied intensively in the literature BrouwerGasserHerty11 ; GattiPata06 ; HsiaoLiu92 ; LopezGomez97 ; RauchTaylor74 ; Zuazua88 . Many interesting results are available even for more general problems including the isentropic Euler equations with damping and rather general hyperbolic systems DafermosPan09 ; MarRub00 . In this note, we contribute to this active research field by establishing the following theoretical results:

(R1)

For $\epsilon\to 0$, the solutions $(p^{\epsilon},m^{\epsilon})$ of (1)–(2) converge to the solution $(p^{0},m^{0})$ of the corresponding parabolic limit problem and

[TABLE]

with a constant $C$ that is uniform in $\epsilon$ and independent of time $t\geq 0$.

(R2)

Assume that the boundary values are kept constant. Then for any $0\leq\epsilon\leq 1$ the solutions $(p^{\epsilon},m^{\epsilon})$ converge to the same steady state $(\bar{p},\bar{m})$ and

[TABLE]

with constants $C$ and $\gamma>0$ that are independent of $t\geq 0$ and $\epsilon$.

Our proofs are based on careful energy estimates that explicitly take into account the dependence on the parameter $\epsilon$. As a consequence, the results not only hold for single pipes but can be extended without much difficulty to pipeline networks.

Due to the many important applications, the systematic approximation of parameter dependent hyperbolic problems and, in particular, the preservation of asymptotic stability have been investigated intensively as well BuetDepresFranck12 ; DuMaTuBe15 ; ErvedozaZuazua09 ; GalHerHurLeRou06 ; GosseToscani02 . For the discretization of the model problem (1)–(2) we here consider a mixed finite element method in space combined with an implicit Runge-Kutta time-stepping scheme. The resulting method can be shown to exactly conserve mass and to be slightly dissipative in energy, thus capturing the relevant physical behavior EggerKugler15 . In this paper, we additionally establish the following properties:

(R3)

The scheme is asymptotic preserving, i.e., the solutions $(p_{h,\tau}^{\epsilon},m_{h,\tau}^{\epsilon})$ converge with $\epsilon\to 0$ to the solution $(p_{h,\tau}^{0},m_{h,\tau}^{0})$ of the parabolic limit problem, and

[TABLE]

with $C$ independent of $\epsilon$ and of the discretization parameters $h$ and $\tau$.

(R4)

The method is uniformly exponentially stable, i.e., for constant boundary data the solutions $(p_{h,\tau}^{\epsilon},m_{h,\tau}^{\epsilon})$ converge towards steady state $(\bar{p}_{h},\bar{m}_{h})$ and

[TABLE]

with $C$ and $\gamma>0$ independent of $\epsilon$ and the discretization parameters $h$, $\tau$.

The numerical method is also well-balanced in the sense that it automatically provides a stable approximation $(\bar{p}_{h},\bar{m}_{h})$ for the corresponding stationary problem. Since the proposed discretization strategy is of variational and dissipative nature, the above assertions can be proven with only slight modification of the energy arguments used on the continuous level. In summary, we thus obtain uniformly stable and accurate approximations for the parameter dependent problem (1)–(2) that capture all relevant physical and mathematical properties of the underlying system.

The remainder of this note is organized as follows: In Section 2, we prove the assertions (R1) and (R2) for the case of a single pipe. Section 3 is then concerned with the numerical approximation and the proof of assertions (R3) and (R4) for a single pipe. In Section 4, we briefly indicate how the results can be generalized with minor modifications to pipe networks. In Section 5, we discuss in detail a specific test problem and present numerical results that illustrate the sharpness of our estimates and also indicate directions for possible improvements.

2 Analysis on a single pipe

Let us start with describing in more detail the model problem under investigation. The pipe shall be represented by the unit interval and we consider

[TABLE]

We assume that $0<\underline{a}\leq a(x)\leq\overline{a}$ and that the pressure at the boundary is given by

[TABLE]

For ease of presentation $g_{0}$, $g_{1}$ are assumed to be independent of time here. Other boundary conditions could be considered with obvious modifications. From standard results of semigroup theory, one can easily deduce the following.

Lemma 1

Let $p_{0},m_{0}\in H^{1}(0,1)$ be given with $p_{0}(0)=g_{0}$ and $p_{1}(1)=g_{1}$. Then for any $\epsilon>0$ problem (3)–(5) has a unique classical solution

[TABLE]

satisfying initial conditions $p^{\epsilon}(x,0)=p_{0}(x)$ and $m^{\epsilon}(x,0)=m_{0}(x)$ for all $x\in(0,1)$.

The parabolic problem (3)–(5) with $\epsilon=0$ also has a unique solution

[TABLE]

satisfying the initial condition $p^{0}(x,0)=p_{0}(x)$ for all $x\in(0,1)$.

Note that only one single initial condition is required in the parabolic limit. By elementary arguments one can verify that the corresponding stationary problem

[TABLE]

is independent of $\epsilon$ and has a unique solution $(\bar{p},\bar{m})\in H^{1}(0,1)\times H^{1}(0,1)$ as well. Using standard energy arguments and the linearity of the time dependent and of the stationary problem, one can then establish the following assertions.

Lemma 2

Let $(p^{\epsilon},m^{\epsilon})$ and $(\bar{p},\bar{m})$ denote solutions of (3)–(5) and (6)–(8), respectively. Then for any $\epsilon\geq 0$ and any $t\geq 0$, there holds

[TABLE]

For $\epsilon>0$, one can additionally bound the time derivatives of $(p^{\epsilon},m^{\epsilon})$ by

[TABLE]

Here and below, $\|\cdot\|$ and $(\cdot,\cdot)$ denote the norm and the scalar product on $L^{2}(0,1)$. In addition, the functions $p^{\epsilon}$, $m^{\epsilon}$ are understood as functions of time with values in Hilbert spaces. The fact that the second estimate degenerates as $\epsilon\to 0$ resembles the fact that the second initial condition becomes superfluous in the parabolic limit.

Proof

∎Due to linearity of the problem, we may assume without loss of generality that $g_{0}=g_{1}=0$ and hence $\bar{p}\equiv\bar{m}\equiv 0$. From (3)–(4) we then get

[TABLE]

Using integration-by-parts for the second term in the last line, the homogeneous boundary conditions for $p^{\epsilon}$, and the lower bound for the parameter $a$, we get

[TABLE]

The first estimate now follows by integration with respect to time. Next assume that $(p^{\epsilon},m^{\epsilon})\in C^{2}(\mathbb{R}^{+};L^{2}(0,1)\times L^{2}(0,1))$. Then by formal differentiation of the problem one can see, that the time derivative $(\partial_{t}p^{\epsilon},\partial_{t}m^{\epsilon})$ also solves (3)–(5) with homogeneous boundary conditions. The previous estimate thus yields

[TABLE]

The differential equations (3) and (4) can be used to replace the terms on the right hand side which proves the second estimate for the case of smooth solutions. The general case finally follows by a density argument. ∎

A combination of these energy estimates allows us to provide a precise formulation and to prove the first assertion about solutions of the continuous problem.

Theorem 2.1

Let $\epsilon>0$ and let $(p^{\epsilon},m^{\epsilon})$ and $(p^{0},m^{0})$ denote the unique solutions of problem (3)–(5) with initial values $p^{\epsilon}(0)=p^{0}(0)=p^{0}$ and $m^{\epsilon}(0)=m_{0}$. Then

[TABLE]

Proof

∎Let $r^{\epsilon}=p^{\epsilon}-p^{0}$ and $w^{\epsilon}=m^{\epsilon}-m^{0}$ denote the differences between the solutions of the hyperbolic and the parabolic problem. Then by linearity of the equations, one can deduce that $r^{\epsilon}=0$ at the boundary and that

[TABLE]

Applying similar arguments as in the proof of the previous lemma then leads to

[TABLE]

Multiplication by two and integration with respect to time further yields

[TABLE]

Since $p^{\epsilon}$ and $p^{0}$ satisfy the same initial conditions, we have $r^{\epsilon}(0)=0$, and the remaining integral on the right hand side can be estimated by Lemma 2. ∎

The estimates of Lemma 2 provide uniform bounds for the distance to steady state. A refined analysis reveals that in fact exponential convergence takes place.

Theorem 2.2

Let $(p^{\epsilon},m^{\epsilon})$ denote a solution of (3)–(5) for some $0\leq\epsilon\leq 1$. Further let $(\bar{p},\bar{m})$ be the unique solution of the corresponding stationary problem. Then

[TABLE]

which holds for all $0\leq s\leq t$ and with some constants $C,\gamma>0$ independent of $\epsilon$.

Proof

∎Set $\tau=t/\epsilon$ and $\sigma=s/\epsilon$ and define $\pi^{\epsilon}(\tau)=p^{\epsilon}(t)$ and $\mu^{\epsilon}(\tau)=\epsilon m^{\epsilon}(t)$. Then by elementary calculations, one can see that

[TABLE]

The exponential convergence for this problem has been established in EggerKugler15 and a direct application of Theorem 3.3 in EggerKugler15 yields

[TABLE]

Using $\tau=t/\epsilon$ and $\sigma=s/\epsilon$ and the definition of $\pi^{\epsilon}$ and $\mu^{\epsilon}$ then directly yields the estimate for $\epsilon>0$. The result for $\epsilon=0$ follows directly but also from the uniformity of those for $\epsilon>0$ and the convergence to the parabolic limit. ∎

3 A mixed finite element Runge-Kutta scheme

For the discretization of problem (3)–(5), we now consider a mixed finite element method in space and the implicit Euler method in time. More general Galerkin and time-integration schemes could be analyzed in a similar manner. Let $T_{h}=\{e\}$ denote a uniform mesh of the interval $(0,1)$ into elements $e$ of size $h$ and denote by

[TABLE]

the spaces of piecewise constant and piecewise linear and continuous functions, respectively. Furthermore, let $\tau>0$ be the time step size, define $t^{k}=k\tau$, and denote by $\bar{\partial}_{\tau}u(t^{k})=\frac{1}{\tau}[u(t^{k})-u(t^{k-1})]$ the backward difference quotient. We then consider

Problem 1

Let $p^{\epsilon}_{h,\tau}(0)$ and $m^{\epsilon}_{h,\tau}(0)$ be the $L^{2}$ projections of the initial data onto the finite element spaces. For $k\geq 1$ find $(p_{h,\tau}^{\epsilon}(t^{k}),m_{h,\tau}^{\epsilon}(t^{k}))\in Q_{h}\times V_{h}$, such that

[TABLE]

holds for all test functions $q_{h}\in Q_{h}$ and all $v_{h}\in V_{h}$.

Recall that $(\cdot,\cdot)$ denotes the scalar product of $L^{2}(0,1)$. Existence of a unique discrete solution $(p^{\epsilon}_{h,\tau},m^{\epsilon}_{h,\tau})$ to Problem 1 and of a unique solution $(\bar{p}_{h},\bar{m}_{h})$ of the corresponding stationary problem can be deduced from the results in EggerKugler15 .

Lemma 3

For any $\epsilon\geq 0$, Problem 1 admits a unique solution $(p_{h,\tau}^{\epsilon},m_{h,\tau}^{\epsilon})$ and

[TABLE]

for all $k\geq 0$, where $(\bar{p}_{h},\bar{m}_{h})\in Q_{h}\times V_{h}$ denotes the unique solution of the corresponding stationary problem. For $\epsilon>0$, we additionally have

[TABLE]

with constant $C$ that is independent of $\epsilon$ and the discretization parameters $h$ and $\tau$.

Proof

∎Without loss of generality, we may set $g_{0}=g_{1}=0$ and hence $\bar{p}_{h}\equiv\bar{m}_{h}\equiv 0$. For ease of notation, let us abbreviate $p^{k}:=p_{h,\tau}^{\epsilon}(t^{k})$ and $m^{k}:=m_{h,\tau}^{\epsilon}(t^{k})$. Then by elementary calculations, one can verify that

[TABLE]

Using the discrete problem and the lower bounds for the parameter, we thus obtain

[TABLE]

The first estimate now follows by recursion and by noting that $\|p^{0}\|\leq\|p_{0}\|$ and $\|m^{0}\|\leq\|m_{0}\|$, since the initial iterates were defined as $L^{2}$ orthogonal projections of the initial values onto the respective subspaces. By linearity of the problem, one can then deduce in a similar manner that

[TABLE]

Using the discrete problem for $k=1$, we further get

[TABLE]

Using Young’s inequality, the bounds for the parameter $a$, and the stability of the $L^{2}$ projection in the $H^{1}$ norm, we may conclude that

[TABLE]

which together with the energy estimate from above completes the proof. ∎

Similarly as on the continuous level, a combination of the previous estimates now immediately allows to show convergence of the solutions $(p^{\epsilon}_{h,\tau},m^{\epsilon}_{h,\tau})$ of the discrete hyperbolic problem to that of the discrete parabolic problem when $\epsilon\to 0$.

Theorem 3.1

Let $(p^{\epsilon}_{h,\tau},m^{\epsilon}_{h,\tau})$ and $(p^{0}_{h,\tau},m^{0}_{h,\tau})$ denote solutions of Problem 1 for $\epsilon>0$ and $\epsilon=0$, respectively. Further assume that $p^{\epsilon}_{h,\tau}(0)=p^{0}_{h,\tau}(0)$. Then

[TABLE]

with constant $C$ independent of $\epsilon$ and of the discretization parameters $h$ and $\tau$.

Proof

∎Define $r^{k}=p_{h,\tau}^{\epsilon}(t^{k})-p_{h,\tau}^{0}(t^{k})$ and $w^{k}=m_{h,\tau}^{\epsilon}(t^{k})-w^{0}_{h,\tau}(t^{k})$. Then by linearity of the discrete problem, one can see that

[TABLE]

for all $q_{h}\in Q_{h}$ and $v_{h}\in V_{h}$ and for all $k\geq 0$. Testing with $q_{h}=w^{k}$ and $v_{h}=m^{k}$ and proceeding similarly as in the previous lemmas leads to the energy estimate

[TABLE]

The assertion now follows by noting that $r^{0}\equiv 0$ and application of the second estimate of the previous lemma to estimate the last term in this expression. ∎

Similarly as on the continuous level, one can again prove uniform exponential convergence of discrete solutions to steady states.

Theorem 3.2

Let $(p^{\epsilon}_{h,\tau},m^{\epsilon}_{h,\tau})$ denote a solution of Problem 1 and let $(\bar{p}_{h},\bar{m}_{h})$ let be the unique solution of the corresponding stationary problem. Then

[TABLE]

for all $0\leq j\leq k$ with constants $C,\gamma>0$ that are independent of $\epsilon$, $h$, and $\tau$.

Proof

∎Using a rescaling like in the proof of Theorem 2.2, the result for $\epsilon>0$ can be deduced directly from Theorem 7.4 in EggerKugler15 . The estimate for $\epsilon=0$ follows from the uniformity of the estimates and convergence to the parabolic limit.∎

4 Extension to pipe networks

The results of the previous sections can be extended to the following class of hyperbolic problems on networks: Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a finite directed graph representing the topology of the network. On every single pipe $e$, the dynamics shall again be described by the linear damped hyperbolic system

[TABLE]

At any junction $v$ of several pipes $e\in\mathcal{E}(v)$ of the network, we require that

[TABLE]

Here $n_{e}(v)$ takes the value minus or plus one, depending on whether the pipe $e$ start or ends at the junction $v$. At the boundary vertices $v$ of the network, we require

[TABLE]

Using the arguments developed in EggerKugler16 , all results stated in Theorem 2.1–3.2 hold verbatim also for the system (9)–(13). Details are left to the interested reader.

5 Numerical validation

We now illustrate our theoretical results by considering in detail a particular model problem. For constant damping parameter $a\equiv 1$, initial data $p_{0}=sin(\pi x)$, $m_{0}\equiv 0$, and boundary values $g_{0}=g_{1}\equiv 0$, the solution of problem (3)–(5) is given by

[TABLE]

and

[TABLE]

with parameter $s(\epsilon)=\sqrt{1-4\pi^{2}\epsilon^{2}}$. By Taylor expansion w.r.t. $\epsilon$, we deduce that

[TABLE]

and

[TABLE]

For $\epsilon=0$, we simply obtain $p^{0}(x,t)=e^{-\pi^{2}t}\sin(\pi x)$ and $m^{0}(x,t)=\pi e^{-\pi^{2}t}\cos(\pi x)$ and the steady state for this problem is given by $\bar{p},\bar{m}\equiv 0$.

From the explicit solution formulas, one can then immediately see that exponential convergence towards the steady state takes place with $t\to\infty$ for all $0\leq\epsilon\leq 1$ with a rate that is independent of $\epsilon$ which was the assertion of Theorem 2.2. In Table 1, we depict numerical results obtained with the numerical scheme discussed in Section 3.

As predicted by Theorem 3.2, the exponential convergence towards steady state with $t\to\infty$ is uniform in $\epsilon$ also for the discrete schemes. Mesh independence of the exponential decay rate was already demonstrated in EggerKugler16 .

Let us next have a closer look on the convergence to the parabolic limit. Using the analytical solution formulas and Taylor expansion w.r.t. $\epsilon$, one can deduce that

[TABLE]

and

[TABLE]

This shows that $\|p^{\epsilon}-p^{0}\|^{2}=\mathcal{O}(\epsilon^{4})$ and $\int_{0}^{t}\|m^{\epsilon}-m^{0}\|^{2}=\mathcal{O}(\epsilon^{2})$ which yields exactly the asymptotic behavior predicted in Theorem 2.1. In Table 2, we display the corresponding results obtained with the proposed discretization scheme.

Also here we can exactly observe the convergence rate predicted by Theorem 3.1. Note that the second term in the error measure is strictly increasing w.r.t. time, which together with the exponential convergence to steady states explains that the error is almost independent of $t$ here.

In Table 3, we report about further numerical tests to illustrate the independence of the results on the discretization parameters.

Again, the observations are in perfect agreement with the theoretical predictions made in Theorem 3.1.

Let us finally note that the previous formulas reveal that the error between the solutions of the hyperbolic and the parabolic problem actually behaves like

[TABLE]

This shows that the estimate of Theorem 3.1 is dominated by the error in the mass flux within the initial layer $0\leq t\preceq\epsilon$ which again resembles the fact that the second initial condition gets superfluous in the parabolic limit. This behavior can also be observed for the numerical approximations obtained with the method discussed in Section 3. A theoretical explanation of this fact would require a refined analysis which is left for future research.

Acknowledgements.

The authors would like to gratefully acknowledge financial support by the German Research Foundation (DFG) via grants IRTG 1529, GSC 233, and TRR 154.