Finite Volume approximations of the Euler system with variable congestion

TL;DR

This paper develops an asymptotic preserving finite volume scheme for simulating the Euler system with variable congestion, effectively modeling crowd dynamics with singular pressure effects in multiple dimensions.

Contribution

It introduces a novel AP scheme for the Euler system with congestion, including a second order accurate version and validation in 1D and 2D scenarios.

Findings

The scheme accurately captures crowd behavior in simulations.

The second order scheme improves precision over first order.

Numerical results demonstrate realistic crowd dynamics modeling.

Abstract

We are interested in the numerical simulations of the Euler system with variable congestion encoded by a singular pressure. This model describes for instance the macroscopic motion of a crowd with individual congestion preferences. We propose an asymptotic preserving (AP) scheme based on a conservative formulation of the system in terms of density, momentum and density fraction. A second order accuracy version of the scheme is also presented. We validate the scheme on one-dimensional test-cases and extended here to higher order accuracy. We finally carry out two dimensional numerical simulations and show that the model exhibit typical crowd dynamics.

Figures

Tables

Table 1. Table 1: L 1 subscript 𝐿 1 L_{1} error between the numerical solutions to Riemann problem ( 24 ) and exact solution at time t = 0.1 𝑡 0.1 t=0.1 . Numercial solution computed using the ( ϱ , 𝒒 , Z ) italic-ϱ 𝒒 𝑍 (\varrho,\boldsymbol{q},Z) -method. Numerical parameters: Δ ​ x = 1 × 10 − 3 Δ 𝑥 1 superscript 10 3 \Delta x=1\times 10^{-3} , Δ ​ t = 0.1 ​ Δ ​ x Δ 𝑡 0.1 Δ 𝑥 \Delta t=0.1\,\Delta x , α = 2 𝛼 2 \alpha=2 , γ = 2 𝛾 2 \gamma=2 .

		$\varrho$	$q$	$Z$	${\varrho }^{\ast }$
$\epsilon ={10}^{-2}$	order 2 in $x$	$8.66\times {10}^{-4}$	$1.28\times {10}^{-3}$	$3.03\times {10}^{-4}$	$5.70\times {10}^{-4}$
	order 2	$1.17\times {10}^{-3}$	$3.52\times {10}^{-3}$	$5.89\times {10}^{-4}$	$5.77\times {10}^{-4}$
$\epsilon ={10}^{-4}$	order 2 in $x$	$9.75\times {10}^{-4}$	$2.11\times {10}^{-3}$	$3.70\times {10}^{-4}$	$5.71\times {10}^{-4}$
	order 2	$9.89\times {10}^{-4}$	$3.04\times {10}^{-3}$	$3.84\times {10}^{-4}$	$5.77\times {10}^{-4}$

Equations

$${\partial }_t\varrho +\nabla \cdot (\varrho \mathbf{v})=0,$$

$${\partial }_t(\varrho \mathbf{v})+\nabla \cdot (\varrho \mathbf{v}\otimes \mathbf{v})+\nabla \pi +\nabla p\left(\frac{\varrho }{{\varrho }^{\ast }}\right)=0,$$

$${\partial }_t{\varrho }^{\ast }+\mathbf{v}\cdot \nabla {\varrho }^{\ast }=0,$$

$$0\le \varrho \le {\varrho }^{\ast },$$

$$\pi ({\varrho }^{\ast }-\varrho )=0,\pi \ge 0,$$

$$\begin{array}{r}\varrho (0,x)={\varrho }_0(x)\ge 0,\mathbf{v}(0,x)={\mathbf{v}}_0(x),{\varrho }^{\ast }(0,x)={\varrho }_0^{\ast }(x),{\varrho }_0<{\varrho }_0^{\ast },\end{array}$$

$$\begin{array}{r}p\left(\frac{\varrho }{{\varrho }^{\ast }}\right)={\left(\frac{\varrho }{{\varrho }^{\ast }}\right)}^{\gamma },\gamma >1,\end{array}$$

$$\begin{array}{r}\nabla \cdot \mathbf{v}=0\text{in}\{\varrho ={\varrho }^{\ast }\},\end{array}$$

$$\begin{array}{r}{\pi }_{\epsilon }\left(\frac{\varrho }{{\varrho }^{\ast }}\right)=\epsilon {\left(\frac{\frac{\varrho }{{\varrho }^{\ast }}}{1-\frac{\varrho }{{\varrho }^{\ast }}}\right)}^{\alpha },\alpha >0.\end{array}$$

$${\partial }_t\varrho +\nabla \cdot (\varrho \mathbf{v})=0,$$

$${\partial }_t(\varrho \mathbf{v})+\nabla \cdot (\varrho \mathbf{v}\otimes \mathbf{v})+\nabla \pi =0,$$

$$0\le \varrho \le 1$$

$$\pi (\varrho -1)=0,\pi \ge 0.$$

$$\begin{array}{r}{\partial }_t\varrho +\nabla \cdot \mathbf{q}=0,\end{array}$$

$$\begin{array}{r}{\partial }_t\mathbf{q}+\nabla \cdot \left(\frac{\mathbf{q}\otimes \mathbf{q}}{\varrho }+{\pi }_{\epsilon }(Z)\mathbf{I}+p(Z)\mathbf{I}\right)=0,\end{array}$$

$$\begin{array}{r}{\partial }_{t}Z+\nabla \cdot \left(Z\frac{\mathbf{q}}{\varrho }\right)=0,\end{array}$$

$$\begin{array}{r}\varrho (0,x)={\varrho }_0(x),\mathbf{q}(0,x)={\mathbf{q}}_0(x),Z(0,x)=Z_0(x),\end{array}$$

$$\begin{array}{rl}& {\lambda }_1^{\epsilon }(\varrho ,q_1,Z)=\frac{q_1}{\varrho }-\sqrt{\frac{Z}{\varrho }{p}_{\epsilon }^{\prime }(Z)},\\ & {\lambda }_2^{\epsilon }(\varrho ,q_1,Z)=\frac{q_1}{\varrho },\\ & {\lambda }_3^{\epsilon }(\varrho ,q_1,Z)=\frac{q_1}{\varrho }+\sqrt{\frac{Z}{\varrho }{p}_{\epsilon }^{\prime }(Z)},\end{array}$$

$$\frac{{\varrho }^{n+1}-{\varrho }^n}{\Delta t}+{\nabla }_x\cdot {\mathbf{q}}^{n+1}=0,$$

$$\frac{{\mathbf{q}}^{n+1}-{\mathbf{q}}^n}{\Delta t}+{\nabla }_x\cdot \left(\frac{{\mathbf{q}}^n\otimes {\mathbf{q}}^n}{{\varrho }^n}+p(Z^n)\mathbf{I}\right)+{\nabla }_x({\pi }_{\epsilon }(Z^{n+1}))=0,$$

$$\frac{Z^{n+1}-Z^n}{\Delta t}+{\nabla }_x\cdot \left(Z^n\frac{{\mathbf{q}}^{n+1}}{{\varrho }^n}\right)=0.$$

$$\frac{Z^{n+1}-Z^n}{\Delta t}+{\nabla }_x\cdot \left(Z^n\frac{{\mathbf{q}}^n}{{\varrho }^n}\right)$$

