Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model

We introduce a second order model for traﬃc ﬂow with moving bottlenecks. The model consists of the 2 × 2 Aw-Rascle-Zhang system with a point-wise ﬂow constraint whose trajectory is governed by an ordinary diﬀerential equation. We deﬁne two Riemann solvers, characterize the corresponding invariant domains and propose numerical strategies, which are eﬀective in capturing the non-classical shocks due to the constraint activation.


Introduction
Conservation laws with pointwise unilateral constraints on the flux have been studied widely in the last decade.Their peculiarity and main challenge from the analytical and numerical points of view is the possible presence of a non-classical shock (i.e.violating the classical Kružkov [18] or Lax [19] entropy admissibility conditions) at the constraint location.
Scalar conservation laws with fixed flux constraints have been introduced in [12], where an existence result is established for data with bounded total variation (BV ).Further results have been obtained by [4,11,13].In particular, [4] provides the first description of the numerical treatment of the flux constraint in a finite volume setting.In the articles above, the problem formulation is intended to provide a general mathematical framework to model local constraints in traffic flow, such as toll gates, construction works and traffic lights.The approach has been extended to systems in the case of the Aw-Rascle-Zhang (ARZ) second order model [6,23].[16] provides the definition and analysis of two Riemann solvers, as well as the corresponding numerical strategies to compute approximate solutions by a modified Godunov scheme.Existence results for the ARZ model with fixed constraints are provided in [3,5,17].
More recently, the approach has been extended to moving constraints that, in the framework of traffic flow models, are meant to describe the dynamics of moving bottlenecks caused by large slow moving vehicles like trucks or buses.The resulting model consists of a strongly coupled PDE-ODE system, where the PDE is a constrained scalar conservation law describing the global traffic evolution and the ODE gives the trajectory of the slow vehicle, which in turn is affected by downstream traffic conditions through the speed law.The existence of solutions with BV initial data has been proved in [15], while two numerical strategies have been proposed in [9,14] to compute approximate solutions, which correctly capture the non-classical shock arising when the constraint is enforced.
In this article, we consider the ARZ traffic flow model subject to moving constraints.We define two constrained Riemann solvers, characterize the corresponding invariant domains and propose numerical strategies, which are effective in capturing the non-classical shocks due to the constraint activation.The present contribution is a significant step forward compared to the previous contributions since we propose here a both theoretical and numerical comprehensive treatment of a system of conservation laws with a moving constraint.Note also that one of the two constrained Riemann solvers is non conservative.
The numerical approximation of non-classical shocks is well-known to be very challenging since these discontinuities are very sensitive to the underlying numerical diffusion of the scheme, see for instance [20] for more details.In [7], the authors define a conservative scheme which is based on in-cell discontinuous reconstructions of the non-classical shocks for approximating the solutions of a scalar conservation law with a non convex flux.The striking feature of the strategy is to allow for a perfect control of the numerical diffusion associated with the non-classical discontinuities.More precisely, it allows for the exact computation of such isolated simple waves.In [9], the authors succeeded in extending this approach based on in-cell reconstructions to a constrained but still scalar conservation law.Here, we extend it to the system case of the ARZ constrained traffic flow model, which is much more difficult to deal with since in-cell reconstructions have to be proposed for the two components of the unknowns vector.Note in addition that a non-conservative Riemann solver associated with the ARZ model will be studied in details.At last, let us mention that the extension of the in-cell discontinuous reconstruction strategy to the system case has already been considered for non-constrained problems in the recent paper [2] for the non diffusive computation of classical solutions of the barotropic gas dynamics equations in Lagrangian coordinates, and in [1] for the computation of non-classical solutions associated with a model of elastodynamics.

