Transport of congestion in two-phase compressible/incompressible flows

We study the existence of weak solutions to the two-phase model of crowd motion. The model encompasses the flow in the uncongested regime (compressible) and the congested one (incompressible) with the free boundary separating the two phases. The congested regime appears when the density in the uncongested regime $\varrho^*(t,x)$ achieves a threshold value $\varrho^*(t,x)$ that describes the comfort zone of individuals. This quantity is prescribed initially and transported along with the flow. We prove that this system can be approximated by the fully compressible Navier-Stokes system with a singular pressure, supplemented with transport equation for the congestion density. We also present the application of this approximation for the purposes of numerical simulations in the one-dimensional domain.


Introduction
Mathematical modelling of crowd dynamics is a very challenging problem. The recently rapidly growing literature on this subject covers several parallel approaches. We can distinguish, for example, the mean-field game models [14,25], in which the individuals behave as the players following some strategy, or optimizing certain cost; the microscopic models which describe precise position and velocity of an individual (Individual-Based-Models) using Newtonian framework [20][21][22]35]; or the macroscopic models formulated in the language developed for the fluids [8,9,23,24,34]. The behaviour of the crowd in the later is characterised by some averaged quantities such as the number density or mean velocity. The macroscopic models, although less precise than the microscopic ones, are computationally more affordable. Moreover, they allow for asymptotic studies that proved to be useful for understanding various aspects like: swarming or pattern formation observed in the experiments.
Our aim is to analyze the free-boundary two-phase fluid system that could be used to model the congestions in the large group of individuals in a bounded area. Our individuals do not follow any neighbour trying to align their velocity or reach a certain evacuation point. They are just the agents that have their individual preferences for how close they let the closest neighbour to approach and they carry this information with them in the course of motion. We prescribe their initial velocity that determines their desired direction of motion and check how the individual preferences as well as the initial distribution of the agents determines creation of congestions.
Our system writes as follows: with the unknowns: = (t, x) -the mass density, u = u(t, x) -the velocity vector field, * = * (t, x) -the congestion density, also referred to as the barrier or the threshold density, and π -the congestion pressure appearing only when = * .
The barotropic pressure is an explicit function of * p * = * γ , γ > 1, (2) and plays the role of the background pressure. The stress tensor S is a known function of u, characteristic for the Newtonian fluid, namely S = S(u) = 2µ D(u) + λdivu I, µ > 0, 2µ + λ > 0, where D(u) = (∇u + ∇ T u)/2 denotes the symmetric part of the gradient of u, and I = I 3 is the identity matrix. It is justified to call the above system the two-phase system because for (t, x) < * (t, x) it behaves as the compressible Navier-Stokeas system with the barotropic pressure, while when the congestion is achieved, i.e. for (t, x) = * (t, x), the system behaves like the incompressible Navier-Stokes equations. We thus observe a switching between two phases: compressible and incompressible depending on the size of the density ratio * . The fluid systems with congestion constraints have been recently intensively studied, especially in the hyperbolic regime [1,3,6]. The first analytical result for system (1) with * = 1 is due to P.-L. Lions and N. Masmoudi [26], who showed that it can be obtained as a limit of compressible Navier-Stokes equations with barotropic pressure γ with γ → ∞, similar studies were performed recently for the model of tumour growth [33].
In the system (1) * models preferences of the individuals, it is given initially and then transported with the flow. Therefore, * depends on time and position, but more importantly it depends on initial configuration * 0 . The form of * relaxes the restrictions from the models studied in [7,32], where the threshold density * was either assumed to be constant or independent of time. This allows to cover more physical applications. Including the transport of the congestion density * allows also to study the system (5) with the contribution from the pressure in the form of the pure gradient, without factor * as it was done in [32].
We will consider the system (1) in the 3-dimensional domain Ω with the smooth boundary ∂Ω, and the Dirichlet boundary conditions for the velocity vector field The initial conditions are given by: and we assume that they satisfy: in Ω, m 0 ∈ L 2 (Ω), Moreover, we assume that in the region of the absence of the individuals 0 (x) = 0, the congestion density is equal to a constant value, being the characteristic mean preference of the group: * The main result of this paper is the existence of solutions to the system (1) under the aforementioned assumptions on the constitutive relations and the initial condition, in the sense of the following definition.
The main theorem of the paper states as follows.
Theorem 1 Let the initial conditions 0 , m 0 , * 0 satisfy the conditions above. Then the system (1) with p and S given by (2), (3) respectively, has a weak solution in the sense of Definition 1.