$$-\Delta t{\nabla }_x\cdot \left(\frac{Z^n}{{\varrho }^n}{\nabla }_x\cdot \left(\frac{{\mathbf{q}}^n\otimes {\mathbf{q}}^n}{{\varrho }^n}+p(Z^n)\mathbf{I}\right)+\frac{Z^n}{{\varrho }^n}{\nabla }_x({\pi }_{\epsilon }(Z^{n+1}))\right)=0,$$

$$Z^{n+1}-\Delta t^2{\nabla }_x\cdot \left(\frac{Z^n}{{\varrho }^n}{\nabla }_x\left({\pi }_{\epsilon }(Z^{n+1})\right)\right)=\varphi ({\varrho }^n,q^n,Z^n),$$

$$\varphi ({\varrho }^n,q^n,Z^n)$$

$$=Z^n+\Delta t^2{\nabla }_x\cdot \left(\frac{Z^n}{{\varrho }^n}{\nabla }_x\cdot \left(\frac{{\mathbf{q}}^n\otimes {\mathbf{q}}^n}{{\varrho }^n}+p(Z^n)\mathbf{I}\right)\right)-\Delta t{\nabla }_x\cdot \left(Z^n\frac{{\mathbf{q}}^n}{{\varrho }^n}\right).$$

$$\frac{{\varrho }_i^{n+1}-{\varrho }_i^n}{\Delta t}+\frac{1}{\Delta x}(F_{i+\frac{1}{2}}^{n+1}-F_{i-\frac{1}{2}}^{n+1})=0,$$

$$\frac{q_i^{n+1}-q_i^n}{\Delta t}+\frac{1}{\Delta x}(G_{i+\frac{1}{2}}^n-G_{i-\frac{1}{2}}^n)+\frac{{\pi }_{\epsilon }(Z_{i+1}^{n+1})-{\pi }_{\epsilon }(Z_{i-1}^{n+1})}{2\Delta x}=0,$$

$$\frac{Z_i^{n+1}-Z_i^n}{\Delta t}+\frac{1}{\Delta x}(H_{i+\frac{1}{2}}^{n+1}-H_{i-\frac{1}{2}}^{n+1})=0.$$

$$\displaystyle F_{i+\frac{1}{2}}^{n+1}=\frac{1}{2}\big{(}q^{n+1}_{i+1}+q^{n+1}_{i}\big{)}-(D_{\varrho})_{i+\frac{1}{2}}^{n},$$

$$\displaystyle G_{i+\frac{1}{2}}^{n}=\frac{1}{2}\Big{(}\frac{(q_{i+1}^{n})^{2}}{\varrho_{i+1}^{n}}+\frac{(q_{i}^{n})^{2}}{\varrho_{i}^{n}}+p(Z_{i+1}^{n})+p(Z_{i}^{n})\Big{)}-(D_{q})_{i+\frac{1}{2}}^{n},$$

$$\displaystyle H_{i+\frac{1}{2}}^{n+1}=\frac{1}{2}\Big{(}\frac{Z^{n}_{i+1}}{\varrho^{n}_{i+1}}q^{n+1}_{i+1}+\frac{Z^{n}_{i}}{\varrho^{n}_{i}}q^{n+1}_{i}\Big{)}-(D_{Z})_{i+\frac{1}{2}}^{n}.$$

$$(D_{w})_{i+\frac{1}{2}}^{n}=\frac{1}{2}c^{n}_{i+\frac{1}{2}}\big{(}w_{i+1}^{n}-w_{i}^{n}\big{)},$$

$$c^{n}_{i+\frac{1}{2}}=\max\left\{\left|\lambda_{k}^{0}\big{(}\varrho^{n}_{i+1},q^{n}_{i+1},Z^{n}_{i+1}\big{)}\right|,\left|\lambda_{k}^{0}\big{(}\varrho^{n}_{i},q^{n}_{i},Z^{n}_{i}\big{)}\right|,\quad k=1,2,3\right\},$$

$$Z_i^{n+1}-Z_i^n+\frac{\Delta t}{\Delta x}({\bar{H}}_{i+1/2}^n-{\bar{H}}_{i-1/2}^n)$$

$$\displaystyle-\frac{\Delta t^{2}}{\Delta x^{2}}\frac{1}{2}\Big{(}\frac{Z^{n}_{i+1}}{\varrho^{n}_{i+1}}(G_{i+\frac{3}{2}}^{n}-G_{i+\frac{1}{2}}^{n})-\frac{Z^{n}_{i-1,j}}{\varrho^{n}_{i-1}}(G_{i-\frac{1}{2}}^{n}-G_{i-\frac{3}{2}}^{n})\Big{)}$$

$$\displaystyle-\frac{\Delta t^{2}}{\Delta x^{2}}\frac{1}{2}\Big{(}\frac{Z^{n}_{i+1}}{\varrho^{n}_{i+1}}\big{(}\pi_{\varepsilon}(Z^{n+1}_{i+2})-\pi_{\varepsilon}(Z^{n+1}_{i})\big{)}-\frac{Z^{n}_{i-1}}{\varrho^{n}_{i-1}}\big{(}\pi_{\varepsilon}(Z^{n+1}_{i})-\pi_{\varepsilon}(Z^{n+1}_{i-2})\big{)}\Big{)}=0,$$

$$\bar{H}_{i+\frac{1}{2}}^{n}=\frac{1}{2}\Big{(}\frac{Z^{n}_{i+1}}{\varrho^{n}_{i+1}}q^{n}_{i+1}+\frac{Z^{n}_{i}}{\varrho^{n}_{i}}q^{n}_{i}\Big{)}-(D_{Z})_{i+\frac{1}{2}}^{n}.$$

$$Z_i^{n+1}(({\pi }_{\epsilon }{)}_i^{n+1})-\frac{\Delta t^2}{\Delta x^2}\frac{1}{2}$$

$$\displaystyle-\frac{Z^{n}_{i-1}}{\varrho^{n}_{i-1}}\big{[}(\pi_{\varepsilon})^{n+1}_{i}-(\pi_{\varepsilon})^{n+1}_{i-2}\big{]}\Big{)}=\phi(\varrho^{n},q^{n},Z^{n})_{i},$$

$$\varphi ({\varrho }^n,q^n,Z^n{)}_i=$$

Full text

Keywords: fluid model of crowd, Euler equations, free boundary, singular pressure, finite volume, Asymptotic-Preserving scheme

1. Department of Mathematics, Imperial College London,

London SW7 2AZ, United Kingdom.

2. Interdisciplinary Center for Scientific Computing, Heidelberg University,

Im Neunheimer Feld 205, D-69120 Heidelberg, Germany.

3. Institut de Recherche Mathématique Avancée, UMR 7501,

Université de Strasbourg et CNRS, 7 rue René Descartes,

67000 Strasbourg, France & Inria TONUS

1 Introduction

In this work we study two phase compressible/incompressible Euler system with variable congestion:

[TABLE]

with the initial data

[TABLE]

