Abstract
The macroscopic models for solving the pedestrian flow problem can be generally classified into two categories as follows: first-order continuum models and high-order continuum models. In first-order continuum models, the density satisfies the mass conservation law, the speed is defined by a flow-density relationship, and the desired directional motion of pedestrians is determined by an Eikonal-type equation. In contrast, in high-order models, the velocity is governed by the momentum conservation law. In this study, we summarize existing first-order and high-order models and rewrite them in the form of unified scalar or system hyperbolic conservation laws. Next, we apply high-order discontinuous Galerkin methods with a positivity-preserving limiter on unstructured triangular meshes to solve the conservation law and a second-order fast-sweeping scheme to solve the Eikonal equations. Our method can efficiently model real-life complex computational regions and avoid nonphysical solutions and simulation blow-ups. Finally, numerical examples are presented to demonstrate the accuracy and effectiveness of the proposed solution algorithm. The numerical results validate the reliability of the proposed numerical method and highlight the advantages of triangular meshes.
1. Introduction
Due to the importance of the security and safety management of pedestrian traffic, the pedestrian and crowd dynamic problem has received considerable attention in recent decades. Models for simulating pedestrian flows can be classified into microscopic and macroscopic models. Commonly used microscopic models such as social force models, discrete element models, and cellular automaton-based models simulate pedestrian movement by quantitatively describing the physical and psychological interactions between individuals [1–3]. These models enable the simulation of pedestrian dynamics and the reproduction of empirical observations, such as the decision-making process [4], uni- and bidirectional fundamental diagrams [5], and pedestrian movements at bottlenecks [6] at low densities. However, at high crowd density levels, microscopic models are not as computationally efficient as macroscopic models [7]. In addition, unlike microscopic models, dynamic macroscopic models can simulate the evolution of collective pedestrian behaviors, thus facilitating the analyses of key variables and parameters of crowd dynamics.
Macroscopic models can enhance the understanding of pedestrian flow dynamics by efficiently simulating large-scale and high-density crowds. The variables (e.g., flow intensity, density, and velocity) in a macroscopic model can be represented by smooth mathematical functions, and their relationship can be described by a set of partial differential equations. Among these models, the Lighthill–Whitham–Richards (LWR) model [8, 9] is one of the simplest continuum pedestrian models. Although this model was originally independently proposed by Lighthill and Whitham [8] and Richards [9] (in 1950) to solve traffic models, it is widely used to study pedestrian dynamics due to the similarity between vehicular and pedestrian flows. In this model, a mass conservation law is used to describe the pedestrian flow, assuming an equilibrium state of the flow-density relationship, similar to that reported by Greenshields [10]. The traditional LWR model has several limitations, and thus, many extended LWR models have been designed to solve problems such as the multiclass problem [11]. Specifically, the traditional LWR model is defined in a one-dimensional (1D) space, in which the walk direction is unidirectional and known. The route-choice strategy can be identified using the flow-density relationship. In a two-dimensional (2D) pedestrian flow problem, although the speed is known, the walk direction is no longer a single entity, and many (even infinite) feasible walking directions may exist. Consequently, a mathematical formulation must be obtained to describe the route-choice strategy. To describe the route-choice strategy in a 2D pedestrian flow problem, Hughes [12] proposed a systematic modeling framework based on the continuum theory and used an Eikonal equation to implement the route-choice strategy. Subsequently, Hoogendoorn and Bovy [13, 14] developed a dynamic user-optimal pedestrian flow model to study user-optimal dynamic traffic assignment problems in continuous time and space, in which pedestrians are assumed to have perfect information regarding the modeled domain that can be used to choose the route that minimizes the actual walk cost. In this model, a time-dependent Hamilton–Jacobi equation is used to describe the route-choice strategy. Huang et al. [15] reformulated Hughes’ model and discovered that the pedestrian route-choice strategy in Hughes’ model satisfies the reactive dynamic user equilibrium (RDUE) principle, in which a pedestrian chooses a route to minimize their instantaneous travel cost to the destination. Xiong et al. [16] extended Huang’s reactive DUE model to study bidirectional pedestrian flow problems. Later, Du et al. [17] proposed a predictive dynamic user equilibrium (PDUE) model, in which a pedestrian chooses a route to minimize their actual travel cost to the destination.
In the abovementioned 1D or 2D LWR models or extended LWR models, the speed is defined by a flow-density relationship. Such models are typically known as first-order or equilibrium models, and numerical simulations have demonstrated their usefulness in planning and designing walking facilities. Although equilibrium models can sensibly predict pedestrian dynamics, they cannot simulate the stop-and-go waves and forward propagation of disturbances in heavy traffic [18]. Thus, nonequilibrium models (second- or high-order models) have been designed. The Payne–Whitham (PW) model [19, 20] is the earliest known nonequilibrium model. The speed in the PW model satisfies the momentum equation rather than a flow-density relationship, as in the case of LWR models. The PW model has been used to study the wrong-way travel problem [21, 22], and many PW-like models have been developed, using different traffic sound speeds to solve different problems. In general, PW and PW-like models have the same structure, incorporating a mass conservation law and a momentum equation. Such models are strictly hyperbolic when the characteristic speed exceeds the vehicle speed, which may lead to gas-like behavior [23]. Aw and Rascle [24] and Zhang [25] developed a nonequilibrium model devoid of this behavior. To describe the property of crowd pressure and identify the force chains, Liang et al. [26] proposed an extended PW (PWP) model that explicitly considers the crowd pressure and effects of pushing forces and high densities to describe pedestrian movement.
The abovementioned first- and higher-order models can be rewritten into a unified form that consists of a scalar or system conservation law coupled with an Eikonal equation (or Hamilton–Jacobi equation). In studies of numerical solution of these pedestrian flow models, finite difference (FD) [11, 14, 26–29] and finite volume (FV) methods [30] have been commonly used. These methods can be directly used to solve many pedestrian flow problems on structured grid points, rectangular meshes, or some simple unstructured meshes. However, FD methods cannot address problems involving real and complex computation regions, for which unstructured grids must be used. Moreover, FV methods are difficult to implement because of the challenges associated with reconstructions on complex meshes. Direct numerical simulations may lead to significant oscillations, nonphysical solutions, or even blow-ups when solving shock problems.
Discontinuous Galerkin (DG) methods are a class of finite element methods that have been widely used to solve hyperbolic equations. These methods use completely discontinuous basis functions and are thus more flexible than traditional continuous finite element methods. In addition, DG methods are characterized by a high-order accuracy, high parallel efficiency, flexibility for -adaptivity, and use of arbitrary geometries and meshes. The first DG method, introduced in 1973 by Reed and Hill [31], was used to solve a time-independent linear hyperbolic equation. Subsequently, Cockburn et al. [32–35] established a framework to solve nonlinear time-dependent problems by combining the explicit, nonlinearly stable high-order Runge–Kutta time discretization scheme [36] with the spatial DG discretization. In this method, exact or approximate Riemann solvers are used as interface fluxes, and total variation-bounded nonlinear limiters [37] are used to achieve nonoscillatory properties for strong shocks. There are some extended DG methods, such as local DG [38], direct DG [39], hybrid DG [40], and oscillation-free DG [41], for partial differential equations (PDEs) that contain higher-order spatial derivatives or to reduce the oscillation. This method is widely applied in the domains of aeroacoustics, gas dynamics, granular flows, magnetohydrodynamics, modeling of shallow water, oceanography, oil recovery simulations, turbulent flows, viscoelastic flows, and weather forecasting. In this study, we implement the high-order DG method over unstructured triangular meshes to solve the conservation law and establish a systematic numerical algorithm for modeling pedestrian dynamics.
Physically, the pedestrian density must be nonnegative in a continuum model. However, in numerical experiments, this property does not hold for many numerical schemes. The failure to preserve a positive density leads to ill-posedness of the models and can easily cause blow-ups, especially when high-order numerical schemes are applied to low-density regions. In many applications, the negative values are simply truncated to zero or small positive values. Although this approach may be effective in certain problems, brute truncation will eventually result in blow-ups because the total mass increases every time a negative density is set as zero. To ensure stable computing, schemes that are demonstrably positivity-preserving while also preserving the conservation property must be designed. In reference [42], the authors first constructed an arbitrary high-order accurate maximum-principle-satisfying DG scheme for one-dimensional and multidimensional scalar conservation laws on rectangular meshes. Subsequently, the authors extended this method to a high-order accurate positivity-preserving scheme for Euler systems on rectangular and triangular meshes [43, 44].
To fill the research gap of robust numerical algorithms on triangular meshes for the reviewed pedestrian models, this study applies a genuinely high-order DG scheme with a positivity-preserving limiter on unstructured triangular meshes to solve the 2D conservation law, enabling the realization of robust simulations in complex computational domains. In addition, we use the second-order fast-sweeping method to solve the Eikonal equation that describes the route choice in first-order models or equilibrium speed in higher-order models. In numerical examples, we first provide a numerical example that has an exact solution to demonstrate the accuracy of the high-order scheme. Then, we test different cases of the PWP model. First, to validate the correctness of the schemes, we test an evacuation case involving a rectangular obstacle. The results are consistent with those obtained from the FD weighted essentially nonoscillatory method [26]. The simulations based on the two methods yield consistent numerical results, which confirm that the DG method can be applied to nonsmooth and unstructured meshes without the loss of accuracy or conservation properties. Second, an evacuation case with circular obstacles is designed to validate the higher-order model and highlight the advantages of the proposed approach. The averaged local flow-density relationship derived from the simulation produces a second-peak phenomenon, consistent with the findings of empirical studies [18, 45]. Moreover, the occurrence of stop-and-go waves is observed, which verifies the applicability of the proposed model to unstable crowd conditions. In addition to the density distribution, the aggregated pushing pressure calculated using the proposed model can serve as another risk-level indicator to strategically prevent crowd disasters. Thus, the high-order DG scheme with a positivity-preserving limiter is demonstrated to be accurate and effective.
The remainder of the paper is organized as follows. Section 2 reviews different continuum pedestrian dynamic models and describes the method of rewiring them in a unified form. Section 3 describes the numerical algorithm, including the DG scheme with a positivity-preserving limiter for hyperbolic equations, fast-sweeping method for Eikonal equations, and total variation diminishing (TVD) Runge–Kutta temporal discretization. Section 4 presents the numerical examples designed to test the properties and effectiveness of the proposed methods. Section 5 presents the concluding remarks.
2. Pedestrian Dynamic Model
When studying the pedestrian and crowd dynamics through the continuum modeling approach, the pedestrian dynamic model can be formulated as a set of partial differential equations (PDEs) or PDEs with relaxation. In this section, we review several commonly used continuum pedestrian models, including first-order and high-order models, and rewrite them into a unified form. The considered first-order models are the traditional LWR model and 2D extended LWR models (such as RDUE/PDUE). Among high-order models, different second-order models are introduced.
2.1. First-Order Models
We first introduce the first-order models, in which pedestrians always choose routes that minimize their cost. The LWR model is one of the simplest continuum pedestrian dynamic models, in which the pedestrian flow is considered to be in a 1D space. The traffic state variables such as the flow, speed, and density satisfy the following mass conservation law:where is the pedestrian density and is the travel demand. is the preferred speed, which is a given nonnegative and nonincreasing function with respect to and defined in ; is the maximum density, corresponding to a traffic jam. Thus, the flow flux is defined as follows:In the LWR model, because is a given function of , the flow flux satisfies several fundamental diagrams, such as the Greenshields model.
The traditional LWR model is defined in a 1D space, in which the pedestrian walks forward to the destination and the speed (norm of the velocity vector) function is given. Thus, the velocity is known. Although this model is useful for analytical problems, it cannot be applied to real-world problems because pedestrians typically walk in 2D scenarios. Therefore, we review the 2D-extended LWR models coupled with different route strategy assumptions. In 2D space, the following conservation law is satisfied:where is the 2D-preferred velocity vector. Although the speed intensity can be computed through a given speed-density relationship , as in the 1D case, pedestrians can move freely in any direction. Consequently, additional route-choice strategies are required to describe the travel direction of . Commonly used route strategy models include RDUE and PDUE models.
In RDUE models, pedestrians always choose routes that minimize their instantaneous walk cost and change their movement directions in a reactive manner. The route-choice strategy satisfieswhere // indicates that two vectors are parallel. Here, is defined as the instantaneous cost potential, which can be computed using the following Eikonal equation:where is the local cost function.
In the PDUE model, pedestrians are assumed to have perfect information regarding the modeling domain that can be used to predict future traffic information. Thus, the users can choose routes that minimize their actual walk cost. The route-choice strategy satisfies the form presented in equation (4). However, the total travel cost must be redefined as the actual travel cost and computed using a time-dependent Hamilton–Jacobi equation as follows:
According to the route-choice strategy presented in equation (4), the density in both the RDUE and PDUE models is governed by the following 2D conservation law:where is the flow flux, which is defined as follows:
2.2. High-Order Models
In first-order models, the magnitude of the movement speed is a given function of the density (corresponding to the fundamental diagram), and the movement direction is computed using a suitable route-choice strategy. The existing fundamental diagrams are typically based on empirical observations and indicate the average behavior of moving pedestrians. However, evident individual differences and unstable movement patterns have been observed from empirical data. To describe these behavioral characteristics of pedestrian movement, the second-order PW model with relaxation has been designed. Similar to LWR models, the first equation in the PW model still expresses the mass conservation, and the density, flux, and travel demand satisfy the following conservation law:where are the velocities in the and directions, respectively. In contrast to the LWR models, the velocities are no longer computed by the given speed functions and route-choice algorithms and satisfy the following momentum equations:where is the relaxation time; are the equilibrium speed in the and directions, respectively; and is the traffic sound speed. Because is a negative parameter, the traditional PW model may produce negative speeds (i.e., wrong-way travel). To address this problem, many nonequilibrium PW models have been designed, in which the constant traffic sound speed is replaced by a function that depends on the density . The new momentum equations are defined as follows:
The mathematical structure of these improved models is similar to that of the PW model, and thus, such models are named PW-like models.
To eliminate the gas-like behavior in PW-like models, Aw and Rascle [24] and Zhang [25] developed a nonequilibrium model shown to be devoid of the gas-like behavior. In this model, the momentum equation in PW-like models is replaced with the following equations:where are the pressures in the and directions, respectively, which are increasing functions of the density.
The abovementioned higher-order models focus on the simulation of normal pedestrian movement in low-density scenarios. To address large-scale and dense problems such as crowd disasters, Liang et al. [26] proposed a higher-order continuum model that considers panic effects and explicitly takes into account the effect of the aggregated pushing force. The model formulation iswith suitable initial boundary conditions. Here, is the average mass of a pedestrian, and is the compression force, which is determined by the pedestrian density as follows:where is a given function that indicates the relationship between the average compression force and density. is the pushing force, which is generated by panic in a dense flow and satisfies the following equation:where is a given function that describes the relationship between the pushing capacity and density; is a given parameter that describes the level of panic in the crowd; and represents the relaxation factor defined as follows:where represents the critical density when the pushing capacity is higher than 0, and represents the maximum density.
2.3. Unified Form of Different Models
In the previous sections, we reviewed several classical pedestrian dynamic models including the first-order models (equilibrium models) and high-order models (nonequilibrium models). All these models can be rewritten in the following unified form. Note that we present the 2D cases of the models, except for the traditional LWR model as follows:
Table 1 presents the variables , , , and for the different models.
The equilibrium velocity vector is computed using equations (4) and (5) (or equation (6)).
Remark. Many other continuum pedestrian models are available, but they have not been introduced in this section due to limited space. Most of these models can also be written as the conservation law, with different definitions of the conservation variables, pressure , source terms, and equilibrium speeds.
3. Solution Algorithms
Equation (17) represents a 2D system of hyperbolic conservation laws, and we define . In this model, the source term is computed implicitly. At each fixed time, equation (5) must be solved to determine the equilibrium speed . Moreover, in the PWP model, equation (15) must be solved to obtain the gradient of and compute the source term . Equations (5) and (15) can be written as the following Eikonal-type equation with suitable boundary conditions:
Several numerical algorithms have been designed to solve the conservation law system. The DG method is one of the most popular algorithms and can be used to study problems with complex geometries by using unstructured meshes. In this study, we use the high-order DG method with a positivity-preserving limiter on triangular meshes to solve equation (17). At each time level of the DG method, we use the second-order fast-sweeping method on triangular meshes to solve the Eikonal equation, equation (18).
The formulation of the fast-sweeping method is presented in Section 3.1. The high-order DG method and positivity-preserving limiter are discussed in Sections 3.2 and 3.3, respectively.
3.1. Fast-Sweeping Method for Solving Eikonal Equations
We apply the fast-sweeping method to solve equation (18) on triangular meshes. Details of first- and second-order fast-sweeping methods on triangular meshes can be found in References [46, 47], respectively. In this study, the second-order fast-sweeping method is used.
We discretize the computational domain into unstructured triangular cells . As shown in Figure 1, for a given triangle , the three angles are and the lengths of the three edges are . The fast-sweeping method is an iterative method. Assuming that and at nodes and and their related derivatives are known, we aim to update at node .
Before introducing the local solver for updating , we examine the directional derivatives along edges and at node . The second-order approximations are defined as follows:
The following definitions are considered:
According to the connection between the gradient and directional derivatives,where , and superscript represents transposition. According to equation (21), the gradient at node can be approximated as follows:
If this equation is combined with equation (18), the value of at node satisfieswhere is the source term at location at time .
As presented in Reference [46], the computed solution from equation (23) can be accepted only if the following upwind criteria are satisfied:
Otherwise, the gradient satisfies
Thus, by solving equation (23) and combining it with equation (25), the following expression can be obtained:
This form pertains to the Godunov numerical Hamiltonian for Eikonal equations. According to equation (26),where .
Next, we introduce the sweeping process, which is a Gauss–Seidel iteration process. The sweeping order is a key parameter in the fast-sweeping method. The systematic order must efficiently cover all directions. The order of a rectangular mesh can be derived in a straightforward manner, and the following four sweeping directions are considered:where and are indices in the and directions, respectively. Notably, these natural orderings cannot be directly applied to unstructured meshes. A new systematic ordering must be identified to cover all directions. Following Reference [46], we choose at least three noncollinear reference points and then sort all nodes according to the distance to each reference point. In this manner, we can sweep all nodes based on their ascending and descending orders.
The process of the fast-sweeping method on a triangular mesh can be summarized as follows.(1)Initialization(a)Choose the reference points , and sort all nodes according to the distance to each reference point. Define and as the arrays of all nodes according to the distance to the -th reference point in decreasing and increasing orders, respectively.(b)Based on the boundary condition, assign the exact values and at the inflow boundary. Large values (e.g., ) are assigned as the initial estimates at all other grid points.(2)Gauss–Seidel iteration(c)Denote the vector of the initial solutions as and set .(d)Sweeping(i)Choose the first sweeping direction and set .(ii)Loop through all nodes in set and apply the local solver to update the value.(iii)Loop through all nodes in set and apply the local solver to update the value.(iv)If , proceed to Step 2(e); otherwise, set and return to Step ii.(e)If , stop; otherwise, return to Step 2(d). Here, is a convergence threshold and is the norm.
3.2. DG Methods
Equation (17) is solved using the DG method. The triangulation of the domain is , and is any triangular element. Equation (17) is multiplied by a test function and then integrated over each cell as follows:
Integrating by parts yieldswhere is the outward unit normal of the cell boundary . The line integral in this equation is typically discretized by a Gaussian quadrature with a sufficiently high order of accuracy as follows:where is the quadrature point and is the quadrature weight.
Next, we approximate solution by using a numerical solution within the finite element solution space , defined as follows:where represents the space consisting of polynomials of degree up to in element . In the finite element space , the functions can be totally discontinuous across the element interfaces. We replace function in equation (30) with . Moreover, we select test functions from the space . Function in equation (30) is replaced with . As and are allowed to have discontinuities across the cell interface, the flux term and test function on the cell boundary are not well defined. In the DG methods, on the cell boundary must be replaced with a numerical flux. Although many types of numerical fluxes are available, the Lax–Friedrichs flux is used in this study as follows:where and are the values of on the cell boundary obtained from the inside and outside of triangle , respectively, and . The test function on the cell boundary is set as , derived from the inside of . The DG method for solving equation (20) can be summarized as follows: find the unique solution such that . For any triangular element in the triangulation ,
Similar to that in equation (31), this cell boundary integral can be approximated using the Gaussian quadrature rule. Similarly, the term is computed by the following numerical quadrature:where is the quadrature point within and is the quadrature weight.
The basis functions for the space are for , where . Thus, the numerical solution can be represented as follows:Substituting this form into equation (34) and taking yieldwhere , and is a matrix with the entry , defined as follows:and is a vector with the -th entry, defined as follows:
Equation (37) can be rewritten as the following ordinary differential equation (ODE):
To discretize the time to solve this ODE system, we use the third-order TVD Runge–Kutta method [36, 48], which can maintain the stability of the spatial discretization as follows:
3.3. Positivity-Preserving Limiter
This section describes the positivity-preserving schemes for triangular meshes developed in Reference [42]. As TVD time discretizations are convex combinations of the Euler forward operators, we need to consider only the first-order Euler forward time discretization. In this study, we use the global Lax–Friedrichs scheme, but the results can be easily extended to other monotonic schemes.
For the test function in equation (34), we obtain a scheme that is satisfied by the cell averages of the DG method as follows:where is the cell average of numerical solutions over triangle at the time level and and are the values of on the -th boundary obtained from the inside and outside of the triangle , respectively. For the first component of , i.e., the density ,where and are the first components of and , respectively. Let be the -th degree polynomial on triangle . Then, the line integral must be approximated by the least Gauss quadrature point. This scheme can be rewritten as follows:where is the Gauss quadrature weight on the interval and . and represent the values of at the -th Gauss quadrature point on the -th edge, taken from inside and outside the triangle , respectively; and is the length of the -th edge.
Next, we decompose the cell average as a convex combination of the point values of the DG polynomial . Following Reference [44], we denote the set of quadrature points as follows:where are the position vectors of the three vertices of triangle ; and and are the Gauss–Lobatto and Gauss quadrature points on the interval , respectively, where is the smallest integer such that . If , the set of quadrature points is as shown in Figure 2. The cell average of can be decomposed as follows:where is the value of the polynomial at the quadrature points on edge , is the value of the polynomial at the quadrature points in the interior of , and the number of interior points is . The coefficient is . Similarly, the integral of the source term can be decomposed as follows:where is the value of at the quadrature points on edge and is the value of at the quadrature points in the interior of .
Next, we consider equation (42), which is satisfied by the cell averages of the DG method. If is nonnegative for any in the Courant–Friedrichs–Lewy condition with and , then is also nonnegative. Here, is the length of the -th edge of triangle ; is the quadrature weight of the -point Gauss–Lobatto rule on for the first quadrature point; is the value of at point ; and is a function of and and , which is positive for positive density values. Its explicit form depends on the specific source term q [49].
At the time level , the original DG polynomial does not satisfy the sufficient condition in this theorem in certain cases. Thus, we need to modify the polynomial such that it is nonnegative for all . For any triangle, we assume that and define as the DG polynomial for in triangle . Next, we introduce the linear scaling limiter to derive the modified polynomial as follows:with
According to equations (48) and (49), if the density is always nonnegative, the parameter always takes a value of . In this case, and the limiter does not work. However, it can be easily derived that for all . Therefore, we use the modified polynomial to replace in equation (42). According to Theorem 1, the cell average at the time level is also nonnegative.
4. Numerical Examples
In this section, we first provide a numerical test with a simplified setup to the model which has an exact solution, to demonstrate the accuracy of the high-order positivity-preserving DG scheme. Then, the PWP model is considered as an example to evaluate the proposed numerical algorithms.
4.1. Example 1
We solve the following model equation:where and and are the source terms, chosen so that an explicit exact solution is available. In this example, we solve the above system on a domain of , and the initial data take . The exact solutions for this problem are when the source terms are .
The errors and orders of accuracy are shown in Table 2. We can clearly see that the desired second-order of accuracy is achieved with the error of ρ on the triangular meshes.
4.2. Numerical Example of the PWP Model
We consider two panic evacuation scenarios on a rectangular railway platform with obstacles of different shapes. Specifically, the platform is considered to have rectangular and circular obstacles in Examples 2.1 and 2.2, respectively.
In both examples, the railway platform is initially empty, with at . Over time, pedestrians enter the platform from the left boundary and move from the left to the right within the platform. The boundary conditions for both examples are set as follows. The origin is at , . The inflow at the origin iswhereand the unit of is . The destination is , where the cost potential is zero. The upper () and lower () boundaries are solid wall boundaries, and the projection of the speed vector on the normal direction of a boundary is set as zero. Table 3 specifies the other parameters and functions.
4.2.1. Example 2.1
A rectangular obstacle is located at (unit: m). The proposed methods are used to solve the PWP model on triangular meshes. Figure 3 illustrates a sample triangular mesh with cells. The characteristics of the numerical results obtained by the DG method are consistent with those of the findings of our previous study [26].
Figure 4 compares the density cut along the lines y = 40 m for the PWP model simulated by the DG scheme without and with the positivity-preserving limiter at (see subfigures 4(a) and 4(b)). The DG scheme without the positivity-preserving limiter produces negative solutions, whereas the positivity of the density is maintained when using the DG scheme with the positivity-preserving limiter. The solution will blow up if the PWP model is continually computed using the DG scheme without the positivity-preserving limiter.
(a)
(b)
Figure 5 shows the density curves in the 1D section corresponding to m (sub-Figure 5(a)) and (subfigure 5(b)) at for different mesh sizes, as obtained from the DG method. The numbers of elements in the triangular meshes are , and 9893. By gradually refining the mesh, the deviation between the density curves can be reduced. The convergence of both numerical solutions is thus validated.
(a)
(b)
Figure 6 plots the density distributions during the evacuation at different time points, determined using the second-order DG method on triangular meshes with N = 9893. Before , the pedestrian flow is characterized by calm behaviors; the pedestrians try to avoid the obstacle and move at the free-flow speed (see subfigure 6(a)). At , congestion is induced near the corner of the rectangle because the obstacle decreases the available space, leading to a maximum density of approximately . The passage is further blocked by the congested pedestrian flow (see subfigure 6(b)). After , panic begins to propagate in the crowd (see subfigures 6(c)–6(f)). In this condition, the density increases more slowly than in the normal condition as it is already extremely high. However, pedestrians continue to enter the flow, generating a high pushing pressure in the dense crowd, as shown in Figure 7. The pushing pressure propagates from the low-density region to the high-density region (close to the obstacle) and reaches a value of approximately (see subfigures 7(a) and 7(b)). The pressure increases rapidly in the dense crowd, and the density level remains high. In other words, the crowd pressure characteristics, which have been observed in many crowd disasters [49–51], are correctly rendered.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
The averaged local flow-density relationship obtained by the DG method on triangular meshes also reveals the second-peak phenomenon, consistent with the findings reported by Reference [26]. The averaged local flow and corresponding density are calculated in the simulation of a panic evacuation between by using the method presented in Reference [26]. In Figure 8, these entities correspond to the “actual flow” and are compared with the values derived from the fundamental diagram. A notable second peak occurs at approximately due to the panic sentiment and aggregated pushing pressure. Helbing et al. [18] first observed this phenomenon when they conducted a postdisaster analysis of a Hajj crowd disaster and related it to the forward and backward moving shock waves in the dense crowd. In our model, we generate the waves by considering the aggregated pushing pressure created by the panic, thereby reproducing this complex phenomenon. Therefore, the realistic nature of the results of the continuum model is validated.
4.2.2. Example 2.2
Example 2.2 corresponds to a more realistic computational domain than that in Example 2.1. Three circular obstacles are placed at with radii of (unit: m), respectively. This complex region cannot be easily divided by structured meshes, and thus, the traditional FD or FV methods cannot be applied. Therefore, we apply the DG methods over unstructured triangular meshes, as illustrated in Figure 9.
Figure 10 plots the density distribution during the evacuation at different time points. Before , the pedestrian flow is characterized by calm behaviors; the pedestrians try to avoid the obstacles and follow the shortest path to their destinations (see subfigure 10(a)). However, the three circular obstacles cause congestion upstream and downstream of the obstacles, with the maximum density being approximately (see subfigure 10(b)). This maximum density is smaller than that in example 2.1 (with rectangular obstacles), which suggests that the influence of circular boundaries on pedestrian movement is weaker. This finding is supported by existing simulation results [52, 53]. After , panic begins to propagate in the crowd (see subfigures 10(c)–10(f)). The density gradually increases to approximately , following a similar pattern to that in example 2.1. The pedestrians continue to enter the flow, generating a high pushing pressure in the dense crowd, as shown in Figure 11. The propagated pushing pressure increases to approximately , considerably lower than that in example 2.1 (see subfigures 11(a) and 11(b)).
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
4.2.3. Illustration of Instability Using the DG Method
At low sonic speeds ( in this study), the simulation results exhibit stop-and-go waves. At , three stop-and-go waves are observed near the largest circular obstacle, as illustrated in Figure 12(a). These waves are generated near the bottleneck and propagate backward to the low-density area until they vanish. Moreover, when the pedestrians panic, small and dense waves are generated and propagated continuously in the high-density regions, as illustrated in Figure 12(b). This phenomenon is similar to the turbulence observed by Fruin [52], who indicated that shockwaves can propagate through dense crowds. The DG scheme can successfully reproduce the instability of crowd dynamics induced at low sonic speeds, which validates the continuum model.
(a)
(b)
5. Conclusion
Several commonly used continuum pedestrian flow models are reviewed. According to the equilibrium state of speed, these models can be classified into first- and high-order models. The first-order models include the traditional LWR model and other extended LWR models, coupled with RDUE and PDUE route strategies. High-order models include the PW model, PW-like model, Aw–Zhang model, and PWP model. Both types of models can be rewritten in a unified form as a scalar or system conservation law with different conservation variables, flow fluxes, and source terms. To address problems in real and complex computational regions, a high-order DG scheme with a positivity-preserving limiter is applied to solve the conservation law, and a second-order fast-sweeping scheme is used to solve the Eikonal equations. The proposed solution algorithm can efficiently solve all types of pedestrian dynamic models considered in this study. We provide a numerical example that has the exact solution to demonstrate the accuracy of the high-order scheme and validate the proposed numerical method through numerical examples involving evacuation cases with a rectangular obstacle and three circular obstacles based on the PWP model. The numerical results (which indicate the occurrence of the second-peak phenomenon and stop-and-go waves) are consistent with those reported in our previous work. Thus, the proposed numerical algorithm is accurate and effective.
Notably, because of the limitations of the fast-sweeping method for solving the Eikonal equation on triangular meshes, the complete algorithm is no more than second-order accurate. Future research can be aimed at developing numerical methods with higher-order accuracy to solve the Eikonal equation on triangular meshes to ensure that the algorithm can reach third-order accuracy or higher with a higher-order DG scheme. In this manner, the solution efficiency and applicability to more complicated real-world cases can be enhanced.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Liangze Yang and Haoyang Liang were responsible for conceptualization, methodology, software, and writing of the original draft. S. C. Wong and Jie Du were responsible for conceptualization, methodology, review writing, and Editing.
Acknowledgments
The work described in this paper was supported by grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (project nos. 17201318, 17204919), Guangdong–Hong Kong–Macau Joint Laboratory Program of the 2020 Guangdong New Innovative Strategic Research Fund, Guangdong Science and Technology Department (project no. 2020B1212030009), and Tsinghua University Initiative Scientific Research Program. In addition, the second author was supported by a postgraduate scholarship from the University of Hong Kong, and the fourth author was supported by the Francis S Y Bong Professorship in Engineering.