A reactive dynamic continuum user equilibrium model for bi-directional pedestrian flows

Acta Mathematica Scientia 2009,29B(6):1541–1555 http://actams.wipm.ac.cn

A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL FOR BI-DIRECTIONAL PEDESTRIAN FLOWS∗ Dedicated to Professor James Glimm on the occasion of his 75th birthday Yanqun Jiang1 Tao Xiong1 S.C. Wong2† Chi-Wang Shu3 Mengping Zhang1 Peng Zhang4 William H.K. Lam5 1.Department of Mathematics, University of Science and Technology of China, Hefei 230026, China 2.Department of Civil Engineering, The University of Hong Kong, Hong Kong, China 3.Division of Applied Mathematics, Brown University, Providence, RI 02912, U.S.A. 4.Shanghai Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, China 5.Department of Civil and Structural Engineering, The Hong Kong Polytechnic University, Hong Kong, China E-mail: [email protected]; [email protected]; [email protected]; [email protected]; [email protected]; [email protected]; [email protected]

Abstract In this paper, a reactive dynamic user equilibrium model is extended to simulate two groups of pedestrians traveling on crossing paths in a continuous walking facility. Each group makes path choices to minimize the travel cost to its destination in a reactive manner based on instantaneous information. The model consists of a conservation law equation coupled with an Eikonal-type equation for each group. The velocity-density relationship of pedestrian movement is obtained via an experimental method. The model is solved using a finite volume method for the conservation law equation and a fast-marching method for the Eikonal-type equation on unstructured grids. The numerical results verify the rationality of the model and the validity of the numerical method. Based on this continuum model, a number of results, e.g., the formation of strips or moving clusters composed of pedestrians walking to the same destination, are also observed. Key words pedestrian flows; conservation law; Eikonal-type equation; density-velocity relationship; finite volume method; fast marching method; unstructured grids 2000 MR Subject Classification

76M12

∗ Received October 22, 2009. The work described in this paper was jointly supported by grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (HKU 7183/06E), the University of Hong Kong (10207394) and the National Natural Science Foundation of China (70629001 and 10771134). † Corresponding author: S.C. Wong.

1542

1

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

Introduction

Walking is a sustainable mode of transportation that is environmentally friendly, highly flexible, and healthy, which also serves as an interface with motorized transportation facilities [1, 2]. However, there are inherent dangers associated with large-scale public walking facilities, such as transfer stations, airports, and train and subway stations. Many tragic crushing incidents have occurred around the world due to overcrowding in such facilities. For example, in Sheffield, 96 Liverpool fans were killed at the Hillsborough football grounds when police open gates to alleviate crowding, and the Mena Valley disaster saw 251 Muslims killed in a stampede during the annual Hajj pilgrimage [3–5]. In most such instances, the fatalities and injuries are caused by the behavior of the crowd itself, and this issue is receiving increased attention from scientists. Many crowd-flow models (e.g., [3–19]) have been developed to investigate the complex dynamics and behavior of pedestrians in various situations. A number of quantitative observations of pedestrian motion characteristics, including selforganization phenomena in pedestrian crossing streams and oscillations of the passing direction at bottlenecks [4, 7–9, 12, 20–22], have been successfully depicted by computer simulation using micro-simulation models and by empirical or experimental research that makes use of video analysis. In microscopic models, each individual is represented separately [23, 24]. These include the cellular automation model [7, 8, 17, 25, 26], the social force mdoel [10, 11, 27, 28], and the lattice gas model [6, 29–31]. Such an approach captures each pedestrian’s unique situation well, such as his or her route choice, which are particularly well suited for use with small crowds. Macroscopic models, in contrast, focus on the overall behavior of pedestrian flows and are more applicable to investigations of extremely large crowds, especially when examining aspects of motion in which individual differences are less important. The present research investigates the flux or density distribution and the impact factor of the speed effect and path choices between two interactive streams of pedestrians, i.e., a bidirectional pedestrian flow, traveling on crossing paths in a continuous walking facility. Based on the observation of groups of pedestrians with fluid-like properties moving in a continuous space, a continuum model is adopted to simulate pedestrian movement. One prior study [15] provided a systematic framework for a dynamic macroscopic model without considering the user equilibrium concept. In contrast to the predictive user equilibrium model [13, 14], the reactive dynamic user equilibrium model describes the movement of pedestrians who do not have predictive information when they are making a path-choice decision [19, 32]. This means that pedestrians have to rely on the instantaneous information available to them and make their choices in a reactive manner to minimize the walking cost to their destination. The model presented in this paper is designed for bi-directional pedestrian flows and consists of a two-dimensional (2D) scalar hyperbolic conservation law equation coupled with an Eikonal-type equation for each group. It is extended from a single pedestrian-type reactive dynamic user equilibrium model [16, 19, 32]. The model system is solved using a finite volume method (FVM) and a fast marching method (FMM) on an unstructured grid. The former is simple and capable of flexibly dealing with body boundaries with complicated structures. In addition, it provides conservative schemes with shock-capturing properties and is extensively applied to solve the Euler and Navier-Stokes equations (e.g., [33–35]). The FMM, which is

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1543