where the unknowns are: $\varrho=\varrho(t,x)$ – the mass density, $\boldsymbol{v}=\boldsymbol{v}(t,x)$ – the velocity, $\varrho^{*}=\varrho^{*}(t,x)$ – the congestion density, and $\pi$ – the congestion pressure. The barotropic pressure $p$ is an explicit function of the density fraction $\frac{\varrho}{\varrho^{*}}$

[TABLE]

and plays the role of the background pressure.

The congestion pressure $\pi$ appears only when the density $\varrho$ achieves its maximal value, the congestion density $\varrho^{*}$. Therefore $\varrho^{*}$ is sometimes referred to as the barrier or the threshold density. It was observed in [25], and then generalized in [15], that the restriction on the density (1d) is equivalent with the condition

[TABLE]

if only $\varrho,\boldsymbol{v},\varrho^{*}$ are sufficiently regular solutions of the continuity equation (1a) and the transport equation (1c). For that reason, system (1) can be seen as a free boundary problem for the interface between the compressible (uncongested) regime $\{\varrho<\varrho^{*}\}$ and the incompressible (congested) regime $\{\varrho=\varrho^{*}\}$.

The main purpose of this work is to analyze (1) numerically, i.e. to propose the numerical scheme capturing the phase transition. To this end we use the fact that (1) can be obtained as a limit when $\varepsilon\to 0$ of the compressible Euler system with the congestion pressure $\pi$ replaced by a singular approximation $\pi_{\varepsilon}$

[TABLE]

Note that for fixed $\varepsilon>0$, $\pi_{\varepsilon}\to\infty$ when $\varrho\to\varrho^{*}$. Therefore, at least formally, for $\varepsilon\to 0$, $\pi_{\varepsilon}$ converges to a measure supported on the set of singularity, i.e. $\{(x,t)\in\Omega\times(0,T):\varrho(x,t)=\varrho^{*}(x,t)\}$. The rigorous proof of this fact is an open problem, at least for the Euler type of systems. There have been, however, several results for a viscous version of the model, see [10] for the one-dimensional case, [31] for multi-dimensional domains and space-dependent congestion $\varrho^{*}(x)$ and [15] for the case of congestion density satisfying the transport equation (1c). An alternative approximation leading to a similar two-phase system was considered first by P.-L. Lions and N. Masmoudi [25], and more recently for the model of tumour growth [32]. The advantage of approximation (5) considered here lies in the fact that for each $\varepsilon$ fixed, the solutions to the approximate system stay in the physical regime, i.e. $\varrho\leq\varrho^{*}$. This feature is especially important for the numerical purposes, see for example [27] for further discussion on this subject.

System (1) is a generalization of the pressureless Euler system with the maximal density constraint

[TABLE]

introduced originally by Bouchut et al. [8], who also proposed the first numerical scheme based on an approach developed earlier for the pressureless systems, see for example [9], and the projection argument. The model was studied later on by Berthelin [3, 4] by passing to the limit in the so-called sticky-blocks dynamics, see also [35], and a very interesting recent paper [30] using the Lagrangian approach for the monotone rearrangement of the solution to prove the existence of solutions to (6) with additional memory effects.

The pressureless Euler equations with the density constraint were originally introduced in order to describe the motion of particles of finite size. Our model extends this concept by including the variance of the size of particles. In system (1) $\varrho^{*}$ is given initially and is transported along with the flow.

One can also think of $\varrho^{*}$ as a congestion preference of individuals moving in the crowd (cars, pedestrians), which is one of the factors determining their final trajectory and the speed of motion. The macroscopic modelling of crowd is one of possible approaches and it allows to determine the averaged quantities such as the density and the mean velocity rather than the precise position of an individual. One of the first models of this kind based on classical mechanics was introduced by Henderson [22]. More sophisticated model was introduced by Hughes [23] where the author considers the continuity equation equipped with a phenomenological constitutive relation between the velocity and the density. For a survey of the crowd models we refer the reader to [11, 1, 33, 28, 24] and to the review paper [2].

As far as the numerical methods are concerned, the macroscopic models of pedestrian flow with condition preventing the overcrowding were studied, for example in [34]. The influence of the maximal density constraint was investigated also in the context of vehicular traffic in [7, 5]. The strategy that we want to adapt in this paper, i.e. to use the singularities of the pressure similar to (5) has been developed in the past for a number of Euler-like systems for the traffic models [5, 7, 6], collective dynamics [14, 13], or granular flow [26, 29]. In our previous work [15], we have drafted the numerical scheme for system (1) in the one-dimensional case. We used a splitting algorithm at each time step that consists of three sub-steps. At first, the hyperbolic part is solved with the AP-preserving method presented in [14]. Next the diffusion is solved by means of cell-centered finite volume scheme, and the transport of the congested density is resolved with the upwind scheme.

The extension of this method to two-dimensions is one of the main results of the present paper. We also propose an alternative scheme using different formulation in terms of the conservative variables: the density $\varrho$, the momentum $\boldsymbol{q}=\varrho\boldsymbol{v}$, and the density fraction $Z=\frac{\varrho}{\varrho^{*}}$:

[TABLE]

where $Z_{0}=\frac{\varrho_{0}}{\varrho^{*}_{0}}$, and ${\boldsymbol{q}}_{0}=\varrho_{0}\boldsymbol{v}_{0}$. $\boldsymbol{I}$ denotes the identity tensor. This is a stricly hyperbolic system whose wave speeds in the $x_{1}$-direction are given by:

[TABLE]

where $p_{\varepsilon}=p+\pi_{\varepsilon}$, and $q_{1}$ denotes the component of $\boldsymbol{q}$ in the $x_{1}$ direction. Consequently, in region where the density $\varrho$ is closely congested, i.e. $Z$ is close to 1, the characteristic speeds of the system are extremely large. This corresponds to the nearly incompressible dynamics.

The paper is organised as follows. In Section 2 we present our numerical schemes using the two formulations (1) and (7). They are referred to as $(\varrho,\boldsymbol{q})$-method/SL and $(\varrho,\boldsymbol{q},Z)$-method, respectively. In Section 2.1 we describe the first-order semi-discretization in time and the full discretization for the $(\varrho,\boldsymbol{q},Z)$-method. Then, in Section 2.2, we discuss the second order scheme for the $(\varrho,\boldsymbol{q},Z)$-method. At last, in Section 2.3 we present the $(\varrho,\boldsymbol{q})$-method/SL for the system written in terms of the physical variables (1). Section 3 is devoted to validation of the schemes on the Riemann problem whose solutions are described in A. Finally, in Section 4 we discuss the two-dimensional numerical results: in Section 4.1 we present how these schemes work for three different initial congestion densities, and in Section 4.2 we present an application of $(\varrho,\boldsymbol{q})$-method/SL to model crowd behaviour in the evacuation scenario.

2 Numerical schemes

In this section, we first introduce a numerical scheme based on system (7) using the conservative variables. In order to use large time steps not restricted by too drastic CFL condition, implicit-explicit (IMEX) type methods need to be designed. The scheme can be solved through the following steps: first an elliptic equation on the density fraction $Z$ is solved, and then we update $\boldsymbol{q}$ and $\varrho$, respectively.