The ARZ model with moving constraints
Let V and R be respectively the maximal speed and the maximal density of the vehicles on a stretch of road and let α ∈ (0, 1) be the coefficient expressing the reduction of the road capacity at the bus position, so that αR is the maximal density at x = y(t) for which the velocity of the vehicles is v = 0. We denote by F α the upper bound of the flux at the bus position and we consider the following PDE-ODE system: with initial conditions ) The phase space is defined by the domain away from the bus position, and at the bus position.Above, ρ = ρ(t, x) and v = v(t, x) denote respectively the density and the mean velocity of traffic.Note that the quantity w = v + p(ρ), usually referred to as Lagrangian marker, is transported at velocity v and depends on the density through a pressure law p ∈ C 2 ([0, +∞[; [0, +∞[) satisfying the following hypotheses: p (ρ) > 0 for every ρ > 0, (2.3b) Basics.We recall in this paragraph the basic properties of system (2.1a).If we denote the flux function by the eigenvalues of the Jacobian matrix Df in the (ρ, v) coordinates are given by the eigenvectors are and we observe that which means that the first characteristic field is genuinely non-linear and the second is linearly degenerate.We also recall that system (2.1a) belongs to the Temple class [21], and its Riemann invariants are given by The Lax curves through a point (ρ, v) in the (ρ, v)-plane are therefore given by In case of no flux restriction due to the presence of the bus, namely when v ≤ V b , the solution to the Riemann problem associated with (2.1) is expected to be the same as the one associated with (2.1a) (with no constraint) and briefly described above.On the contrary, when the bus acts on the flow, namely when v > V b around the bus position, the general idea to construct the solution to the Riemann problem is similar in the sense that the left and right initial states will be joined by simple waves associated with the two characteristic fields, but one has to check carefully that the phase space D α is preserved, namely that 0 ≤ v ≤ V and v + p(ρ) ≤ p(αR) at the bus position.The first constraint is already satisfied when the bus is not present, so that one has to take care of the second inequality which is also equivalent to Under this form, it is clear that the constraint reads as a constraint on the relative flux at the bus position.In order to go further in that direction, let us determine under which condition the quantity v + p(ρ) equals a constant K such that K ≤ p(αR), or equivalently ) and non-decreasing in K since ρ ≥ 0.Moreover, we have φ(0; K) = 0 for all K. Therefore, the maximal possible value F α of the relative flux corresponds to K = p(αR) and with ρ α such that φ (ρ α ; p(αR)) = 0, that is to say Easy calculations allow to find that the largest admissible flux is given by F α = ρ 2 α p (ρ α ).Therefore, the classical solution will remain admissible provided that the relative flux does not exceed the upper bound F α .This criterion will be the key ingredient to determine the two possible Riemann solutions devised in the next section.As we will see, these solutions necessarily develop non-classical shocks, in addition to the usual simple waves associated with the first and second characteristic fields.To conclude this section, Figure 1 below illustrates the above considerations with the notations w α := p(αR) and w max := p(R).

The Riemann problem
Let (ρ l , v l ) and (ρ r , v r ) be two points in the domain D. We consider the Riemann problem for (2.1) corresponding to the initial data We look for self-similar solutions of (2.1), (3.1).Therefore, we assume that the bus speed will be constant: ẏ(t) = V b for all t > 0. Let I be the set If I = ∅, let (ρ, v) and (ρ 1 , v1 ) be the points defined by These are respectively the points with maximal and minimal density of the Lax curve of the first family passing through (ρ l , v l ) which satisfy the condition (2.1b) on the flux.Moreover, we define the point (ρ 2 , v2 ) as This is the point of maximal density of the Lax curve of the second family passing through (ρ r , v r ) for which (2.1b) is satisfied.Finally we denote by (ρ m , v m ) the middle state of the classical Riemann solver for (2.1a), (3.1a), (3.1b): All these points are depicted in Figure 2. Let RS be the standard Riemann solver for (2.1a), see [6], and let ρ((ρ l , v l ), (ρ r , v r ))(•) and v((ρ l , v l ), (ρ r , v r ))(•) be respectively the ρ and v components of the classical solution RS((ρ l , v l ), (ρ r , v r ))(•).
Proof.Let us consider the plane (ρ, ρv) and let w l = v l +p(ρ l ) be the value of the Riemann invariant w at (ρ l , v l ).In these coordinates the set of points (ρ, v) such that v + p(ρ) = w l is the graph of the function ψ(ρ) := ρw l − ρp(ρ).
This function is strictly concave, by (2.3c).Then the cardinality of the set I is at most 2. Since ), there will be exactly two points (ρ, v) and (ρ 1 , v1 ) belonging to the curve ψ and such that Lemma 2 Let w l ∈ [0, p(R)] be fixed and let ρv = F α + V b ρ be the constraint.Let us consider the function ρ → ψ(ρ) = ρw l − ρp(ρ).Let (ρ σ , v σ ) be a point such that ρ σ v σ = ψ(ρ σ ).Under the hypotheses of Lemma 1, we have Proof.By the hypotheses of Lemma 1, there exists a point (ρ m , v m ) such that Assume by contradiction that ψ(ρ σ ) ≥ F α + V b ρ σ , then by the strict concavity of ψ, we would have which contradicts the definition of (ρ, v).Similarly for ρ σ < ρ1 .

The first Riemann solver RS α 1
Let us introduce the first Riemann solver for the constrained problem (2.1).

Definition 1
The Riemann solver is defined as follows. 1 Note that, due to the Rankine Hugoniot relations, the inequalities above are well defined also along jump discontinuities, being equivalent along the right and left traces of the solution.
The first case refers to a situation in which the traffic is influenced by the bus and the bus travels with its own velocity; in the second case the bus and the traffic do not influence each other; the third case represents a road where the traffic is so slow that the bus travels with the speed of the downstream traffic.
Remark that the solution RS α 1 is conservative for both density and momentum of the vehicles.Moreover, in case 1., the solution given by RS α 1 does not satisfy the Lax entropy condition between the states (ρ, v) and (ρ 1 , v1 ), because we have ρ1 < ρ.Therefore, we say that (ρ, v) and (ρ 1 , v1 ) are connected by a non-classical shock.

The second Riemann solver RS α 2
Let us introduce the second Riemann solver for the constrained problem (2.1).

Definition 2 The Riemann solver
The Riemann solver RS α 2 conserves only the number of the vehicles.Indeed, along the line x = y(t) the Rankine-Hugoniot condition holds for the first flux component, because both (ρ, v) and (ρ 2 , v2 ) belong to the line ρv = F α + V ρ, but not for the second component.

Invariant domains
The study of invariant domains for RS α 1 and RS α 2 follows closely [16,Section 3].First of all, let us recall that the sets are invariant for (2.1a) for any 0 < v 1 < v 2 and 0 < w 1 < w 2 with v 2 < w 2 .Moreover, we define the function which gives the value of the invariant w at the point The following Lemma is a direct consequence of the strict convexity of the function ρ → ρp(ρ), see (2.3c).
Lemma 3 Let us suppose that the hypotheses (2.3) hold.For any , then the Riemann solvers RS α 1 and RS α 2 coincide with the standard Riemann solver RS.Therefore, the domain By the hypotheses we have and for any (ρ l , v l ) and (ρ r , v r ) in D v 1 ,v 2 ,w 1 ,w 2 , the classical solution RS((ρ l , v l ), (ρ r , v r )) satisfies the constraint for x = V b t.Hence the solutions given by the two Riemann solvers RS α 1 and RS α 2 coincide with the classical solution.
The next Theorems characterize the invariant domains of the Riemann solvers RS α 1 and RS α 2 when they are different from the standard Riemann solver.Figure 3 gives two examples of these domains.Let us define the sets The proofs of Theorems 1 and 2 will be developed in the next subsections.
,w 2 is invariant for the standard Riemann solver RS, we have Therefore RS α 1 and RS α 2 coincide with RS and D v 1 ,v 2 ,w 1 ,w 2 is invariant for both, proving points (i) of Theorems 1 and 2.

D
Since v + p(ρ) = v l + p(ρ l ), we find , hence the function h α is well defined in v and its value is Let ṽ ∈ ]V b , +∞[ be the minimum of the function h α which exists by Lemma 3. If Since (ρ 1 , v1 ) satisfies the same hypotheses of (ρ, v), the proof is similar.

Let us suppose that
Hence the right trace of hence the inequality h α (v 2 ) ≥ w 2 must hold.

Proof of Theorem 2 (ii)
The proof is the same of part 2 of Section 4.2.1.
Proof.Assume by contradiction that there exists ṽ We note that 5a).Indeed by definition we have This is in contradiction with the invariance of the domain w 1 ,w 2 be the solution to the system (see Figure 5b) and let (ρ l , v l ) ∈ D v 1 ,v 2 ,w 1 ,w 2 be the solution to The points (ρ l , v l ) and (ρ r , v r ) are connected by the standard Riemann solver with a rarefaction, because v l = v 1 < v 2 = v r .We claim that ρ l ≥ ρ and ρ r < ρ1 .Indeed, let ṽ be the minimum of the function h α .We have v 1 ≤ ṽ, otherwise for every If v ≥ ṽ, then we clearly have v 1 ≤ v.By the definition of (ρ l , v l ), we have In both cases we have v ≥ v l which is equivalent to , which is a contradiction.Therefore v 2 > ṽ.If v1 ≤ ṽ, then clearly v 2 > v1 , while if v1 > ṽ, since This proves the claim.Since the classical solution RS((ρ l , v l ), (ρ r , v r )) does not satisfy the constraint for Lemma 2, the right trace of Therefore (ρ 2 , v2 ) does not belong to D v 1 ,v 2 ,w 1 ,w 2 , in contradiction with the hypothesis of invariance of this domain for the Riemann solver RS α 2 . ρ Representation of the points used in Lemma 5.The invariant domain is the colored area.If, by contradiction, there was ṽ ∈ [v 1 , v 2 ] such that h α (ṽ) < w 1 , the point (ρ 2 , v2 ) would not be in the domain D v1,v2,w1,w2 .
Representation of the points used in the proof of Lemma 6.The invariant domain is the colored area.If, by contradiction, the inequality h α (v 2 ) > w 2 held, the point (ρ 2 , v2 ) would not be in the set D v1,v2,w1,w2 .Proof.Since D v 1 ,v 2 ,w 1 ,w 2 is invariant for RS, we have only to show that the left and right traces of )) coincides with RS((ρ l , v l ), (ρ r , v r )), there is nothing to prove.Otherwise, let (ρ 2 , v2 ) and (ρ, v) be respectively the right and left traces of the solution at We can restate the results that we have obtained until now as follows.