Remark 1
The same result holds for the lower dimensions, d=1,2.
The core of the proof of Theorem 1 is to show that the system (1) can be obtained as a limit when ε → 0 of the following approximation where the π ε stands for the singular pressure of the form A similar form of the pressure [4] or [11]. Singularities of the pressure of this type were also previously studied in the context of traffic models [2,4,5], collective dynamics [10,11], or granular flow [28,31]. Our first goal is to prove the existence of solutions to a certain reformulation of the system (5) for ε > 0 fixed. Similar system, but with the barotropic form of the pressure π = γ was studied recently in [27], see also the stability result from [29] and the low Mach number analysis from [17]. It was proved that for certain values of parameter γ, there is an equivalence between three different definitions of weak solutions to (5) based on the reformulations of the transport equation for the entropic variable.
Here we will essentially use the formulation involving a new unknown Z, that for , * smooth enough can be identified with the density fraction Z = * . We can formally check, that dividing (5a) by * , multiplying (5c) by − * 2 , and summing the resulting expressions, the system (5) can be transformed to the following one with the initial data In consistency with the assumptions on 0 , m 0 , * 0 from the previous section, we postulate that 0 , m 0 , Z 0 satisfy The first condition from (9) should be understood as the restriction of the initial congestion density, having in mind our notation Z = * we see that the condition above means that * 0 cannot be zero on the regions the density 0 in positive. The left condition means that the congestion density is initially bounded. The restriction on the integral Ω Z dx < |Ω| means that the assumption 0 = * 0 holds on a set of non-zero Lebesgue measure.
Definition 2 (Weak solution of the approximate system) Suppose that the initial conditions satisfy (9). We say that the triplet ( , u, Z) is a weak solution of problem (7) with the initial and boundary conditions (8), (4) if , and for any T > 0 we have: for all test functions ϕ ∈ C 1 ([0, T ] × Ω); for all test functions ϕ ∈ C 1 ([0, T ] × Ω); (iv) the energy inequality holds for a.a T > 0, where Our second goal is to show the convergence of the weak solutions to the system (7), to the solutions of the limit system: In the limit ε, the uncongested (compressible) flow changes to the incompressible when the density hits the value * . When this is the case the dynamics of the system is modified abruptly, meaning that the transition from the uncongested motion ( < * ) to the congested motion ( = * ) is very sudden. The weak solutions to the limit system are defined as follows: Definition 3 (Weak solution of the limit system with Z) A quadruple ( , u, Z, π) is called a weak solution to (16) with (4) and (8) if equations (16a), (16b), (16c) are satisfied in the sense of distributions, the constraints (16d) satisfied a.e. in (0, T ) × Ω, the divergence free condition (16e) is satisfied a.e. in {Z = 1}, and the following regularity properties hold Moreover, π is sufficiently regular so that the condition is satisfied in the sense of distributions on (0, T ) × Ω.
The convergence result reads as follows.
be a a sequence of weak solutions to the approximate system (7), (6). Then, for ε → 0, the sequence ( ε , u ε , Z ε ) converges to the weak solution of (16) in the sense of Definition 3 More precisely, for any p < ∞, and u ε → u weakly in L 2 (0, T ; W 1,2 (Ω, R 3 )). Moreover, The paper is organized in the following manner. In Section 2, we present details of approximation and prove the existence of solutions to the system (7) for ε fixed. Then, in Section 3, we recover the two-phase system (16) by letting ε → 0. After this, in Section 4 we recover the solution to the original two-phase system (1). Finally, in Section 5 we briefly describe the numerical scheme and present computational examples that illustrate the behaviour of approximate solutions to the system (1).

The existence of solution for ε fixed
When ε is fixed, say ε = 1, the system (7) shares the features of the system considered in [32] and in the recent paper [27]; in our proof of existence of solutions we will recall some elements of these two approaches in order to avoid repetitions.