Such scheme is compared with an extension of the method introduced in [15], where the congestion density is advected separately from the update of $\varrho$ and $\boldsymbol{q}$. For the sake of completeness, a description of the scheme is given in Section 2.3.

2.1 The first order $(\varrho,\boldsymbol{q},Z)$-method

Discretization in time

We adopt the previous work [14] to introduce a method treating implicitly the stiff congestion pressure $\pi_{\varepsilon}(Z)$. We consider a constant time step $\Delta t>0$ and $\varrho^{n}$, $\boldsymbol{q}^{n}$, $Z^{n}$, $\varrho^{\ast\,n}$ denote the approximate solution at time $t^{n}=n\Delta t$, $\forall n\in\mathbb{N}$. We thus consider the following semi-implicit time discretization:

[TABLE]

Note that in the flux term in equation (9c), the momentum is taken implicitly. Inserting (9b) into (9c), we obtain:

[TABLE]

This is an elliptic equation on the unknown $Z^{n+1}$, that can be written as:

[TABLE]

where

[TABLE]

The $n$-th time step of the scheme is decomposed into three parts: first get $Z^{n+1}$ when solving (10), then compute $\boldsymbol{q}^{n+1}$ thanks to (9b) and then $\varrho^{n+1}$ from (9a).

Discretization in space

We only derive the fully discrete scheme in the one-dimensional case; the two-dimensional formula are given in B. We consider the computational domain $[0,1]$ and a spatial space step $\Delta x=1/N_{x}>0$, with $N_{x}\in\mathbb{N}$: the mesh points are thus $x_{i}=i\Delta x$, $\forall i\in\left\{0,\ldots,N_{x}\right\}$. Let $\varrho^{n}_{i}$, $\boldsymbol{q}^{n}_{i}$, $Z^{n}_{i}$, $\varrho^{\ast\,n}_{i}$ denote the approximate solution at time $t^{n}$ on mesh cell $[x_{i},x_{i+1}]$. The spatial discretization have to capture correctly the entropic solutions of the hyperbolic system. To derive the fully discrete scheme, we thus make the same algebra on the following fully discrete system:

[TABLE]

where the stiff pressure is discretized by the centered finite difference and the numerical fluxes $F^{n+1}$, $G^{n}$, $H^{n+1}$ (we denote implicit-explicit fluxes by current timestep $n+1$ and fully explicit fluxes by previous timestep $n$) are splitted into centered part and the upwinded part:

[TABLE]

The upwinded parts are given explicitly. They can be given by the diagonal Rusanov (or local Lax-Friedrichs) upwindings:

[TABLE]

for any conserved quantities $w$, where $c^{n}_{i+\frac{1}{2}}$ is the maximal characteristic speed (in absolute value):

[TABLE]

where $\lambda_{k}^{0}$ are given by eq. (8) with $\varepsilon=0$ (no congestion pressure). These correspond to the eigenvalues of the hyperbolic system taken explicitly in (9). One could also consider less diffusive numerical fluxes like the Polynomial upwind scheme [17].

Like in the semi-discrete case, we now obtain the fully discrete elliptic equation on $Z$ by replacing the implicit momentum terms appearing in the flux $H$ (14) by their expressions given by the momentum equation (11b). We get:

[TABLE]

where $\bar{H}^{n}$ denotes the same expression as (14) where all quantities are taken explicitly:

[TABLE]

In practice, in order to preserve the constraint $Z\leqslant 1$, this elliptic equation is solved with respect to the congestion pressure variable $\pi_{\varepsilon}$:

[TABLE]

where the right-hand side is given by:

[TABLE]

and $Z(\pi_{\varepsilon})$ is the inverse function of $\pi_{\varepsilon}(Z)$. This equation is supplemented by periodic or Dirichlet boundary conditions and the non-linear equation is solved using the Newton iterations.

The $(n\!+\!1)$-th time step of the algorithm thus consists in getting $Z^{n+1}$ by solving (17)-(18) and then obtaining $q^{n+1}$ from (11b) and $\varrho^{n+1}$ from (11a).

Since the singular pressure $\pi_{\varepsilon}$ is treated implicitly, the scheme remains stable even for small $\varepsilon$. The stability condition only depends on the wave speeds of the explicit part of the scheme, that is under the Courant-Friedrichs-Levy (CFL) condition:

[TABLE]

where $\lambda_{j}^{0}$, given by eq. (8), denotes the eigenvalues of the hyperbolic system with no congestion pressure ($\varepsilon=0$). The scheme is asymptotically stable with respect to $\varepsilon$.

2.2 The second order $(\varrho,\boldsymbol{q},Z)$-method

Discretization in time

The second-order discretization in time is based on the combined Runge-Kutta 2 / Crank-Nicolson (RK2CN) method as described in [12]: it consists into replacing Euler explicit by Runge-Kutta 2 solver and Euler Implicit by Crank-Nicolson solver in semi-discretization (9). We here only detail the semi-discretized scheme. However, to be unambiguous, we will denote by $\mathcal{D}_{\varrho}$, $\mathcal{D}_{q}$, and $\mathcal{D}_{Z}$ the numerical diffusion terms resulting from the upwinding terms and the divergence operators will be replaced by centered fluxes. We thus consider the following scheme:

First step (half time step): get $\varrho^{n+1/2}$, $\boldsymbol{q}^{n+1/2}$ and $Z^{n+1/2}$ from

[TABLE]

Second step (full time step): get $\varrho^{n+1}$, $\boldsymbol{q}^{n+1}$ and $Z^{n+1}$ from

[TABLE]

Like in the first-oder scheme, equations (20b)-(20c) and (21b)-(21c) result in elliptic equations for $Z$. In practice, the scheme may fail capturing discontinuities, in particular when small values of $\varepsilon$ are concerned. Indeed, the semi-implicit pressure \big{(}\pi_{\varepsilon}(Z^{n})+\pi_{\varepsilon}(Z^{n+1})\big{)}/2 in (21b) is constrained to be larger than $\pi_{\varepsilon}(Z^{n})/2$ preventing from having large discontinuities in pressure. One way to overcome this difficulty is to dynamically replace this semi-implicit pressure by an implicit pressure $\pi_{\varepsilon}(Z^{n+1})$ as soon as the non-linear solver of the elliptic equation detects a pressure lower than half the explicit one.

Discretization in space

To get second order accuracy in space, we consider a MUSCL strategy. For any conserved quantity $v$, it consists in introducing at each mesh interface left and right values $w_{L}$ and $w_{R}$:

[TABLE]

where the minmod function is defined as:

[TABLE]

Then all explicit terms in fluxes (12)-(13)-(14) depend on $(\varrho_{i,R}^{n},q_{i,R}^{n},Z_{i,R}^{n})$ and $(\varrho_{i+1,L}^{n},q_{i+1,L}^{n},Z_{i+1,L}^{n})$ instead of $(\varrho_{i}^{n},q_{i}^{n},Z_{i}^{n})$ and $(\varrho_{i+1}^{n},q_{i+1}^{n},Z_{i+1}^{n})$. Implicit terms are unchanged in order to be able to get the elliptic equation.

2.3 Congested Euler/Semi-Lagrangian scheme ($(\varrho,{\boldsymbol{q}})$-method/SL)