The domain
Proof.The domain D v 1 ,V b ,w 1 ,w 2 and D V b ,v 2 ,w 1 ,w 2 coincide respectively with If the solutions given by the two Riemann solvers RS α 1 and RS α 2 coincide with the classic one, then the result is trivial.Moreover we observe that, for every pair of initial data ((ρ l , v l ), (ρ r , v r )) ∈ (R + × R + ) 2 , the points (ρ 1 , v1 ), (ρ 2 , v2 ) and (ρ, v) belong to D V v 1 ,v 2 ,w 1 ,w 2 .Therefore, if (ρ l , v l ), (ρ r , v r ) ∈ D V v 1 ,v 2 ,w 1 ,w 2 , the solutions given by the Riemann solvers RS α 1 and RS α 2 are contained in Therefore, the solutions given by the Riemann solvers RS α 1 and RS α 2 coincide with the classical one.If v r = V b , then the points (ρ l , v l ) and (ρ r , v r ) are both in Hence the solutions RS α 1 ((ρ l , v l ), (ρ r , v r )) and RS α 2 ((ρ l , v l ), (ρ r , v r )) are entirely contained in one of the two sub-domains and therefore also in we can conclude by Theorem 1 (ii) and Theorem 2 (ii).
, and the standard Riemann solver connects these two points with a rarefaction.If (ρ m , v m ) satisfies the constraint, then the same holds for each point (ρ σ , v σ ) of the rarefaction, because ρ σ ≥ ρ m and Lemma 2. Hence RS α 1 and RS α 2 coincide with the standard Riemann solver and the solutions are contained in D v 1 ,v 2 ,w 1 ,w 2 .If (ρ m , v m ) does not satisfy the constraint, let (ρ * , v * ) be the solution to the system By consistency, we can decompose the solution to the Riemann problem (2.1), (3.1) given by RS α 1 , as 1 ((ρ * , v * ), (ρ r , v r ))(x/t) if x > λ 1 (ρ * , v * )t, and similarly for the second Riemann solver If the classical solution does not satisfy the constraint in x = V b t, then the left trace of RS α i for i = 1, 2 at λ = V b is (ρ, v).The solutions given by the two Riemann solvers RS α i for i = 1, 2 present a rarefaction between (ρ l , v l ) and (ρ * , v * ) and another rarefaction linking (ρ * , v * ) to (ρ, v).Moreover, since D U v 1 ,v 2 ,w 1 ,w 2 and D V v 1 ,v 2 ,w 1 ,w 2 are invariant for both Riemann solvers, the point (ρ * , v * ) belongs to the intersection of the two domains and since This concludes the proofs of Theorems 1 and 2.