Formulation of the approximate problem
The first level of approximation introduces truncation parameter δ in the singular pressure and artificial pressure κ K , with K sufficiently large to be determined in the course of the proof. We consider with κ, δ > 0, and π δ , p κ given by We drop the subindex δ when no confusion can arise, and we introduce the notion of weak solution for the system (17).
Definition 4 (Weak solution of the approximate system ) Suppose that the initial conditions satisfy (9). We say that the triplet ( , u, Z) is a weak solution to the problem (17) with the initial and boundary conditions (8) and (4) if , and for any T > 0 we have: for all test functions ϕ ∈ C 1 ([0, T ] × Ω); for all test functions ϕ ∈ C 1 ([0, T ] × Ω); (iv) the energy inequality holds for a.a. T > 0, where We have the following existence result for solutions defined by Definition 4 (see also [27], Theorem 2).
Then there exists a weak solution ( , u, Z) to problem (17), (18) with boundary conditions (4), in the sense of Definition 4.
Moreover, (Z, u) solves (17c) in the renormalized sense, i.e. (Z, u), extended by zero outside of Ω, satisfies in the sense of distributions on (0, In addition, The proof of Theorem 3 requires further modification of the system (17) with two additional approximation levels involving the parabolic regularisation of two continuity equations and the Galerkin approximation of the velocity. This approximation allows, in particular, to deduce inequalities (27). At this point existence of regular solution can be obtained following the arguments presented in [27]. Also the compactness arguments needed to recover system (17) are analogous, therefore we skip this part and focus only on the a-priori estimates needed to perform the limit passages δ → 0, κ → 0.
Having the existence of solutions to the system (17)-(18), we show that this solution can be used to recover the weak solution to the system (7), where only the parameter ε is present.
More precisely, for any p < ∞, and Moreover, In the next two sections we prove Theorem 4. Starting from the weak solution to the system (17), we first derive uniform bounds and then we let δ → 0. The modification of the reasoning needed to perform the limit passage κ → 0 is explained at the end.

Uniform estimates
In this section we obtain the uniform estimates for the weak solutions to the system (17)- (18). This involves three a-priori estimates: the standard energy estimate, and two estimates involving the application of the Bogovskii operator, one of which gives us the uniform bound for the singular pressure and the other the uniform integrability of the pressure. All of the aforementioned estimates are essential for performing the limit passage δ, κ → 0, as well as the last limit passage ε → 0. However, as we shall see later on, the uniform integrability of the pressure will no longer be valid in this case.

The energy estimate
We present a formal computation that can be made rigorous at the level of the Galerkin approximation of the velocity. Multiplying the momentum equation (17b) by u and integrating by parts with respect to space, yield the energy equality Note that the second term in our case is different than in [32], there is no additional * in front of the gradient. This however allows us to proceed straightforwardly, for reeder's convenience we repeat here the derivation where we used the equation (17c), we denoted Q κ,δ (Z) = π δ (Z)+p κ (Z) Z , and Γ κ,δ is a solution of the following ODE: Using the definition of Q κ,δ and π κ,δ we can express Γ κ,δ as in (24). Integrating (28) with respect to time, and using the definition of the stress tensor (3) we get the following uniform estimates

The integrability of the pressure
The energy estimate obtained above is insufficient to control the L 1 norm of the pressure, because the singularity appearing in Γ κ,δ (Z) at Z = 1 is of lower order than the singularity for π κ,δ (Z) (take f.i. α = 2, β = 1 in (6) and use (24)). Therefore, further estimates are needed. The first of them is obtained by testing the momentum equation by the function where φ is smooth and compactly supported in the interval (0, T ), and B is the Bogovskii operator, whose main properties can be found, for example, in [30,Lemma 3.17], and [32,Appendix]. In particular, B is a linear operator such that for any f ∈ L p (Ω), 1 < p < ∞, and Moreover Remark 2 Note that, alike for , the integral of Z over the space is an invariant of the motion, therefore we have according to (9), at least for sufficiently regular solutions. In what follows we denote Ω Z 0 (x) dx = M Z .
Remark 3 Function ϕ must be of a certain regularity in order to use it as a test function in (20). In fact, it follows from the construction of the solution in [27], that taking K sufficiently large guarantees admissibility of this function.
Using in (20) the test function (30) results in the following equality whose r.h.s. can be bounded using the uniform estimates (29) together with (27), provided that K is sufficiently large, say K > 4. Therefore, one gets the following estimate, which is now uniform with respect to κ and δ We then consider two complementary subsets of (0, T ) × Ω: The l.h.s. of (32) can be easily controlled on Σ 1 , because on this subset Z δ stays far away from the singularity, for Σ 2 we have and so, from (32) it follows that as well as The estimate (32) implies that the L K+1 norm of Z δ is bounded uniformly in δ, but not uniformly in κ, which suggests that the passage to the limit δ → 0 should be performed as first.

The equi-integrability of the singular pressure
Because π δ is a nonlinear function of Z, identification of the limit in this term, after letting δ → 0 cannot be justified only by the uniform L 1 bound. Following the idea from [19], we can prove that the pressure π δ enjoys some additional estimate near to the singularity Z = 1. Indeed, for K > 6, β > 5/2 and uniformly with respect to δ one has The proof follows by testing the momentum equation (20) by the function of the form where φ is smooth and compactly supported in the interval (0, T ). For the details of this estimate we refer to [32], for β > 3 and to [18] for β > 5 2 . One of the difficulties in the proof of analogue of (36) presented in [32] concerned the renormalization of the equation for * . Here this problem does not appear anymore, since (Z δ , u δ ) is by definition a distributional solution of the renormalized transport equation (25).
2.3 Passage to the limit δ → 0

Convergences following from the uniform estimates
Using the uniform estimates (27), (29), and the Hölder inequality, we can deduce that up to a subsequence Using the continuity equation (17a) and (17c), the the two first convergences can be strengthened to which together with the weak convergence of the velocity gradient, after using the momentum equation (17b) implies that This in turn, assuming that K is sufficiently large so that the imbedding of L 2K for some q > 1. Moreover, the uniform bound (33) together with the growth condition β > 5/2 in (18) imply that and thus also on the account of (27). Finally, from the uniform bounds (33) we can extract the subsequences such that for some π(Z), Zπ(Z) that need to be determined. Note, however, that (36) together with De La Vallée-Poussin criterion allow us to deduce the equi-integrability of the pressure, therefore the first limit can be strengthened to π δ (Z δ ) → π(Z) weakly in L 1 ((0, T ) × Ω).
As for the second convergence, we cannot say that much immediately. However, using the fact that π δ 1 (·) ≥ π δ 2 (·) provided δ 1 ≤ δ 2 , we can estimate for any smooth, nonnegative, compactly supported function φ(x, t) where the last inequality is a consequence of the fact that · → π δ * (·) is non-decreasing function for fixed δ * . Using again equi-integrability of the pressure and letting δ * → 0 we verify that the r.h.s. of (45) vanishes and thus in the sense of distributions. In order to say something more, we need to investigate the strong convergence of the sequence Z δ , which is the purpose of the next section.

Strong convergence of Z δ
It was observed in [27], that the strong convergence of Z δ does not imply the strong convergence of δ . The proof of the strong convergence of Z δ requires an analogue of the effective flux equality for the barotropic compressible Navier-Stokes equations. In our case it can be written as follows.
Lemma 5 Let δ , u δ , Z δ be the sequence of approximate solutions enjoying the properties from above. Then, at least for a subsequence for any ψ ∈ C ∞ c ((0, T )) and φ ∈ C ∞ c (Ω). The proof of this fact can be seen as a special case of an analogous result proven in [27, cf. Lemma 11]. On the account of (46) and the monotonicity of p κ (·) we obtain from (47) for any ψφ ≥ 0. Recall that Z δ satisfies the renormalized continuity equation (25) with b specified in (26). By density argument and standard approximation technique, we may extend the validity of (25) The renormalization technique applied to the barotropic Navier-Stokes system is due to DiPerna and Lions [13], and the above extension can be found, for example, in [16].
We can now write the renormalized continuity equation with b(Z δ ) = Z δ log Z δ : Passing to the limit δ → 0 + , we hence obtain Writing an analogous equation for the limit function Z and subtracting it from the above, we obtain satisfied in the sense of distributions on (0, T )×R 3 . Integrating the above equation with respect to time and space, and using the convexity of the function s → s log s we get from (49) that which is an opposite to (48). Therefore, recalling (49) we see that Z log Z(t, x) = Z log Z(t, x) almost everywhere in (t, x) ∈ (0, T )×Ω, which yields the strong convergence of Z δ in L p ((0, T )× Ω) for any p < K + 1. With this at hand, we verify that (44) can be replaced by In order to complete the proof of Theorem 4, we have to show that we are allowed to let κ → 0. Note that all the uniform estimates obtained above stay in force independently of κ. Indeed, the only issue here is to see that the uniform L K+1 bound on Z κ , that was previously obtained from (33) follows now directly from (42). On the account of (27), the same is true also for the sequence κ . The proof of Theorem 4 is now complete. 2 3 Passage to the limit ε → 0 The purpose of this section is to prove our main Theorem 1. For technical reasons we first perform the limit ε → 0 in the auxiliary system (7) proving Theorem 2 and then we prove equivalence between systems (16) and (1) in the certain class of solutions.

Convergence following from the uniform estimates
The estimates performed in the previous section give rise to several estimates that are uniform with respect to ε. Indeed, passing to the limit in the energy estimate we obtain (50) Passing to the limit in (42) and in (27) we obtain in particular both sequences Z ε , ε are uniformly bounded in L p ((0, T ) × Ω) for p ≤ ∞. Therefore, by means of the arguments from the previous section we get, up to the subsequence, that This information allows us to pass to the limit in all terms of the system (7), apart from the nonlinear pressure terms p(Z) and π ε (Z). Repeating the Bogovskii type estimate with the test function (30) we obtain however, the estimate (36) does not hold anymore. Therefore the convergence in the sense of measures is the most we can hope for, we have For the background pressure, due to (51), we have for any p < ∞. At this point, we can identify the second limit in (54) using the explicit form of the pressure (6). We have thus letting ε → 0 and observing that the last term converges to zero strongly, we obtain the relation in the sense of the measures. The recovery of the constraint condition (1 − Z)π = 0 and the identification of the limit p(Z) = p(Z) require stronger information about convergence of Z ε .

Strong convergence of Z ε
The nowadays well known technique of proving the strong convergence of the density in the compressible barotropic Navier-Stokes equations involves the study of propagation of the so called oscillation defect measure [15]. Such level of precision is not needed in our case, because, due to the singularity of the pressure argument Z, we have sufficiently high integrability of the pressure in order to apply the DiPerna-Lions technique [13]. Nevertheless, a variant of effective viscous flux equality is still needed. It can be derived the same way as in Lemma 5, after observing that the inverse divergence operator ∇∆ −1 [1 Ω Z] is regular enough to be used as a test function in the limiting momentum equation.
With this information the statement of Lemma 5 adapted to the ε-labelled sequences gives rise to the equality for any ψ ∈ C ∞ c ((0, T )) and φ ∈ C ∞ c (Ω). Let us now explain the meaning of the product Zπ on the r.h.s. of (57). To this end, one needs to come back to the limiting momentum equation and use the bounds (50), (51) to justify that π is in fact more regular than it follows just from (54). Indeed, we have Moreover from the equation for Z, we easily get Regularizing in space and time the limits Z and π using the standard multipliers ω n , Z n = Z * ω n , π n = π * ω n , we can clarify the meaning of Zπ by writing Zπ = Z n π n + (Z − Z n )π n + Z(π − π n ), and by passing to the limit with the support of mollifying kernel, see [32] for more details.
From (57) it follows that where to get to the r.h.s. of the above we have used subsequently: monotonicity of p(·), (56), and the limit of (51). Since both pairs (Z ε , u ε ) and (Z, u) satisfy the renormalized continuity equation, we can use the renormalization in the form b(z) = z log z to justify that Note however, that this property is not transferred to the sequence ε , for which we only have (52). Nevertheless, using this information and formula (60) we can justify that which together with (56) implies (16f). It remains to show the condition (16e), or rather its compatibility with the other conditions in system (16). This follows from the following lemma proven by Lions and Masmoudi in [26], that we recall here without the proof.

Recovery of the original system
In Section 3 we proved that system (16) possesses a weak solution in the sense of Definition 3. Our next aim is to prove that this solution can be identified with the solution to the original problem (1). In other words, we need to deduce the existence of * ∈ L ∞ ((0, T ) × Ω) satisfying the transport equation, such that the measure in the momentum equation vanishes for = * . We proceed similarly to [27]. Note that When extended by 0 outside Ω the couples ( , u) and (Z, u) satisfy the renormalized continuity equations, due to uniform L ∞ bounds for both ε and Z ε . Therefore, we may test the equations (16a), (16c) by ω n (x − ·), where ω n is a standard mollifier, which leads to ∂ t n + div( n u) = r 1 n , satisfied a.e. in (0, T ) × R 3 , where by a n we denoted a * ω n . It follows from the Friedrichs commutator lemma, see e.g. [16,Lemma 10.12], that r 1 n , and r 2 n converge to 0 strongly in L 1 ((0, T ) × R 3 ) as n → ∞.

Numerical scheme
Numerical simulation of two-phase flows with free boundary requires to design a method that captures phase transition and the limit behaviour. In our case, the main difficulty is to propose a scheme that is independent of singular pressure parameter ε (6). This property is referred to as the Asymptotic Preserving (AP) property, see e.g. [12]. The passage with ε → 0 resembles the low Mach number limit problem, where one observes incompressible behaviour in regions where the Mach number approaches 0. However, in the contrary to the low Mach number limit, in our model the singularity is embedded in the definition of the singular pressure π.
We adapt numerical method from [11] where the Euler system with constant maximal density constraint has been studied. One dimensional version of the Direct method is modified and extended to capture variable density constraint and the viscosity term in the momentum equation. In what follows we focus only on the new elements of our approach, for the detailed description of the other parts we refer to [11] and references therein. Following this rule, we present our technique on the time semi-discrete level, the discretization in space is omitted for brevity.

Discretization scheme
To find an approximate solution to the system (5) we propose a splitting algorithm. Time is discretized by one step finite difference (with fixed time step ∆t) and a finite volume method is used in space. At each time step the set of the equations is decomposed into three parts which are solved subsequently in three sub-steps.

Step 1: Hyperbolic part
The numerical solution of the Euler part of the system follows the strategy presented in [11]. The flux in the mass balance and the singular pressure are treated implicitly: The system (65) is reformulated on a discrete level in terms of singular pressure. Into (65a) we substitute implicit mass flux from (65b) and obtain an elliptic equation for the singular pressure where the right hand side of (66) reads φ( n , * n , u n ) = ρ n − (∆t)div(ρ n u n ) + (∆t) 2 div div(ρ n u n ⊗ u n ) + ∇p n * n .
The singular pressure π ε is computed by solving (66) by means of the Newton method with numerical Jacobian. In the next step we invert singular pressure to get the density. The purpose of this approach is to ensure that the density constrain is satisfied ( ≤ * ), which now follows from the definition of π ε . After the new density is obtained we directly update the momentum. This approach is called the Direct method, see [11,Section 4.1]. The second approach presented in the literature is referred to as the Gauge method [12] that is based on the decomposition of the momentum u = a + ∇ϕ, diva = 0 into a divergence free part a, and the irrotational part ϕ. As reported in [11] the Direct method indicates oscillations of the velocity in congested part, while the Gauge method is diffusive in uncongested region. Since the first method does not introduce any additional numerical dissipation, we adapt it for this work. For detailed description of the space discretization we refer to [11].
We would like to emphasize that (65) is a strictly hyperbolic problem, with characteristic wave speeds λ 1,2 = u ± ∂p ∂( / * ) . By the definition, the Courant-Friedrichs-Lewy (CFL) condition for the explicit part is equal to max(|λ 1,2 |) ≤ σ ∆x ∆t , with the Courant number σ. Splitting for implicit singular and explicit background pressure (65b) provides that CFL condition is satisfied uniformly in ε.

Step 2: Diffusion
For the sake of numerical simulations, we consider (3) in a simplified form S(u) = 2µ∆u. We treat the diffusion term implicitly to avoid additional stability restrictions: The presence of the diffusion term is important from analytical reasons only. In fact, the presented numerical scheme has been designed to solve the Euler system, therefore one can take arbitrary viscosity, such that µ ≥ 0. The viscosity coefficient is fixed independently of ε and small enough to recover compressible/incompressible transition. Equation (67) is discretized in space by cell-centered finite volume scheme.

Step 3: Congestion transport
The transport of the congested density is undoubtedly a main new feature of the model (5) and so of the presented numerical scheme. Having the new velocity u n+1 we compute the congested density as follows * n+1 − * n ∆t + u n+1 ∇ * n = 0, where a cell-centered finite volume scheme together with upwind is used in space.

Numerical results
In this section we present four numerical examples that demonstrate behaviour of the proposed model in one-dimensional periodic setting. As a consequence of finite volume framework the proposed scheme conserves mass. As for domain we take the unit interval with the mesh size ∆x = 10 −3 and the time-step ∆t = 10 −4 . In the following we choose singular pressure parameter ε = 10 −4 with the exponents α = β = 2, (6), and background pressure (2) with the exponent γ = 2, if not stated differently. The test cases are: • Case 1 (constant congestion): • Case 2: (x, 0) = 0.6, * (x, 0) = 0.9 + 0.05(cos(10πx) − cos(6πx) + cos(134πx) + cos(24πx)), Case 1 illustrates shock and rarefaction of the density for constant initial congestion density. The initial value of the congestion density stays the same for all times, due to transport. The congestions and rarefactions are created solely due to opposite initial velocities, exactly as in the analogous case from [11]. In Figure 1 we moreover present the behaviour of unknowns for different values of parameter in the singular pressure: ε ∈ {10 −2 , 10 −4 , 10 −6 }. These numerical results show that the algorithm indeed satisfies the Asymptotic Preserving property. Note that with ε decreasing to zero we approach incompressible limit more thoroughly. However, the exact value of the maximal density constraint can never be reached by the numerically computed density.
Case 2 and Case 3 show the main feature of presented model, namely variable congestion density. For both cases the initial maximal density is set to "smooth hat". In the first of them the initial velocity describes the velocities of two groups of individuals that want to go in the opposite directions. The individuals close to the contact line are more willing to compress (as the congestion density is higher). We see that initially the individuals at the rear of the groups press so intensively onto the front members, so that the whole "hat" is tightly filled. When this happens, we see a similar effect as for elastic collision: part of individuals start to move in the opposite direction to initially intended.
Case 3 describes a situation, when the well organized crowd moving in one direction with the same velocity approaches a barrier ahead, being the group of individuals that move much slower and prefer to keep bigger distances between each other. We observe how the faster individuals behind push the slower group to speed up, by filling the all available gaps between the individuals (this is where the congestion occurs at position x ≈ 0.8). This kind of behaviour could be observed, for example, at airports or in the groups of marathon runners.
Case 4 illustrates shock and rarefaction of the density when the maximal density constrain consists of a sum of cosines (periodic setting) with different frequencies. This example mimics randomness in the individual preferences of the members of population. We observe that congested regions "freeze" maximal density due to the zero velocity, which is consistent with the theoretical prediction. As expected from the properties of the limiting system (1f), for ε 1, the singular pressure π ε is activated only in the congested region. In the theoretical part of the paper, this pressure is merely a nonnegative measure in the limit and our simulations seem to confirm this lack of regularity.
Another interesting feature observed in this case is the traveling wave-like behaviour of the density of the crowd. Note that taking time derivative of (1a) and substituting ∂ t (ρu) from (1b) we obtain wave-like equation for density. We observe this effect on the Figure 4 between time t = 0.25 and t = 0.5 at x = 0.2, where two "crowds" interfere with each other. This leads to reaching the congestion density and propagation of the congestion in the opposite directions.
Looking at the Figure 4 it seems that the scheme for the transport of the congestion density is quite diffusive. The high spatial frequency oscillations present at the beginning are very quickly washed away and there only subsists the small frequency components. Thus, the numerical examples presented above should be treated just as the illustration of the behaviour of solutions to the approximation (5). The thorough numerical discussion as well as the study of twodimensional case is a purpose of our further research.