an ordered upwind method, is appropriate for boundary value problems and is very effective in solving the continuous Eikonal problem, which is discussed in detail in a number of other studies [36, 37]. By simulating this continuum model, a number of results, e.g., the formation of strips or moving clusters composed of pedestrians walking to the same destination, that are obtained with microscopic models and experimental investigations [22], are also observed here. The rest of the paper is organized as follows. In Section 2, the governing models of pedestrian flow are described. The numerical methods for solving the models are given in Section 3. Section 4 demonstrates a numerical experiment and presents related analysis. Section 5 concludes.

2

Problem Formulation

In this section, a series of formulas for bi-directional pedestrian flows are given. For ease of reference, the following nomenclature is provided. (1) Ω represents a 2D continuous walking facility (in m2 ). (2) Group c represents the c-th type of pedestrians marching toward the c-th destination of the walking facility. (3) Γ is the boundary of Ω (in m), Γco the original segment from which Group c enters into the walking facility (in m), Γcd the destination segment from which Group c leaves the walking facility (in m), and Γch the wall segment from which nobody in Group c is allowed to enter or leave the walking facility (in m). (4) T is the time horizon (in s). (5) q c (x, y, t) expresses the number of pedestrians who cross a unit width of Γco and describes the time-varying demand of Group c (in ped/m/s). (6) ρc (x, y, t) is a time-varying function that denotes the pedestrian density of Group c (in ped/m2 ). (7) Fc (x, y, t) := (f1c (x, y, t), f2c (x, y, t)) represents the flow vector for Group c, where f1c (x, y, t) is the flow flux in the x-direction and f2c (x, y, t) is the flow flux in the y-direction (in ped/m/s). By default, Fc (x, y, t) · n = 0, where n denotes the unit outward normal of Γch . (8) v c (x, y, t) is the local walking speed of Group c (in m/s). (9) τ c (x, y, t) is the local cost of Group c (in s). (10) Φc (x, y, t) is the instantaneous walking time potential of Group c from location (x, y) to Γcd , at which Φc = 0 (in s).   Finally, Γ = Γco Γcd Γch , c ∈ {a, b}, (x, y) ∈ Ω and t ∈ T . 2.1 Flow conservation equation Similar to most physical systems, the density, velocity and flow of pedestrians follow the physical principle of conservation. The conservation of mass for bi-directional pedestrian flows is given by ρct (x, y, t) + ∇ · Fc (x, y, t) = 0, (2.1) c

∂f c (x,y,t)

∂f c (x,y,t)

where ρct (x, y, t) = ∂ρ (x,y,t) and ∇ · Fc (x, y, t) = 1 ∂x + 2 ∂y . ∂t 2.2 Walking speed function The local walking speeds of the two groups depend on the densities and traveling directions of the pedestrians in the walking facility and are determined by (see [38])     v a (˜ (2.2) ρ, Ψ) = vf exp − α(ρa + ρb )2 exp − β(1 − cos(Ψ))(ρb )2 ,     v b (˜ (2.3) ρ, Ψ) = vf exp − α(ρa + ρb )2 exp − β(1 − cos(Ψ))(ρa )2 , where ρ˜ := {ρa (x, y, t), ρb (x, y, t)}. Here, vf is the free-flow walking speed of pedestrians, α, β, are the model parameters, and Ψ := Ψ(x, y, t) is the intersecting angle between the two

1544

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

pedestrian streams of Group a and Group b, respectively, at location (x, y), as shown in Figure 1.

Fig.1

Angle Ψ between Streams a and b

For the case in which there is only single stream a, i.e., ρb = 0, (2.2) and (2.3) can be formulated as   v a (˜ ρ, Ψ) = vf exp − α(ρa )2 ,     v b (˜ ρ, Ψ) = vf exp − α(ρa )2 exp − β(1 − cos Ψ)(ρa )2 . However, v b (˜ ρ, Ψ) in this situation is, effectively, only a direction-dependent potential speed if there is potentially a single pedestrian trying to walk in direction Ψ against Stream a. Therefore, the actual speed of Stream b is not relevant. Also, the speed is identical to v a (˜ ρ, Ψ) if this potential pedestrian happens to walk parallel to Stream a (i.e., Ψ = 0). 2.3 Local cost function The local cost, τ c , as a function of location and time is assumed to depend on walking speed alone and can be specified as τ c (x, y, t) : = τ (v c (˜ ρ, Ψ)) =

1 , v c (˜ ρ, Ψ)

(2.4)

which describes the walking time per unit distance of movement incurred by a pedestrian in Group c in the walking facility. 2.4 Definitional relationship For each Group c ∈ {a, b}, along the direction of the flow vector (or movement), flow intensity, which is determined as the norm of the pedestrian flow, is equal to the product of speed and density,  Fc (x, y, t) = v c (˜ ρ, Ψ)ρc (x, y, t), (2.5) where

1   Fc (x, y, t) = f1c (x, y, t)2 + f2c (x, y, t)2 2 .

2.5

Path choice constraints It is assumed that a pedestrian in each group at location (x, y) ∈ Ω will choose a path that minimizes his or her travel cost to the destination, based on the instantaneous travel cost information that is available at the time of making a decision. Then, we have τ c (x, y, t)