Numerical schemes
The aim of this section is to propose a numerical strategy to compute the solutions to the constrained problem (2.1).We will more precisely propose two variants of the strategy which will correspond to the two Riemann solvers RS α 1 and RS α 2 .Importantly, a particular attention will be given to the numerical treatment of the non-classical shocks arising at the bus position since the numerical approximation of such moving discontinuities is well-known to be challenging.We refer the reader to the introduction of this paper for more details.
Let h and k n be positive real numbers representing the increments for space and time discretizations and let us define the mesh points (t n , x j+1/2 ) by x j+1/2 = jh for every j ∈ Z and t n+1 = t n + k n for every n ∈ N.
We divide the x-axis in a sequence {C j } j∈Z of cells such that C j = [x j−1/2 , x j+1/2 ).We also define the cell centers as For every n ∈ N, we aim at constructing a piecewise constant approximation x → u(t n , x) of the conserved variables u = (ρ, z) = (ρ, ρw) = (ρ, ρ(v + p(ρ))) given by For completeness, we recall that the classical Godunov scheme can be written as for all j ∈ Z, where This classical Godunov scheme will be used except near the non-classical shocks since such discontinuities require specific treatments, which are described in the next two sections for both solvers RS α 1 and RS α 2 .Note that above, and in the following, we denote by RS, RS α 1 and RS α 2 the Riemann solvers defined on the conserved variables (ρ, z).At last, the proposed numerical approximation of the bus trajectory is described in Appendix A and we emphasize that it is exact when the bus velocity is given by V b .

