Time discretization of nonlinear hyperbolic systems on networks

In view of gas networks, the simulation of hyperbolic systems on networks has recently caused large interest. We consider the case with a nonlinear damping term and a small parameter ϵ such that – in the linear case – the system turns parabolic for ϵ =0. Based on this property and an expansion in ϵ, we derive a novel time integration scheme.


Introduction
We consider the propagation of pressure waves in a network of pipes which leads to a coupled system of hyperbolic partial differential equations. Each pipe is modeled by an edge in a graph on which we consider the one-dimensional wave systeṁ εṁ(t, x) + ∂ x p(t, x) + δ m(t, x)|m(t, x)| p(t, x) = 0.
Here, p and m model the pressure and mass flux on a single edge e. For the parameter ε we assume in this paper that 0 ≤ ε 1. In the considered case of a gas network this equals the product of the adiabatic coefficient and the square of the Mach number and is of order 10 −3 , cf. [1]. Note that in the here considered model we neglect differences in temperature and kinetic energy, cf. [2].
In the corresponding linear model, the nonlinear damping term in (1b) is often replaced by δm(t, x). In the network case, we need to add continuity conditions on the pressure and balance laws for the mass flow in each junction. In total, this then leads to a partial differential-algebraic equation (PDAE), cf. [3,4]. In this paper, we translate the expansion strategy from [4] to the nonlinear model given in (1). The expansion is motivated by the fact that -in the linear case -the system is hyperbolic for ε > 0 but parabolic in the limit case ε = 0.

Expansion and discretization
Due to the assumption ε 1, we consider an expansion in ε, meaning that we write where p 0 , m 0 serve as a first-order andp := p 0 + εp 1 ,m := m 0 + εm 1 as second-order approximations. Including the continuity and boundary conditions of the pressure in form of an explicit constraint Bp = g dir and the flow condition for the mass flux implicitly by Cr, we obtain the PDAE systemṡ Therein, K corresponds to the differential operator ∂ x and the nonlinear terms are given by Note that we have used the fact that f (x) = x|x| = sign(x)x 2 has a (classic) derivative f (x) = 2|x|. Further, we emphasize that the operator equations correspond to the weak formulation of the systems and that we consider p to be piecewise in H 1 on the edges of the graph whereas m is piecewise in L 2 .
For the numerical approximation of p and m we now consider the interplay of the ε-expansion and a time-discretization of the respective PDAE systems. In the linear case, one benefits here from the parabolic behavior of the limit equations [4].  For the spatial discretization we consider P 1 finite elements for the pressure and P 0 finite elements for the mass flux, both with mesh size h. This then leads to the semi-discrete system where M p , M m , K, and B denote the standard mass, stiffness, and boundary matrices. The term C T r includes the demand of the consumer on the right-hand side, whereas g dir includes again the boundary data.

Numerical example
To illustrate the performance of the expansion-based time discretization, we consider the case of a single pipe of length 1 with δ = 0.5 and consistent initial data, i.e.,ṁ h (0) = 0. As boundary conditions we set p = 1 on the left-hand side of the pipe and m = e −t on the right. Further, we restrict ourselves to the implicit Euler discretization in time but emphasize that this may be replaced by Runge-Kutta schemes, cf. [5]. Figure 1 shows the convergence history of the errors of p 0 ,p, m 0 , and m. As in the linear case, we can observe first-order convergence for p 0 , m 0 and second-order convergence forp,m. In our second numerical experiment, we consider the convergence order of p − p 0 , denoted by α, as a function of the spatial mesh size h. We choose non-consistent initial data, i.e., m(0) = m 0 (0), and Dirichlet boundary conditions for the pressure on both ends of the pipe. To compute α we fitted the exponent by a least square problem for several values of ε. In Table 1 one can observe that the exponent is slowly decreasing to around 0.767. In particular, the asymptotic limit is far off one, which equals the exponent in the case of consistent initial data. This behavior was also observed in the linear case and transfers to the here considered nonlinear one.