Fc (x, y, t) + ∇Φc (x, y, t) = 0,  Fc (x, y, t) 

where ∇Φc =

 ∂Φc ∂Φc  , . ∂x ∂y

(2.6)

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1545

Based on (2.4)–(2.6), Fc (x, y, t) = −

ρ, Ψ)ρc (x, y, t) v c (˜ ∇Φc (x, y, t) τ c (x, y, t)

= −ρc (x, y, t)(v c (˜ ρ, Ψ))2 ∇Φc (x, y, t),

(2.7)

and  ∇Φc (x, y, t) = τ c (x, y, t).

(2.8)

From (2.7), assuming that ρc (x, y, t) = 0, the following property is obtained [19, 39–43]. Fc (x, y, t)//(−∇Φc (x, y, t)),

(2.9)

where “//” means that the two vectors, Fc (x, y, t) and ∇Φc (x, y, t), are parallel and pointing in the opposite direction. Furthermore, (2.8) can be shown as follows. For any used path p, from location (x, y) ∈ Ω to destination (x0 , y0 ) ∈ Γcd , if we integrate the local travel cost along the path, then the total cost incurred by pedestrians can be obtained by    Fc (x, y, t) τ¯pc (x, y, t) = τ c (x, y, t)ds = τ c (x, y, t) · ds = − ∇Φc (x, y, t) · ds  Fc (x, y, t)  p p p = −(Φc (x0 , y0 , t) − Φc (x, y, t)) = Φc (x, y, t),

(2.10)

Fc (x,y,t) Fc (x,y,t) is ∈ Γcd , t ∈ T.

a unit vector that is parallel to ds along the path, using (2.6)–(2.9), the fact that and Φc (x0 , y0 , t) = 0, ∀(x0 , y0 ) Hence, the total instantaneous travel cost at time t ∈ T is independent of the used paths. In contrast, for an unused path p˜, between location (x, y) ∈ Ω and destination (x0 , y0 ) ∈ Γcd , the total walking time incurred by the pedestrian is    Fc (x, y, t) c c c · ds = − ∇Φc (x, y, t) · ds τ¯p˜ (x, y, t) = τ (x, y, t)ds ≥ τ (x, y, t)  Fc (x, y, t)  p˜ p˜ p˜ = −(Φc (x0 , y0 , t) − Φc (x, y, t)) = Φc (x, y, t),

(2.11)

using (2.6)–(2.9). The inequality in the foregoing derivation occurs because of some segments along path Fc (x,y,t) Fc (x,y,t) p˜ whose normal vectors F c (x,y,t) are not parallel, ds; thus, ds > ( Fc (x,y,t) ) · ds for those segments. Hence, for any unused paths, the total instantaneous travel cost is greater than or equal to that of the used paths. In this way, the model guarantees that pedestrians choose their paths in the walking facility in a user-optimal manner with respect to the instantaneous travel information [44], that is, the pedestrian flow conditions at the time. The model can now be formulated as the following set of differential equations. ⎧ c c ⎪ ⎪ ⎨ ρt (x, y, t) + ∇ · F (x, y, t) = 0, (2.12) Fc (x, y, t) = −ρc (x, y, t)(v c (˜ ρ, Ψ))2 ∇Φc (x, y, t), ⎪ ⎪ ⎩  ∇Φc (x, y, t) = τ c (x, y, t), where τ c (x, y, t) are determined by (2.2)–(2.4), respectively. Note that angle Ψ in (2.2) and (2.3) can be obtained by two potential functions due to (2.7). Then, we have cos(Ψ) =

∇Φa (x, y, t) ∇Φb (x, y, t) · . a  ∇Φ (x, y, t)   ∇Φb (x, y, t) 

1546

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

The system (2.12) is subject to the following initial boundary conditions. Fc (x, y, t) · n(x, y) = q c (x, y, t), ρ (x, y, 0) = c

ρc0 (x, y),

Φ (x, y, t) = 0, c

∀(x, y) ∈ Γco ,

∀(x, y) ∈ Ω, ∀(x, y) ∈

Γcd ,

(2.13) (2.14) (2.15)

where n(x, y) is a unit normal vector that points to the domain boundary, and ρc0 (x, y) is the initial pedestrian density for each group.

3

Solution Algorithms

Here, the cell-centered FVM for an approximation of the 2D conservation law equation and the FMM for the Eikonal-type problem on unstructured meshes, respectively, are introduced to solve the system (2.12). 3.1 Finite Volume Spatial Discretization  Let T h := Ti be an acute triangulation of computational domain Ω with Np nodes, Ti ∈T h

Nt cells and Nf faces. The size of the discretization is defined by h := superTi ∈T h {hTi }, where hTi is the exterior diameter of triangular cell Ti . We label the neighboring triangles of Ti as Tik , with a counter-clockwise convention as sketched in Figure 2.

Fig.2

A typical unstructured triangular grid

A cell-centered FVM is considered, i.e., each cell represents a control volume (CV). The unknowns are stored in the geometric center of the cells. Integrating (2.1) over Ti and applying the Gauss’s theorem yields  ∂ρc dxdy + Fc · nds = 0, ∀ Ti ∈ T h , (3.16) Ti ∂t ∂Ti where ∂Ti denotes three edges of Ti , and n is the outer unit normal pointing to ∂Ti . Assuming that ρci is the average quantity stored in the cell center of Ti , we obtain the following formula. ∂ρci 1 Fc · nds = 0, ∀ Ti ∈ T h , (3.17) + ∂t Ai ∂Ti where Ai refers to the area of Ti . To obtain the semi-discrete form of (3.17), the boundary integral term needs to be evaluated. Let lik denote the k-th edge of Ti and nik be the corresponding normal vector, k = 1, 2, 3. The integral expression is discretized by summing the flux vectors over each edge of the triangular cell. Utilizing a mid-point quadrature formula generates 3

 c  ˜ ik · nik F Fc · nds = |lik |, ∀ Ti ∈ T h , (3.18) m ∂Ti

ik

k=1

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1547

˜ c is the numerical flux through the middle point mik of lik . where F ik

Fig.3

Schematic diagram of two adjacent triangular cells

To obtain the upwind mechanism, high-order accuracy and essential stability of the numerical scheme, a splitting into the normal flux function Fc · n at location Xjk = (xjk , yjk ), k = 1, 2, 3, 4 (see Figure 3) is adopted. For instance, Lax-Friedrichs splitting [45] is used. (Fc · n)± =

1 [(F · n)(ρc ) ± λρc ], 2

(3.19)

where λ := v c (˜ ρ, Ψ) is the upper bound of the spectrum of the Jacobi of the normal flow function. Here, Xjk , k = 2, 3 is the barycentric coordinate of the triangles, and Xjk , k = 1, 4 is the vertex of the triangles. The unknown concentration at a vertex in the mesh is obtained via the inverse distance weighting interpolation method [46]. The normal numerical flux across the cell boundary l between two neighboring cells can be written as ˜ c · n)m = (Fc · n)L + (Fc · n)R . (F (3.20) For a first-order (1st) scheme, let (Fc · n)L = [F · n(ρcj2 )]+ and (Fc · n)R = [F · n(ρcj3 )]− by assuming constant variables in each cell. Here it is postulated that the variables are the adjacent piecewise constant values in each cell. Applying a MUSCL reconstruction technique yields a second-order (2nd) spatial accuracy scheme [47]. The detailed scheme is + + c c (Fc · n)L = G1 ((Fc · n)+ j1 , (F · n)j2 , (F · n)j3 ),

(F · n)R = G2 ((F · c

c

c n)− j2 , (F

·

c n)− j3 , (F

·

n)− j4 ),

(3.21) (3.22)

with 1 G1 (q1 , q2 , q3 ) = q2 + [(1 − )Δ− + (1 + )Δ+ ]q2 , 4 1 G2 (q2 , q3 , q4 ) = q3 − [(1 − )Δ+ + (1 + )Δ− ]q3 , 4

(3.23) (3.24)

c ± + − where (Fc · n)± jk = [(F · n)(ρjk )] , k = 1, 2, 3, 4, and Δ , Δ indicate the forward and backward difference operators, respectively, e.g., Δ+ qi = qi+1 − qi , Δ− qi = qi − qi−1 , i = 2, 3. The value of  is 13 in this paper. Note that for boundary CV, the values of ρc at Xj3 , Xj4 are obtained by extrapolation from the inside computational domain. To avoid great change in the mesh size of neighboring cells, Δ+ q2 in (3.23) and Δ− q3 in (3.24) are modified as follows.

˜ + q2 = Δ

2l1 ˜ − q3 = 2l2 (q3 − q2 ), (q3 − q2 ), Δ l1 + l2 l1 + l2

1548

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

where li , i = 1, 2 is the length from the centroid of the left or right cell to the mid-point of the edge l. Furthermore, a slope limiter functions L(θ− , θ+ ) can be imposed on (3.23) and (3.24) to suppress spurious oscillations in the process of computation. The two formulas are rewritten as L(θ− , θ+ ) ˜ + ]q2 , [(1 − L(θ− , θ+ ))Δ− + (1 + L(θ− , θ+ ))Δ 4 L(θ− , θ+ ) ˜ − ]q3 . [(1 − L(θ− , θ+ ))Δ+ + (1 + L(θ− , θ+ ))Δ G2 (q2 , q3 , q4 ) = q3 − 4

G1 (q1 , q2 , q3 ) = q2 +

(3.25) (3.26)

A double-parameter minmod limiter [48] is adopted due to its robustness. L(θ− , θ+ ) = minmod(1, θ− ) + minmod(1, θ+ ) − 1. q4 −q3 + 1 where θ− = qq23 −q −q2 , θ = q3 −q2 . 3.2 Calculation of ∇Φc on CV surfaces In (2.7), the values of ∇Φc at each CV boundary are unknown. Assuming that Φc is piece-wise linear on each triangle, it can be written as

Φc (x, y)|Ti =

3

ˆk Φck , ∀Ti ∈ T h , N

(3.27)

k=1

ˆk , k = 1, 2, 3 is a Lagrangian interpolation basis function on the triangular meshes and where N c Φk is the value of Φc at the k-th vertex of triangle Ti . Then, ∇Φc is piece-wise constant on each triangle, and 3

ˆk Φc , ∀Ti ∈ T h . ∇Φc (x, y)|Ti = ∇N (3.28) k k=1

The approximate value of Φc at each vertex of mesh is calculated by solving the nonlinear Eikonal equation. Several numerical methods can be used for this purpose, such as the group marching method (GMM), the FMM or the fast sweeping method (FSM), among others. The FMM is chosen, as its yields are consistent, accurate and highly efficient, based on entropysatisfying upwind schemes and fast sorting techniques. Here a brief review of the upwind finite difference discretization for the Eikonal equation is given. Consider simplex X0 X1 X2 in the triangular meshes. We assume that Φc at X0 is to be computed and that Φcl , ∇Φcl , l = 1, 2 at X1 and X2 are given or have already   been N1 X0 −Xl computed. Define unit vectors Nl = |X0 −Xl | , l = 1, 2 and a 2×2 nonsingular matrix N = N2 (X0 , X1 and X2 are not lined up). The directional derivatives in the N1 and N2 directions are approximated by Φc0 − Φcl , l = 1, 2, 1st formula, |X0 − Xl | Φc − Φcl D l Φc = 2 0 − Nl · (∇Φcl ) , l = 1, 2, 2nd formula, |X0 − Xl |

D l Φc =

(3.29) (3.30)

where the superscript “ ” stands for vector or matrix transposition. The approximation directional derivatives are linked to the gradient by DΦc = N · (∇Φc ) + O(hα ),

(3.31)

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1549

 D1 , h = max{|X0 − X1 |, |X0 − X2 |}, and α = 1 or 2, depending on whether where D = D2 (3.29) or (3.30) is used. Substituting ∇Φc , which was obtained by (3.31), into (2.8), we obtain the (quadratic) equation that defines the unknown Φc0 : 

(DΦc ) (N N  )−1 (DΦc ) = (τ0c )2 .

(3.32)

Φc0 can be updated if the upwind criterion is satisfied: the characteristic direction should point to simplex X0 X1 X2 and is equivalent to (N N  )−1 DΦc > 0. Then, the value of Φc0 at X0 is determined by ⎧ ⎨ solution of (3.29) or (3.30) , Φc0 = min ⎩ min {Φcl + |X0 − Xl |τ0c }, l=1,2

(3.33)

if (3.33) holds, others.

(3.34)

With the FMM, the values of Φc at all mesh nodes are considered in an order that is consistent with the way fronts propagate. This leads us to a single-pass algorithm, which is the advantage of this method in solving front propagation problems with a lower degree of computational complexity. Following Ref. [36], the procedure of the algorithm is given as follows. First, the mesh points are divided into three sets, i.e., Accepted, Considered and Far. For example, tag the points on Γco as Accepted (Φc = 0 on Γcd ) and then tag all points adjacent to those in Accepted set as Considered. Finally, tag all other points as Far. The fast-marching algorithm is briefly described as follows. • The point with the lowest value in the Considered set is removed from it and then added to the Accepted set; its unaccepted neighbors are added to the Considered set. • The values of all points in the Considered set are recomputed using (3.34). • Repeat the process until all of the points are in the Accepted set. For more detailed depictions of the algorithm, refer to Refs. [36] and [37]. 3.3 Time discretization To obtain 2nd accuracy in time, the explicit TVD Runge-Kutta time-stepping method [49] is chosen to discretize (3.17) as follows. (1)

(ρc )i

3

= (ρc )ni −

Δt ˜ c (F · n)nmk |lik |, Ai

(3.35)

k=1

(ρc )n+1 = i

3  1 c n Δt ˜ c (1) (ρ )i + (ρc )i − (F · n)(1) mk |lik | . 2 Ai

(3.36)

k=1

Here, time step Δt needs to satisfy the Courant-Friedrichs-Lewy (CFL) condition. For a 2D MUSCL scheme on unstructured grids, the stability condition is  d(Fc · n)  Δt   max   · max |lik | ≤ 1, i,k mini Ai ρc dρc where i = 1, 2, · · · , Nt , k = 1, 2, 3.

1550

4

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

Numerical Example

In this section, we consider a railway platform that is 100m by 50m in size, as shown in Figure 4. Group c enters the platform from segment Γco with a width of 20m and leaves the platform from segment Γcd with a width of 20m.

Fig.4

Overview of the railway platform

The time step is set to be 240s. Initially, the platform is empty. The inlet flow rate, q (x, y, t), on Γco represents time-varying demand and is given as follows. ⎧ t ⎪ × ma , t ∈ [0, 30], ⎪ ⎪ ⎪ 30 ⎪ ⎪ ⎨m , t ∈ [30, 90], a q a (0, y, t) = t − 120 ⎪ ⎪ , t ∈ [90, 120], −ma × ⎪ ⎪ 30 ⎪ ⎪ ⎩ 0, t ∈ [120, 300], c

with ∀y ∈ [15m, 35m], and

q b (x, 0, t) =

⎧ t ⎪ × mb , ⎪ ⎪ ⎪ 30 ⎪ ⎪ ⎨m , b

t − 120 ⎪ ⎪ , −mb × ⎪ ⎪ 30 ⎪ ⎪ ⎩ 0,

t ∈ [0, 30], t ∈ [30, 90], t ∈ [90, 120], t ∈ [120, 300],

with ∀x ∈ [65m, 85m], where ma and mb denote, respectively, the peaks of the inlet flow rate at the two entrances and we take ma = 1.0ped/m/s and mb = 0.5ped/m/s for the calculation. Then, flow vector Fa = (q a , 0) on Γao and Fb = (0, q b ) on Γbo . The free-flow walking speed is taken as 1.034m/s and two model parameters, α, β in (2.2) and (2.3), are taken as 0.075, 0.019, respectively (see [38]). The validation for the convergence of numerical solutions generally involves comparing them with the exact solutions of the problem. However, the exact solutions of the system model (2.12) are unknown. We thus have to use other methods (e.g., by refining the mesh) to convince ourselves about the convergence of the numerical solutions. The numerical experiments on four different sizes of unstructured meshes are tested: Mesh 1 with h = 1.6m; Mesh 2 with h = 0.8m; Mesh 3 with h = 0.4m; and Mesh 4 with h = 0.2m.

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1551

Comparing the numerical solutions obtained using the numerical methods mentioned in Ref. [19], i.e., the 1st order accuracy upwind scheme for the conservation law equation coupled with the 1st order FSM on orthogonal grids, indicates that the two numerical methods achieve almost the same results on different types of grids (see Figures 5 and 6). Here, four different partition sizes that correspond to the aforementioned mesh sizes are applied to generate orthogonal grids. The density curves of Group a along the center line, x = 50m, of the computational domain with two different methods are shown in Figure 5. Figure 6 depicts the density curves of Group b along the center line, y = 25m, of the computational domain with two different methods. From these figures, it can be seen that the densities as functions of time in decreasing discrete sizes converge toward a unique solution. In addition, it is also believed that the numerical solutions are stable and can be found with the finest mesh.

1

1

100*50 200*100 300*150 400*200

0.9 0.8

Mesh 1 Mesh 2 Mesh 3 Mesh 4

0.8

0.7 0.6

ρ

ρ

0.6 0.5 0.4

0.4

0.3 0.2

0.2

0.1 0

10

20

30

0

40

0

10

20

Fig.5

30

40

50

Y

Y

(Color online) Density plot of Group a at x = 50m and at t = 90s

1

1

100*50 200*100 300*150 400*200

0.9 0.8

Mesh 1 Mesh 2 Mesh 3 Mesh 4

0.8

0.7 0.6

ρ

ρ

0.6 0.5 0.4

0.4

0.3 0.2

0.2

0.1 0

20

40

60

X

(a) upwind Fig.6

80

0

0

20

40

60

80

100

X

(b) MUSCL

(Color online) Density plot of Group b at y = 25m and t = 120s

The densities of the bi-directional pedestrian flows in the railway platform are plotted in Figure 7 at different times. This figure clearly shows that the walking area is not fully occupied based on the optimal strategy for making a path-choice decision. It can also be seen that the flows spread over the center of the groups and form two strips that are composed of pedestrians

1552

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

walking in the same direction; then, each pedestrian in his or her group attempts to move toward the center of the group. The two streams of pedestrians walking in crossing directions to their different destinations always tend to separate. The dynamic continuum model of bidirectional pedestrian flows can reproduce self-organized phenomena [8, 12, 20–22], as dynamic lane formation can be also observed in such flows. 0.2 0.4 0.6 0.8

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

4

40

20

20

20

0

1

0

20

0.2 0.4 0.6 0.8

40

X

60

80

0 100

0

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

4

3

3.25 3.5 3.75

4

20

0

0.2 0.4 0.6 0.8

20

40

X

60

80

0 100

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

80

100

80

100

80

100

4

Y

Y

40

20

20

0

20

40

X

60

80

100

0

0

20

(c) t = 90s 0.2 0.4 0.6 0.8

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

40

X

60

(d) t = 90s 4

0.2 0.4 0.6 0.8

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

4

Y

40

Y

40

20

20

0

20

40

X

60

80

100

0

0

20

(e) t = 120s 0.2 0.4 0.6 0.8

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

40

X

60

(f) t = 120s 0.2 0.4 0.6 0.8

4

1

1.25 1.5 1.75

2

2.25 2.5 2.75

3

3.25 3.5 3.75

4

Y

40

Y

40

20

20

0

2.25 2.5 2.75

(b) t = 60s

40

0

2

40

(a) t = 60s

0

1.25 1.5 1.75

Y

40

Y

40

0.2 0.4 0.6 0.8

0

20

40

X

60

(g) t = 180s Fig.7

80

100

0

0

20

40

X

60

(h) t = 180s

(Color online) Density plot (left: Group a; right: Group b)

Figure 8 depicts the flow vector (2.7) of each group at different times, which helps us to observe the distribution and evolution of the flows. Strip formation in the two intersecting pedestrian streams is displayed. At the end of the modeling period, i.e., t = 240s, almost all pedestrians get away from the railway platform.

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1553

At the beginning of the modeling period when the two groups do not interact, the pedestrians move forward toward their goal, with almost none of them changing direction. As time marches on, the two groups arrive in the common area and are then interactive and passed by each other. As a result, the primary walking directions of the two groups vary. In addition, high-density values can be seen near the two exits, as the crowded conditions force the pedestrians’ speed to decrease as they line up to walk through the relatively narrow exits. When pedestrians want to pass through a bottleneck (e.g., coming into being in the common area and at the exits), interaction and negotiation become important to crossing flows and crowd flows [22].

40

40

30

30

Y

50

Y

50

20

20

10

10

0

0

20

40

X

60

80

100

0

0

20

(a) t = 30s

30

Y

40

30

Y

50

40

20

20

10

10

0

20

40

X

60

80

100

0

0

20

(c) t = 90s

40

30

30

Y

40

Y

50

20

20

10

10

0

20

40

X

60

80

100

0

0

20

(e) t = 150s

30

Y

40

30

Y

50

40

20

20

10

10

0

20

40

X

60

(g) t = 210s Fig.8

80

100

40

X

60

80

100

40

X

60

80

100

80

100

(f) t = 180s

50

0

60

(d) t = 120s

50

0

X

(b) t = 60s

50

0

40

80

100

0

0

20

40

X

60

(h) t = 240s

(Color online) Flow vector plot (red: Group a; green: Group b)

1554

5

ACTA MATHEMATICA SCIENTIA

Vol.29 Ser.B

Conclusion

In this work, the reactive dynamic user equilibrium model reported by Hughes [16] and Huang et al. [19] is extended to simulate two groups of pedestrians walking in crossing directions in a continuous walking facility. In this model, strategic path choice is the decision each pedestrian in a group makes to minimize the travel cost to his or her destination in a reactive manner based on the instantaneous available information. By numerically solving the model, a number of results, e.g., the formation of strips composed of pedestrians walking to the same destination, are observed. Lane or strip formation embodies pedestrians’ cooperative strategy when moving in different directions [22]. A relatively high degree of density exists at the bottlenecks, i.e., the two exits, which causes traffic congestion. The numerical results obtained with this continuum model can serve as a fundamental basis for congestion control and traffic management. References [1] Zacharias J. The nonmotorized core of Tianjin. Int J Sustainable Transport, 2007, 1(4): 231–248 [2] Cervero R, Sarmiento O L, Jacoby E, et al. Influences of built environments on walking and cycling: Lessons from Bogota. Int J Sustainable Transport, 2009, 3(4): 203–226 [3] Hughes R L. The flow of human crowds. Annu Rev Fluid Mech, 2003, 35(10): 169–182 [4] Helbing D, Buzna L, Johansson A, et al. Self-organized pedestrian crowd dynamics: Experiments, simulations, and design solutions. Transport Sci, 2005, 39(1): 1–24 [5] Helbing D, Johansson A, Al-Abideen H Z. The dynamics of crowd disasters: An empirical study. Phys Rev E, 2007, 75: 046109 [6] Muramatsu M, Irie T, Nagatani T. Jamming transition in pedestrian counter flow. Physica A, 1999, 267: 487–498 [7] Blue V J, Adler J L. Cellular automata microsimulation of bi-directional pedestrian flows. Transport Res Rec, 2000, 1678: 135–141 [8] Blue V J, Adler J L. Cellular automata microsimulation for modeling bi-directional pedestrian walkways. Transport Res B, 2001, 35: 293–312 [9] Asano M, Sumalee A, Kuwahara M, et al. Dynamic cell-transmission-based pedestrian model with multidirectional flows and strategic route choices. Transport Res Rec, 2007, 2039: 42–49 [10] Helbing D, Moln´ ar P. Social force model for pedestrian dynamics. Phys Rev E, 1995, 51: 4282–4286 [11] Helbing D, Farkas I, Vicsek T. Simulating dynamical features of escape panic. Nature, 2000, 407(6803): 487–490 [12] Helbing D, Moln´ ar P, Farkas I J, et al. Self-organizing pedestrian movement. Environ Planning B, 2001, 28(3): 361–383 [13] Hoogendoorn S P, Bovy P H L. Pedestrian route-choice and activity scheduling theory and models. Transport Res B, 2004, 38(2): 169–190 [14] Hoogendoorn S P, Bovy P H L. Dynamic user-optimal assignment in continuous time and space. Transport Res B, 2004, 38(7): 571–592 [15] Hughes R L. The flow of large crowds of pedestrians. Math Comput Simul, 2000, 53(4): 367–370 [16] Hughes R L. A continuum theory for the flow of pedestrians. Transport Res B, 2002, 36(6): 507–535 [17] Burstedde C K, Klauck K, Schadschneider A, et al. Simulation of pedestrian dynamics using a twodimensional cellular automaton. Physica A, 2001, 295: 507–535 [18] Yu W J, Johansson A. Modeling crowd turbulence by many-particle simulations. Phys Rev E, 2007, 76: 46105 [19] Huang L, Wong S C, Zhang M P, et al. Revisiting Hughes dynamic continuum model for pedestrian flow and the development of an efficient solution algorithm. Transport Res B, 2009, 43(1): 127–141 [20] Hoogendoorn S P, Daamen W. Self-organization in pedestrian flow//Hoogendoorn S P, Luding S, Bovy P H L. Traffic and Granular Flow’03. Berlin: Springer, 2005: 373–382

No.6

Y.Q. Jiang et al: A REACTIVE DYNAMIC CONTINUUM USER EQUILIBRIUM MODEL

1555

[21] Hoogendoorn S P, Daamen W. Pedestrian behavior at bottlenecks. Transport Sci, 2005, 39(2): 147–159 [22] Daamen W, Hoogendoorn S P. Experimental research of pedestrian walking behavior. Transport Res Rec, 2003, 1828: 20–30 [23] Helbing D. Traffic and related self-driven many-particle systems. Rev Mod Phys, 2001, 73(4): 1067–1141 [24] Schadschneider A, Klingsch W, Klupfel H, et al. Complex Dynamics of Traffic Management. New York: Springer, 2001 [25] Kirchner A, Schadschneider A. Simulation of evacuation processes using a bionics-inspired cellular automaton model for pedestrian dynamics. Physica A, 2002, 312: 260–276 [26] Weng W G, Chen T, Yuan H Y. Cellular automaton simulation of pedestrian counter flow with different walk velocities. Phys Rev E, 2006, 74: 036102 [27] Parisi D R, Dorso C O. Microscopic Dynamics of Pedestrian Evacuation. Physica A, 2005, 354: 606–618 [28] Lakoba T I, Kaup D J, Finkelstein N M. Modifications of the Helbing-Moln´ ar-Farkas-Vicsek social-force model for pedestrian evolution. Simulation, 2005, 81(5): 339–352 [29] Muramatsu M, Nagatani T. Jamming transition of pedestrian traffic at a crossing with open boundaries. Physica A, 2000, 286: 377–390 [30] Helbing D, Isobe M, Nagatani T, et al. Lattice gas simulation of experimentally studied evacuation dynamics. Phys Rev E, 2003, 67: 067101 [31] Guo R Y, Huang H J. A mobile lattice gas model for simulating pedestrian evacuation. Physica A, 2008, 387: 580–586 [32] Xia Y H, Wong S C, Zhang M P, et al. An efficient discontinuous Galerkin method on triangular meshes for a pedestrian flow model. Int J Numer Meth Engng, 2008, 76(3): 337–350 [33] Barth T J. Aspects of unstructured grids and finite volume solvers for the Euler and Navier-Stokes equations. AGARD Report 787, 1999, 6.1–6.61 [34] Wang J W, Liu R X. A comparative study of finite volume methods on unstructured meshes for simulation of 2D shallow water wave problems. Math Comput Simul, 2000, 53(3): 171–184 [35] Weller H, Weller H G. A high-order arbitrarily unstructured finite volume model of the global atmosphere: Tests solving the shallow water equations. Int J Numer Meth Fluids, 2008, 56(8): 1589–1596 [36] Sethian J A, Vladimirsky A. Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured mesh. Proc Natl Acad Sci, 2000, 97(11): 5699–5703 [37] Gremaud P A, Kuster C M. Computational study of fast methods for the Eikonal equations. SIAM J Sci Comput, 2006, 27(6): 1803–1816 [38] Wong S C, Leung W L, Chan S H, et al. Bi-directional pedestrian stream model with oblique intersecting angle. ASCE J Transport Engng, 2008, in press [39] Wong S C. Multi-commodity traffic assignment by continuum approximation of network flow with variable demand. Transport Res B, 1998, 32(8): 567–581 [40] Wong S C, Lee C K, Tong C O. Finite element solution for the continuum traffic equilibrium problems. Int J Numer Meth Engng, 1998, 43(7): 1253–1273 [41] Ho H W, Wong S C, Loo B P Y. Combined distribution and assignment model for a continuum traffic equilibrium problem with multiple user classes. Transport Res B, 2006, 40(8): 633–650 [42] Ho H W, Wong S C. Housing allocation problem in a continuum transportation system. Transportmetrica, 2007, 3(1): 21–39 [43] Ho H W, Wong S C. Existence and uniqueness of a solution for the multi-class user equilibrium problem in a continuum transportation system. Transportmetrica, 2007, 3(2): 107–117 [44] Tong C O, Wong S C. A predictive dynamic traffic assignment model in congested capacity-constrained road networks. Transport Res B, 2000, 34(8): 625–644 [45] Toro E F. Shock-Capturing Methods for Free-Surface Shallow Water Flows. John Wiley & Sons, 2001 [46] Shepard D. A two-dimensional interpolation function for irregularly spaced data//Proc 23rd Nat Conf, ACM. New York, 1968: 517–524 [47] Kermani M J, Gerber A G, Stockie J M. Thermodynamically based moisture prediction using roe’s scheme//The 4th Conference of Iranian AeroSpace Society. Amir Kabir Univ of Tech, Tehran, 2003 [48] Zoppou C, Roberts S. Explicit schemes for dam-break simulations. J Hydraul Engng, 2003, 129(1): 11–34 [49] Shu C W, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys, 1998, 77(2): 439–471