The constrained Godunov scheme for RS α 1
First of all, we define for every j ∈ Z the approximation u 0 j of the initial data u 0 on the cell C j as the mean In this section, we consider the Riemann solver RS α 1 and our aim is to modify the classical Godunov scheme near the non-classical shocks that may develop at the moving constraint location.We follow the approach proposed in [1,7,9],x where an in-cell discontinuous reconstruction strategy is developed and allows a sharp computation of the non-classical shocks.By sharp, we mean that such isolated discontinuities are exactly computed with only one point of numerical diffusion, the value of which corresponds to the mean value of the exact solution in the corresponding cell.It is clearly the best result achievable by a conservative scheme.
Let us first consider the case in which the bus is not influenced by the preceding vehicles, so that its position at the time t n is y n := y(t n ) = y 0 + V b t n ∈ C m for some m ∈ Z and the value of the approximate solution at y n is u n m = (ρ n m , z n m ).Following [9] and referring to the definition of the Riemann solver RS α 1 , a first idea to detect whether a non-classical shock is expected to develop around y n or not is to check if the inequality holds.If this inequality actually holds true, one is naturally tempted to consider that a non-classical shock arises as the solution given by RS α 1 to the Riemann problem with initial datum In this case, it is also natural to check whether also the inequality holds true or not.
Let us assume for a moment that (5.2) and ( 5.3) are relevant to detect whether a non-classical shock occurs at the bus position or not (notice from now on that it will not be the case).The proposed strategy to properly compute the non-classical discontinuity which is expected to occur at the bus position consists in making a reconstruction of this non-classical discontinuity within the m-th cell, and to modify the classical Godunov scheme accordingly.More precisely, the proposed reconstruction will lead to a new definition of the left and right numerical fluxes F n m±1/2 , which we describe below.
We introduce in the m-th cell a left state u n m,l = (ρ n m,l , z n m,l ) and a right state u n m,r = (ρ n m,r , z n m,r ) as u n m,l = û and u n m,r = ǔ1 , where û and ǔ1 are the points defined by (3.2), (3.3).Then we replace u n m with the function u n rec = (ρ n rec , z n rec ) defined by where we have used the two points Conditions (5.4) ensure that the reconstructed discontinuity is located in the cell C m ; see Figure 6.
Since the Riemann solver RS α 1 is conservative with respect to u = (ρ, z), we also aim at preserving this property in the process of reconstruction.With this is mind, we require the average value of the reconstructed discontinuity to be equal to u n m , which writes Therefore, we get the following definition of the two constants d ρ,n m and d z,n m , namely Note that these two constants can be different in general, which means that the reconstructed discontinuities of the two components of u can be located at different positions.
The non-classical shock travels with the bus speed V b > 0. As a consequence, if we denote ∆t ρ m+1/2 and ∆t z m+1/2 respectively the time needed by the ρ and the z component of the discontinuity to reach the interface x m+1/2 , we have Therefore, considering the proposed reconstruction and since the waves emerging from the Riemann problem associated with u n m,r and u n m+1 at the interface x m+1/2 propagate to the right, the flux at interface x m+1/2 equals f (u n m,r ) until t n + ∆t ρ,z m+1/2 (depending on the ρ and z components) and f (u n m,l ) afterwards.One then proposes to replace the classical Godunov flux F n m+1/2 by a new numerical flux, whose components are denoted by F 1 (u n m , u n m+1 ) and F 2 (u n m , u n m+1 ) and given by It is also natural to modify the left flux at interface x m−1/2 according to the left state u n m−1 and the new right state u n m,l = û, i.e. we propose to replace F m−1/2 with F (u n m−1 , u n m,l ), which preserves the consistency of the Godunov method.This concludes the description of the reconstruction method.
Let us now briefly go back to the proposed procedure (5.2), (5.3) to detect whether a nonclassical shock appears or not.It turns out that as it stands, it may introduce undesirable oscillations in some cases.The example showed in Figure 7 is one of these situations: at some iterations the exact solution does not satisfy the constraint, while the numerical solution does, leading to an incorrect computation of the non-classical shock and mild oscillations.In view of this example, we propose to remove the first condition (5.2) and to keep only the inequality (5.3) as necessary to enforce the reconstruction procedure, together with (5.4) of course.(5.7) Let us now establish a nice property, which states that the scheme is able to exactly compute isolated non-classical discontinuities and thus justifies the proposed reconstruction procedure.In particular, this property also makes clear why the proposed numerical scheme performs so well in computing the non-classical discontinuities, since such moving shocks are proved to be computed with no extra numerical diffusion by construction.
Proposition 3 The scheme (5.1)-(5.7) is exact when the Riemann initial datum is made of a single nonclassical shock between the left state û and the right state ǔ1 and provided that the bus position is computed exactly (which is the case with the front tracking method introduced in Appendix A).In particular, let d ρ,n m and d z,n m be the two constants defined in (5.5).If u n m−1 = û, u n m+1 = ǔ1 and there exists γ ∈ [0, 1] such that ) Proof.Let us start with an initial condition such that u 0 j = û if j < m and u 0 j = ǔ1 if j ≥ m.It is clear on the one hand that (5.8) holds true with n = 0 and γ = 0, and on the other hand that u exa (t n , x) dx for all x ∈ Z (5.9) with n = 0 and u exa = u 0 .In order to prove that the scheme is exact for such an initial datum, it is thus sufficient to assume that (5.8) and (5.9) hold true for a given time t n , and to show that x j−1/2 u exa (t n+1 , x) dx for all x ∈ Z. (5.10) Note that since the exact solution is a discontinuity propagating with velocity V b , (5.10) implies immediately the validity of a convex combination like (5.8) at time t n+1 .Let us first notice that F n m−1/2 = F (û), and that (5.8) holds true at time t n+1 for j = m, m + 1, since our scheme is equivalent to the usual Godunov's scheme and we are in the constant regions of the exact solution.Next, since condition (5.8) holds, we have Hence we find This means that the reconstructed discontinuity is located at the same position as the discontinuity of the exact solution.As a consequence, the numerical flux F n m+1/2 defined by (5.7) is also nothing but the exact flux passing through the interface x m+1/2 so that u n+1 m and u n+1 m+1 also exactly coincide with the average of the exact solution in the corresponding cell.Which concludes the proof.
Finally, note that if which corresponds to the second case of Definition 1, the bus and the vehicles do not influence each other.Hence, the bus position at the time t n+1 is y n+1 = y n + V b k n and the solution of the ARZ system can be computed with the standard Godunov's method.In the third case of Definition 1, the solution of the ARZ system can also be computed with the standard Godunov's method but the bus position will be updated taking into account the flow velocity v using the front tracking technique described in Appendix A.