Discretization in time

We consider a scheme based on the non-conservative form (1) of the congestion transport. This idea was proposed in [14] in the context of constant congestion and in [15] in the context of variable congestion. The time-discretization reads:

[TABLE]

Inserting(22b) into (22a) results in

[TABLE]

This is an elliptic equation on the density $\varrho^{n+1}$. The $n$-th time step of the scheme is decomposed into three parts: first get $\varrho^{n+1}$ when solving (2.3), then compute $\boldsymbol{q}^{n+1}$ thanks to (22b) and then ${\varrho^{*}}^{n+1}$ from (22c).

Discretization in space

Like for the previous schemes, we restrict the description to the one-dimensional case. Finite volume discretization is used for the spatial discretization of (22a)-(22b) as in section 2.1, see also [14]. A semi-Lagrangian method is used to solve (22c) and thus update the congestion density $\varrho^{\ast}$. The congestion density $\varrho^{\ast n+1}_{i}$ at node $x_{i}$ and time $t^{n+1}$ is computed as follows: first we integrate back the characteristic line over one time step and then we interpolate the maximal density $\varrho^{\ast n}$ at that point. Using Euler scheme for the first step, we obtain:

[TABLE]

where $\Pi\varrho^{\ast n}$ is an interpolation function built from the points $(x_{i},\varrho^{\ast n}_{i})$. We here perform a Lagrange interpolation on the $2r+2$ neighboring points:

[TABLE]

resulting in $2r+1$-th spatial accuracy. First $(r=0)$ and third $(r=1)$ order in space semi-Lagrangian scheme will be used. For more details, we refer to [18].

The second order scheme

Extension of the full scheme to second order accuracy in space is made using the MUSCL strategy for the finite volume fluxes. Extension to second order accuracy in time requires a Crank-Nicolson/Runge Kutta 2 method for $(\varrho,q)$ and a second order in time integration of the characteric line for the semi-Lagrangian scheme (with for instance Taylor expansion) combined to a Strang splitting, see C.

3 One dimensional validation of the schemes

3.1 Riemann test-case

We compare the numerical schemes on one-dimensional Riemann test-cases: the initial data is a discontinuity between two constant states and the solutions are given by the superposition of waves separating constant states. In A, we give the form of these solutions with respect to the relative position of left and right states in the phase space. In the case of colliding states, explicit solutions can be numerically obtained. We thus consider the following Riemann test-case:

[TABLE]

on the domain $[0,1]$. The solution is made of two shock waves and an intermediate contact wave, see (33). The CFL condition (19) can be estimated by:

[TABLE]

For the current Riemann test-case with $\gamma=2$, the time step should satisfy $\Delta t\leqslant 0.4\Delta x$.

Comparison of the schemes ($\varepsilon=10^{-2}$)

In Figure 1, we represent the solution at time $t=0.1$ with the different schemes using $\Delta t=0.1\Delta x$. The $(\varrho,\boldsymbol{q},Z)$-method refers to the method introduced in Section 2.1 for the first order and in Section 2.2 for the second order scheme. The $(\varrho,q)$-method/SL refers to the method described in Section 2.3. For the latter scheme, we use the third order semi-Lagrangian scheme for the transport of the congestion density $\varrho^{\ast}$.

We observe that all the methods correctly capture the exact solution. The $(\varrho,\boldsymbol{q})$-method/SL better captures the contact discontinuity at $x\approx 0.487$ since we use a third order accurate scheme for the transport of $\varrho^{\ast}$. Limiters could be used to avoid overshoot and undershoot at this location.

Oscillations in momentum are brought forth at the discontinuity interface of the shock waves. These oscillations are larger for second order schemes due to dispersion effects. In Figure 2, we provide a zoom on these oscillations and compare the approximate solution to the exact one. The amplitudes of the oscillations are larger for the $(\varrho,\boldsymbol{q})$-method/SL method. This may be the counterpart of the decoupling of the variables $(\varrho,q)$ and $\varrho^{\ast}$: in the computation of the implicit pressure (see eq. (2.3), left-hand side), $\varrho$ and $\varrho^{\ast}$ are not taken at the same time.

Stiff pressure ($\varepsilon=10^{-4}$)

With this value of $\varepsilon$, the intermediate congested state has maximal wave speed equal to $\lambda_{\max}\approx 22$. Hence, taking time step $\Delta t$ equal to $0.1\Delta x$ does not ensure the resolution of the fast waves.

Figure 3 shows the solution at time $t=0.1$ using the $(\varrho,\boldsymbol{q},Z)$-method with second order in space accuracy. In the full second order scheme, the scheme switches automatically to a first order in time version of the scheme due to the large discontinuities in pressure, see Section 2.2. We observe that the waves are well captured. As previously, oscillations in momentum develop at schock discontinuities and we observe that the second order in time version of the scheme leads to large uppershoots. In Table 1, we report the $L_{1}$ error between numerical and exact solution: we point out that the numerical errors are of the same order of magnitude independantly of the value of $\varepsilon$. Quite similar results are obtained using the $(\varrho,\boldsymbol{q})$-method/SL.

3.2 Numerical convergence test-case

We here consider the following smooth initial data:

[TABLE]

on the domain $[0,1]$ and perdiodic boundary conditions. We compute a reference solution at time $t=0.05$ using the second order in space $(\varrho,q,Z)$-method with small space and time steps $\Delta x=5\times 10^{-5}$ and $\Delta t=0.1\,\Delta x$ (see Fig. 4).

Figure 5 shows the $L_{1}$ errors between approximate solutions and the reference solution at time $t=0.05$ when the space step $\Delta x$ goes to [math]. For first order scheme, time step is set to $\Delta t=5\times 10^{-6}$ while for second order schemes, time and space steps satisfy the relation $\Delta t=0.1\,\Delta x$ and both are varying.

We observe that all the schemes exhibit their expected convergence rates. We point out that $(\varrho,\boldsymbol{q},Z)$-method and $(\varrho,\boldsymbol{q})$-method/SL have the same level of numerical errors except for variable $\varrho^{\ast}$: $\varrho^{\ast}$ is better resolved with $(\varrho,\boldsymbol{q})$-method/SL. This is all the more the case when using the third order semi-Lagrangian scheme (on the right two plots of Fig. 5).

4 Two-dimensional numerical results

In this section we present the results of the numerical simulations in two-dimensions. As for domain we take the unit square with the mesh size $\Delta x=10^{-3}$ and the time-step $\Delta t=10^{-4}$. In the following we choose singular pressure (5) with the parameters $\varepsilon=10^{-4}$, $\alpha=2$, and the background pressure (3) with the exponent $\gamma=2$, if not stated differently.

First part is devoted to comparison of $(\varrho,\boldsymbol{q},Z)$-method and $(\varrho,\boldsymbol{q})$-method/SL described in Section 2. Second is an application of $(\varrho,q)$-method/SL to the evacuation scenario. Third order in space semi-Lagrangian scheme is applied.

4.1 Collision of 4 groups with variable congestion

In the unit square periodic domain we specify 4 squares, with the centers in points $(x_{c},y_{c})=\{(0.2,0.5),(0.5,0.2),(0.5,0.8),(0.8,0.5)\}$. The length of the side $l$ of each square equals $0.2$ (for every square we introduce the notation $\text{Square}((x_{c},y_{c}),l)$). We prescribe the initial momentum of $0.5$ pointing into the center of the domain provoking a collision. We consider three test cases varying in the initial congestion density, namely:

• Case 1:

$\varrho^{*}(x,0)=1.0$;

• Case 2:

$\varrho^{*}(x,0)=\begin{cases}0.80\text{ if }x\in\text{ Square}((0.2,0.5),0.2)\\ 1.20\text{ if }x\in\text{ Square}((0.5,0.2),0.2)\\ 0.80\text{ if }x\in\text{ Square}((0.8,0.5),0.2)\\ 1.20\text{ if }x\in\text{ Square}((0.5,0.8),0.2)\\ 1.00\text{ otherwise }\end{cases}$;

• Case 3:

$\varrho^{*}(x,0)=1+0.05(\cos(10\pi x)+\cos(24\pi x))(\cos(6\pi y)+\cos(34\pi y))$.

The results of our simulations for these three cases are presented in Figures 6, 7, and 8, subsequently and in the Movies c1.mp4, c2.mp4, and c3.mp4. We see that in case of constant congestion density (Case 1, Figure 6, Movie c1.mp4) the two schemes provide almost identical outcome. The essential difference appears when $\varrho^{*}_{0}$ varies. We see in Figure 7 (see also Movie c2.mp4) that the initial discontinuities of $\varrho^{*}$ are significantly smoothened by the $(\varrho,\boldsymbol{q},Z)$-method, while the $(\varrho,\boldsymbol{q})$-method/SL preserves the initial shape, which basically confirms our observations from Section 3.2. This is even more visible in Figure 8 (Movie c3.mp4), where the initial oscillations of $\varrho^{*}$ rapidly decay when simulated by the $(\varrho,\boldsymbol{q},Z)$-method.

Another interesting observation following from Figures 7, and 8 (Movies c2.mp4, c3.mp4) when compared to Figure 6 (Movie c1.mp4) is that the preference of the individuals $\varrho^{*}$ is significant factor to determine the density distribution even far away from the congestion zone.

Moreover, comparing Figure 7 (Movie c2.mp4) with Figure 6 (Movie c1.mp4), we see a clear influence of the density constraint on the velocity of the agents. Indeed, for the Case 2, there is a significant disproportion between the velocities in the $x$ and $y$ directions at time $t=0.150$ (see Figure 7 right). This corresponds to the fact that the agents moving toward the center along $y$ axis have ’more space’ to fill since $\varrho^{*}$ for those groups is higher than the one for the groups moving in the $x$ direction. This results in a certain delay between collisions in two directions.

4.2 Application to crowd dynamics

In this section we investigate an influence of the variable density $\varrho^{*}$ on a possible evacuation scenario. For this, we consider an impenetrable room in the shape of unit square, initially filled with uniformly distributed agents. There is an exit located at $x\in[0.4,0.6],\ y=0$ that allows for free outflow. The initial density $\varrho_{0}=0.6$ and the initial momentum is equal to ${\boldsymbol{0}}$. The desire of going to the exit is introduced in the system (1) (5) by adding the relaxation therm in the momentum equation

[TABLE]

where $\boldsymbol{w}$ is the desired velocity, and $\beta$ stands for the relaxation parameter. The desired velocity is given by a unit vector field, that points into the centre of the exit, $\boldsymbol{w}=\left(-x/((x\!-\!0.5)^{2}\!+\!y^{2}),-y/((x-0.5)^{2}\!+\!y^{2})\right)$ .

In the numerical scheme we apply splitting of the momentum equation between the transport and pressure part, and the relaxation (source) part, with the intermediate momentum $\boldsymbol{q}^{*}$. After the momentum is updated we perform implicit relaxation step, for given density $\varrho^{n+1}$,

[TABLE]

We use the $(\varrho,{\boldsymbol{q}})$-method/SL, which requires to solve the transport equation for $\varrho^{*}$. This is especially problematic in the corners of the domain, where the Dirichlet boundary condition is considered. This leads to oscillations of $\varrho$ and $\varrho^{*}$ close to these points. Nevertheless, we may observe, see Figure 9 and Figure 10 (see also Movies exit.mp4 and top.mp4), the so called stop-and-go behaviour, namely distinct high velocity regions in the domain, one in the vicinity of the exit and the second one that propagates in the direction opposite to flow.

This reflects an empirical observation that once a pedestrian arrives to the space of high congestion, he or she slows down or even stops until some space opens up in front. This kind of stop-and-go waves have been described, for example, by Helbing and Johansson in [21]. For the description of the real evacuation experiments we refer to [20], see also [19]. In the last of the mentioned papers the authors provide an experimental demonstration of the so called faster goes slower effect. This means that an increase in the density of pedestrians does not necessarily lead to a larger flow rate. Our simulations show that when the parameter $\varrho^{*}$ is low, the outflow of the individuals is slower. This is especially visible in the third row of Figures 9 and 10 presenting the evacuation scenario for the initial barrier density in the shape of the step function

[TABLE]

This observation can be also confirmed in terms of speed of evacuation.

Indeed, we performed analogous simulations for 3 cases of constant $\varrho^{\ast}_{0}$ equal to $0.9$, $1.0$. $1.1$ show that the speed of emptying the room is bigger the bigger value of $\varrho^{*}_{0}$. To see this we have measured the mass remaining in the room at time $t=1$ and it is equal to $0.51030$, $0.048037$, and $0.457123$, respectively. We have moreover observed that evacuation speed of the room with individuals of the average congestion preference equal to $1$ initially can be improved by placing the individuals with higher $\varrho^{*}_{0}$ closer to the exit. This is illustrated in the Figures 9 and 10 the second row, for which, the initial congestion preference $\varrho^{\ast}_{0}$ equals

[TABLE]

The random distribution of preferences of the individuals with expected value equal to $1$, on the other hand, corresponds to the increase of the evacuation time (see Figures 9 and 10 the bottom row).

5 Conclusion

In this paper, we are interested in the numerical simulation of the Euler system with a singular pressure modeling variable congestion. As the stiffness of the pressure increases ($\varepsilon$ tends to [math]), the model tends to a free boundary transition between compressible (non-congested) and incompressible (congested) dynamics.

To numerically simulate the asymptotic dynamics, we propose an asymptotic preserving scheme based on a conservative formulation of the system and the methodology presented in [14]. We also propose a second order accuracy extension of the scheme following [12]. We then study the one-dimensional solutions to Riemann test-cases, their asymptotic limits and validate the code. We compare the results with those obtained with the scheme proposed in [15]. This latter scheme enables to better approximate the congestion density (at the contact wave) as soon as we use high accuracy in the advection of the congestion density. On the other hand, the former scheme seems to better preserve maximum principle on that variable. On two-dimensional simulations, we finally show the influence of this variable congestion density on the dynamics and show that the model exhibit stop-and-go behavior.

The two schemes generate oscillations in momentum variable at discontinuities between congested and non-congested domain. This feature was already mentioned in [14]. This is all the more the case for the second order accuracy schemes. Specific method should be designed to cure this artefact.