5.2
The constrained Godunov schemes for RS α 2 Let V n be the bus speed at time t n .The Riemann solvers RS α 1 and RS α 2 give the same solution whenever the constraint is satisfied.Hence, if we can apply to RS α 2 the techniques described for RS α 1 .When the constraint is enforced, something special has to be done to properly compute the non-classical shocks, in the spirit of the method proposed in the previous section for RS α  1 .Note that the non-classical shocks now join û and ǔ2 instead of ǔ1 .Actually, we will follow exactly the same approach based on in-cell discontinuous reconstructions but since the Riemann solver RS α 2 is not conservative on the second component of u (unlike RS α 1 ), the scheme will be based on the evolution of ρ and v, instead of ρ and z = ρw.Results not reported here show that considering the evolution of the conservative variable u = (ρ, ρw) does not lead to a good approximation of the non-classical shocks joining û to ǔ2 .
First of all, we introduce the vector of non conservative variables v = (ρ, v) and the natural changes of variables v = v(u) and u = u(v), with obvious notations.We then define for every j ∈ Z the approximation u 0 j of the initial data u 0 on the cell C j as follows, Again, let y n := y(t n ) be the bus position at time t n and let m ∈ Z be such that y n ∈ C m .If the Riemann solver RS α 2 does not give the classical solution, a non-classical shock is expected to appear at x = y(t) and arise as the solution given by RS α 2 to the Riemann problem with initial datum As before, we will make a reconstruction of the discontinuity if the inequality holds.In this case we modify the Godunov's scheme as follows.
We introduce in the m-th cell a left state u n m,l = (ρ n m,l , z n m,l ) and a right state u n m,r = (ρ n m,r , z n m,r ) defined by u n m,l = û and u n m,r = ǔ2 , where û and ǔ2 are given by (3.2), (3.3).We then replace u n m by the function u n rec = (ρ n rec , z n rec ) defined by means of ρ n rec and v n rec (recall that z = ρ(v + p(ρ))) where we have used the two points (5.13) We emphasize here that our in-cell reconstruction procedure is fully based on the ρ and v variables (z = ρw just follows from the relation z = ρ(v + p(ρ)) and therefore generally does not just consists in a single discontinuity) while the positions of the discontinuities in ρ and v are defined in such way that the average values of the reconstructed discontinuities equal ρ n m and v n m respectively.More precisely, solving (5.13) for d ρ,n m and d v,n m , we find (5.14) Clearly, the conditions d ρ,n m ∈ [0, 1] and d v,n m ∈ [0, 1] are necessary to reconstruct the discontinuity in the cell C m .Notice again that these two constants are in general different.
As we did in the previous section for the first Riemann solver, we naturally assume that the discontinuities in (5.12) propagate at the same speed of the non-classical shock.Therefore their positions at time t n+1 are Depending on the values of V n and k n , the ρ and v reconstructed discontinuities at time t n+1 may be located in either the m-th or (m + 1)-th cell, that is to say x n+1,ρ m and x n+1,v m can be less or larger than x m+1/2 .Since our objective is to define the solution at time t n+1 by averaging ρ and the non conservative variable v in the vicinity of the bus position (and to carry on using the conservative variables ρ and z = ρw "away" from the bus position), one is led to distinguish between the following two cases in order to make the ρ and v variables evolve (note that ρ is a conservative variable so that we can use a flux formulation).
Reconstruction of the ρ-component.Since the first variable ρ is still conserved by RS α 2 , we can apply the same strategy used in Section 5.1 for RS and we apply the usual conservation formula Reconstruction of the v-component.We distinguish two cases: (i) If x n+1,v m < x m+1/2 (see Figure 8a), the average of v in the m-th cell at time t n+1 is given by The corresponding value of z is and to update the solution in the (m + 1)-th cell, we use the right trace value ǔ2 to compute the numerical flux at the interface x m+1/2 , i.e.
The other cells are updated in the usual conservative way.
(a) Notations used in case 1v.Like in the previous section, the next proposition states that if the initial datum is a non-classical shock, then the solution given by the in-cell discontinuous reconstruction method on density and velocity is the non-classical shock itself.Proposition 4 The scheme (5.1), (5.11)-(5.18) is exact when the Riemann initial datum is made of a single non-classical shock between the left state û and the right state ǔ2 and provided that the bus position is computed exactly (which is the case with the front tracking method introduced in Appendix A).In particular, let d ρ,n m and d v,n m be the two constants defined in (5.14).If v n m−1 = v, v n m+1 = v2 and there exists γ ∈ [0, 1] such that then d ρ,n m = d v,n m = γ.
Proof.The proof follows the same steps as the one associated with RS α 1 , but is adapted to the non conservative treatment of the v variable.Let us consider an initial condition such that v 0 j = v if j < m and v 0 j = v2 if j ≥ m, so that (5.19) holds true with n = 0 and γ = 0, and x j−1/2 v exa (t n , x) dx for all x ∈ Z (5.20) for n = 0 and v exa = v 0 .In order to prove that the scheme is exact for such an isolated non-classical shock, let us assume that (5.19) and (5.20) hold true for a given time t n , and let us show that x j−1/2 v exa (t n+1 , x) dx for all x ∈ Z. (5.21) 2. Let us suppose that ρ n m−1 ≥ ρ int m−1/2 , so that u n m−1 and u int m−1/2 are connected by a rarefaction wave centered at x m−1/2 .Let us recall that the bus speed at the instant t n is V n = min(V b , v n m ).
Proposition 6 If V n < λ 1 (ρ int m−1/2 , v int m−1/2 ), so that an interaction occurs between the bus and the rarefaction wave centered in x m−1/2 , then the bus speed before the interaction is V n = V b and the bus keeps its maximal speed during all the interaction.Proof.We have V n = v n m if and only if v n m ≤ V b , which is in contradiction with the assumption V n < λ 1 (ρ int m−1/2 , v int m−1/2 ).Therefore we must have V n = V b > v n m .Let (ρ, v) be the density and speed values at the last interaction point of the bus with the rarefaction: we must have λ 1 (ρ, v) = min{V b , v}.
Since v > λ 1 (ρ, v), when the interaction occurs the bus cannot take the speed of the vehicles but it keeps its maximal speed V b .Moreover the bus interacts only with the points (ρ, v) on the rarefaction wave such that Since for all these points we have v > v > V b , the bus keeps its maximal speed V b during the whole interaction.
Therefore in both cases the bus speed remains constant and the bus position at time t n+1 is given by y n+1 = y n + k n min(V b , v n m ).
Flux representation in the phase plane.
Flux representation in the bus reference frame.

Figure 1 :
Figure 1: Representation of the phase plane in the fixed and in the bus reference frame.

Figure 2 :
Figure 2: Notations used in the definition of the Riemann solvers.

Figure 3 :
Figure 3: Example of an invariant domain (the shaded area) for RS α 1 (left) and RS α 2 (right), for v 1 > V b .

Figure 6 :
Figure 6: An example of a discontinuity reconstruction for the Riemann solver RS α 1 .
Notations used in case 2v.
At last, let us briefly recall that the solution to the Riemann problem associated with the usual system (2.1a) (with no constraint) is given by two simple waves separated by a single intermediate state.The first simple wave is a shock wave or a rarefaction wave associated with the first characteristic field and such that the intermediate state belongs to the first Lax curve L 1 passing through the left initial state, while the second one is a contact discontinuity associated with the second characteristic field and such that the intermediate state belongs to the second Lax curve L 2 passing through the right initial state.