Appendix A Solution to the Riemann problem

The one-dimensional Riemann problem for the system (7) is the following initial-value problem:

[TABLE]

where $p_{\varepsilon}(Z)=\pi_{\varepsilon}(Z)+p(Z)$, and

[TABLE]

The purpose of this section is to find possible weak solution to (29) (30). We will also consider the limit of these solutions as $\varepsilon\rightarrow 0$.

As already mentioned in the introduction, the system (29) is strictly hyperbolic provided $p_{\varepsilon}^{\prime}(Z)>0$, see (8). The associated characteristic fields are given by:

[TABLE]

where $v=q/\varrho$ is the velocity. The second characteristic field is linearly degenerate (since $\nabla\lambda_{2}\cdot r_{2}=0$). The two others characteristic field are genuinely non-linear.

We now present the elementary wave solutions of the Riemann problem.

A.1 Elementary waves

Shock discontinuities

A shock wave is a discontinuity between two constant states, $(\varrho,q,Z)$ and $(\hat{\varrho},\hat{q},\hat{Z})$, travelling at a constant speed $\sigma$. We now fix the left (or right) state $(\hat{\varrho},\hat{q},\hat{Z})$ and look for all triples $(\varrho,q,Z)$ that can be connected to $(\hat{\varrho},\hat{q},\hat{Z})$ by the shock discontinuity. Across the shock, the Rankine-Hugoniot conditions must be satisfied meaning that:

[TABLE]

where $[a]:=a-\hat{a}$ denotes the jump of quantity $a$. Treating $\varrho$ as a parameter, we check that the two admissible states are of the form $(\varrho,q_{h,\pm}(\varrho),Z(\varrho))$ with $q_{h,\pm}=\varrho v_{h,\pm}(\varrho)$ and

[TABLE]

The shock speed therefore equals:

[TABLE]

These solutions can also be expressed as functions of $Z$:

[TABLE]

Note that the maximal density ($\varrho^{\ast}=\varrho/Z$) does not jump across a shock discontinuity. Expanding $(\varrho(Z),q_{h,\pm}(Z),Z)$ around $Z=\hat{Z}$, we obtain

[TABLE]

Note that $(\varrho(Z),q_{h,-}(Z),Z)$ is tangent at $(\hat{\varrho},\hat{q},\hat{Z})$ to $r_{1}(\hat{\varrho},\hat{q},\hat{Z})$, therefore $v_{h,-}$ corresponds to the $1$-characteristic field, analogously $v_{h,+}$ corresponds to the $3$-characteristic field. The graph of $Z\mapsto v_{h,-}(Z)$ (resp. $Z\mapsto v_{h,+}(Z)$) is called the 1-Hugoniot curve (resp. 3-Hugoniot curve) issued from $(\hat{v},\hat{Z})$.

To check the admissibility of the discontinuity, we need to check the entropy condition. If $(\hat{v},\hat{Z})$ is the left state, the right states that can be connected to it by an entropic shock wave are those located on the 1-shock curve $\left\{\left(v_{h,-}(Z),Z\right):Z>\hat{Z}\right\}$ or the 3-shock curve $\left\{(v_{h,+}(Z),Z):Z<\hat{Z}\right\}$. Indeed, on these curves the associated eigenvalue is decreasing. If on the other hand, $(\hat{v},\hat{Z})$ is the right state, the left states that can be connected to it by an entropic shock wave are those located on the 1-shock curve $\left\{(v_{h,-}(Z),Z):Z<\hat{Z}\right\}$ or the 3-shock curve $\left\{(v_{h,+}(Z),Z):Z>\hat{Z}\right\}$. Indeed, on these curves the associated eigenvalue is increasing.

Rarefaction waves

The rarefaction waves are continuous self-similar solutions, $(\varrho(t,x),q(t,x),Z(t,x))=(\varrho(x/t),q(x/t),Z(x/t))$, connecting two constant states $(\varrho,q,Z)$ and $(\hat{\varrho},\hat{q},\hat{Z})$. They thus satisfy the following differential equations:

[TABLE]

Denoting $q(s)=\varrho(s)\tilde{v}_{i,\pm}(s)$ and parametrizing by $\varrho$, we obtain:

[TABLE]

From the first and third equation of (31), we have $(\varrho/Z(\varrho))^{\prime}=0$, and so, $\varrho/Z(\varrho)=\hat{\varrho}/Z(\hat{\varrho})$. This means that as in the case of shock discontinuities the maximal density $\varrho^{*}$ does not jump. Denoting $\varrho^{*}=\varrho/Z(\varrho)$ and making the change of coordinates $v_{i,\pm}(Z)=\tilde{v}_{i,\pm}(\varrho)$ with $\varrho=\varrho^{\ast}Z$, we thus have:

[TABLE]

Hence, the states satisfy:

[TABLE]

where $F_{\varepsilon}$ is an antiderivative of $Z\mapsto\frac{1}{Z}\sqrt{\frac{1}{\varrho^{*}}p_{\varepsilon}^{\prime}(Z)}$.

The graph of $Z\mapsto v_{i,+}(Z)$ (resp. $Z\mapsto v_{i,-}(Z)$) is called the 1-integral curve (resp. 3-integral curve) issued from $(\hat{v},\hat{Z})$. If $(\hat{v},\hat{Z})$ is a left state, the right states that can be connected to it by an entropic rarefaction wave are those located on the 1-integral curve $\left\{(v_{i,-}(Z),Z):Z<\hat{Z}\right\}$ or the 3-integral curve $\left\{(v_{i,-}(Z),Z):Z>\hat{Z}\right\}$. Indeed, on these curves the associated eigenvalue is increasing. If $(\hat{v},\hat{Z})$ is a right state, the left states that can be connected to it by an entropic rarefaction wave are those located on the 1-integral curve $\left\{(v_{i,-}(Z),Z):Z>\hat{Z}\right\}$ or the 3-integral curve $\left\{(v_{i,-}(Z),Z):Z<\hat{Z}\right\}$. Indeed, on these curves the associated eigenvalue is decreasing.

Contact discontinuities

Since the second characteristic field is linearly degenerate, there are linear discontinuities that propagate at velocity $\lambda_{2}=\hat{v}$. Let us write the Rankine-Hugoniot conditions:

[TABLE]

From the first relation, we obtain $v=\hat{v}$ and then the second relation states that the pressure jump is zero. By strict monotony of the pressure, it implies that $Z=\hat{Z}$ and the third equation is satisfied. Along this discontinuity, the velocity and the pressure are thus conserved. Note that every density jump is possible.

A.2 Solution to Riemann problem

Let $(\varrho_{\ell},q_{\ell},Z_{\ell})$ and $(\varrho_{r},q_{r},Z_{r})$ be the left and right initial states (30). The solutions to Riemann problems are determined as follows. First, in the $(v,Z)$ plane, find out the intersection state $(v_{m},Z_{m})$ of the 1-st integral/Hugoniot curves issued from $(v_{\ell},Z_{\ell})$ and the 3-rd integral/Hugoniot curves issued from $(v_{r},Z_{r})$. Then, compute the two densities $\varrho_{m,\ell}$ and $\varrho_{m,r}$ so that the congestion density across the two non-linear waves is conserved. Then we connect the two distinct intermediate states by a contact discontinuity. We finally end up with the following solution:

[TABLE]

where $\varrho_{m,\ell}=Z_{m}\varrho_{\ell}/Z_{\ell}$ and $\varrho_{m,r}=Z_{m}\varrho_{r}/Z_{r}$. The nature of the non-linear waves (rarefaction or shock) depends on the relative position of the states $(v_{\ell},Z_{\ell})$, $(v_{r},Z_{r})$ in the $(v,Z)$ plane.

A.3 Limit $\varepsilon\rightarrow 0$

We are now interested in the asymptotic behaviour, when $\varepsilon\to 0$ of the Hugoniot $v_{h,\pm}^{\varepsilon}$ and the integral curves $v_{i,\pm}^{\varepsilon}$ obtained in the previous paragraph for the elementary waves. We have the following result.

Proposition 1**.**

The graph of the Hugoniot curve, $\left\{(Z,v_{h,\pm}^{\varepsilon}(Z)):Z\in[0,1)\right\}$, tends to the union of the set $\left\{(Z,v_{h,\pm}^{0}(Z)):Z\in[0,1)\right\}$ and the horizontal half straight line $\left\{(1,v):v\in[v_{h,\pm}^{0}(1),+\infty)\right\}$.

The graph of the integral curve, $\left\{(Z,v_{i,\pm}^{\varepsilon}(Z)):Z\in[0,1)\right\}$, tends to the union of the set $\left\{(Z,v_{i,\pm}^{0}(Z)):Z\in[0,1)\right\}$ and the horizontal half straight line $\left\{(1,v):v\in[v_{i,\pm}^{0}(1),+\infty)\right\}$.

The proof of this proposition uses the convexity of the pressure and are similar to the one developed in [16].

Regarding the Riemann problem in the limit $\varepsilon\rightarrow 0$, the intersection point of the 1-st integral/Hugoniot curves issued from $(v_{\ell},Z_{\ell})$ and the 3-rd integral/Hugoniot curves issued from $(v_{r},Z_{r})$, denoted by $(v_{m}^{\varepsilon},Z_{m}^{\varepsilon})$, has either a limit $(v_{m}^{0},Z_{m}^{0})$ with $0\leqslant Z_{m}^{0}<1$ or tends to a congested state $(\bar{v},1)$. Then finding a solution can be divided into the following steps:

• (1)

compute the intersection $(v_{m}^{0},Z_{m}^{0})$ of the 1-st integral/Hugoniot curves and 3-rd integral/Hugoniot curves;

• (2a)

if $Z_{m}^{0}<1$, the solution is as described in the previous section, it is a usual Riemann solution of the hyperbolic system with no congestion pressure;

• (2b)

if $Z_{m}^{0}\geq 1$, then the congested state is given by the following proposition.

Proposition 2** (Case $Z_{m}^{0}\geq 1$.).**

The solution consists in three waves:

[TABLE]

where the intermediate velocity $\bar{v}$ and pressure $\bar{p}$ satisfy:

[TABLE]

the intermediate densities are given by:

[TABLE]

and the shock speeds $\sigma_{-}$, $\sigma_{+}$ are given by:

[TABLE]

This proposition can be proven using similar arguments as in [16].

Below, on Figure 11 we present two different solutions to the Riemann problem (29)-(30). Depending on the initial location of the left and right states, the intersection state $(v_{m},Z_{m})$ might be a congested state or not.

Appendix B Fully discrete scheme in dimension 2

We consider the computational domain $[0,1]\times[0,1]$ and spatial space steps $\Delta x=1/N_{x},\Delta y=1/N_{y}>0$, with $N_{x},N_{y}\in\mathbb{N}$: the mesh points are thus $\boldsymbol{x}_{i,j}=(i\Delta x,j\Delta y)$, $\forall(i,j)\in\left\{0,\ldots,N_{x}\right\}\times\left\{0,\ldots,N_{y}\right\}$. Let $\varrho^{n}_{i,j}$, $\boldsymbol{q}^{n}_{i,j}$, $Z^{n}_{i,j}$, $\varrho^{\ast\,n}_{i,j}$ denote the approximate solution at time $t^{n}$ on mesh cell $[i\Delta x,(i+1)\Delta x]\times[j\Delta x,(j+1)\Delta x]$.

The two-dimensional version of (11) reads:

[TABLE]

where fluxes $F^{n+1}$, $\boldsymbol{G}^{n}$, $H^{n+1}$ (in the first spatial direction) are defined:

[TABLE]

with

[TABLE]

Fluxes $\widetilde{F}^{n+1}$, $\widetilde{\boldsymbol{G}}^{n}$, $\widetilde{H}^{n+1}$ in the second spatial direction are defined by:

[TABLE]

with

[TABLE]

The upwindings $D_{\varrho}$, $\boldsymbol{D}_{\boldsymbol{q}}$, $D_{Z}$ are defined similarly as for the one-dimensional case (sse (15)-(16)).

The implicit pressure in (35) is discretized by the centered difference:

[TABLE]

Inserting equation (35) into (36), we obtain:

[TABLE]

where terms $\bar{H}^{n}$ and $\bar{\widetilde{H}}^{n}$ have the same expressions as (39)-(42) but where all quantities are taken explicitly.

Appendix C Second order in time $(\varrho,q)$-method/SL

The second order accuracy scheme for the $(\varrho,q)$-method/SL is based on a Strang splitting between advection of congestion density and advection of $(\varrho,\boldsymbol{q})$. It consists in the following steps:

Compute $\varrho^{\ast\ n+1/2}$ by solving the advection of over $\Delta t/2$

[TABLE] 2. 2.

Compute $(\varrho^{n+1},q^{n+1})$ with the RK2CN scheme as proposed in [12]:

First step (half time step):

[TABLE]

Second step (full time step):

[TABLE]

where $\mathcal{D}_{\varrho}$, $\mathcal{D}_{q}$ denote the numerical diffusion coming from fluxes. 3. 3.

Advection of $\varrho^{\ast}$ on $\Delta t/2$ time step

[TABLE]

A second order in time version of the semi-Lagrangian scheme has to be used. We here consider the second order Taylor approximation of the caracteristic line whose one-dimensional version reads:

[TABLE]

where $v_{i}=q_{i}/\varrho_{i}$ for all $i$ and $a_{i}$ is an upwind finite difference approximation of the first derivative of the velocity: $a_{i}=(v_{i}-v_{i-1})/\Delta x$ if $v_{i}>0$ and $a_{i}=(v_{i+1}-v_{i})/\Delta x$ if $v_{i}\leqslant 0$.

Acknowledgements

P.D. acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/M006883/1, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). P.D. is on leave from CNRS, Institut de Mathématiques de Toulouse, France. P.M. acknowledges the support of MAThematics Center Heidelberg (MATCH). E.Z. was supported by the the Department of Mathematics, Imperial College, through a Chapman Fellowship, and by the Polish Ministry of Science and Higher Education grant Iuventus Plus no. 0888/IP3/2016/74.

Data availability

No new data were collected in the course of this research.