Macroscopic pedestrian flow simulation using Smoothed Particle Hydrodynamics (SPH)

Macroscopic pedestrian models are theoretically simpler than microscopic models, and they can potentially be solved faster while producing reasonable predictions of crowd dynamics. Therefore, they can be very useful for applications such as large-scale simulation, real-time state estimation and crowd management. However, the numerical methods presently used to solve macroscopic pedestrian models, which are mostly grid-based, have some shortcomings that limit their applicability. More specifically, they usually include complex procedures for grid generation and remeshing, and they produce simulation results that may not be sufficiently accurate (for example, because of unclear boundaries between flow states). Smoothed Particle Hydrodynamics (SPH) constitutes an alternative numerical method that could potentially overcome these limitations. SPH is a meshfree method where a crowd is represented by a set of particles that possess material properties and move according to macroscopic laws. Relevant state variables at each particle are approximated using information about the material properties of the neighboring particles and a smoothing function. This paper puts forward for the first time a generic SPH framework for solving macroscopic pedestrian models; in addition, it demonstrates that an SPH-based simulation model can produce meaningful and accurate results by means of three case studies. The first case study shows that the proposed numerical method can approximate well the analytical solution of a simple macroscopic model applied to a queue-discharge scenario. The second case study demonstrates that the proposed numerical method can potentially reproduce density dispersion (a phenomenon observed in real crowds) more accurately than grid-based methods, due to its meshfree, Lagrangian, and particle-based nature. The third case study highlights the need to reformulate the acceleration equation of the basic macroscopic model in order to reproduce lane formation in bi-directional flows (also an observed phenomenon) using the proposed SPH framework, and this paper presents a solution to do so.


Introduction
Macroscopic pedestrian models (e.g., Hoogendoorn et al., 2014;Jiang et al., 2010) represent crowds as continuous media (using quantities such as flow, density and velocity) and describe pedestrian movements using the conservation laws of fluid dynamics.These models have the potential to reduce computation times in comparison to microscopic models (e.g., the social-force model by https://doi.org/10.1016/j.trc.2019.12.017Received 22 July 2019; Received in revised form 12 August 2019; Accepted 21 December 2019 ⁎ Corresponding author.E-mail address: y.yuan@tudelft.nl(Y.Yuan).
These applications require efficient (i.e., both fast and accurate) numerical methods to solve the model equations.However, the numerical methods generally used in macroscopic pedestrian simulation have shortcomings that limit their applicability when it comes to simulating realistic scenarios.These methods are generally grid-based (the grids are Eulerian, as in Hänseler et al. (2014), or Lagrangian, as in van Wageningen-Kessels et al. (2016)).Grid generation can be a rather difficult task, particularly in the case of walking facilities with complex geometries.Moreover, in Eulerian-grid methods, it is difficult to determine the precise location of traffic-state boundaries.In Lagrangian-grid methods, given the freedom of pedestrians to move in many directions, large deformation of mesh cells (representing groups of pedestrians) requires the use of remeshing techniques (e.g., van Wageningen-Kessels et al. (2016)), which adds complexities, and may lead to loss of information.
In computational fluid dynamics, meshfree methods constitute an alternative to grid-based methods.One of the most widely used meshfree methods is Smoothed Particle Hydrodynamics (SPH) (Monaghan, 2005), which was conceived for astrophysical research in the 1970s (Lucy, 1977;Gingold and Monaghan, 1977).In SPH, the state of the system is represented by a set of unconnected Lagrangian particles that possess material properties and move according to the conservation laws.Field variables at each particle are approximated using information about the material properties of the neighboring particles, and a selected smoothing function that has a bell shape.Note that the concept of discretizing continuum media into a finite number of particles was also adopted in another so-called "macro-particle" simulation approach, which was used by Leboeuf et al. (1979) to model plasma flow as an extension of the particle-in-cell (PIC) method developed by Harlow (1964).One of the main differences compared to SPH is that field quantities of fluid materials are accumulated as averages at pre-defined mesh points, and particles carry the quantities associated with fluid elements across a fixed background mesh.Besides, the dynamics of macro-particle quantities (such as motion, position, pressure, temperature) are directly governed by a macroscopic continuum model using quantities at the mesh points.In the macro-particle approach, if the particle size is chosen as the same size of real particles (e.g., single vehicle or pedestrian), this approach can be viewed as a mesoscopic simulation (Chang et al., 1985;Abdelghany et al., 2012).In this work, we would like to retain the categorization of the SPH approach as macroscopic simulation, since individual particles are used as media for numerical solution to the underlying macroscopic continuum model.They represent pieces of continuum and thus cannot be seen as real flow particles (or particle bunches).
Regarding numerical approximation, SPH has several advantages over grid-based methods that make it particularly interesting for macroscopic pedestrian simulation.Firstly, grid generation is unnecessary.In addition, it describes crowd dynamics using a Lagrangian formulation, which makes tracing traffic-state boundaries relatively easy.Finally, since SPH approximations are not affected by the arbitrariness of particle distribution, they can easily handle crowd deformation (i.e., parts of the crowd moving faster and/or in different directions than others).
The SPH method has been extended for applications to problems in fluid and solid mechanics, and a wide industrial context (Shadloo et al., 2016).Although there have been some attempts at solving microscopic pedestrian models using the SPH method (Tantisiriwat et al., 2007;Vetter et al., 2011), currently there exists no general simulation framework for applying this method to macroscopic pedestrian simulation.Here, we put forward a generic SPH framework for solving macroscopic crowd models and we show its potential for macroscopic simulation.The main contributions of this paper are as follows.First, the paper puts forward for the first time a generic framework for macroscopic pedestrian flow simulation underpinned by the SPH numerical method, where new behavioral assumptions can be easily added.Second, the paper presents an application of this framework to solve a specific macroscopic pedestrian flow model.Third, the paper shows that the proposed SPH method is numerically accurate.Fourth, a sensitivity analysis on key simulation variables elicits mathematical insights on the performance of the SPH method for this type of application.Fifth, the potential advantages of SPH over grid-based numerical methods are demonstrated by means of a comparison case study.Finally, the paper shows the method can reproduce key phenomena observed in real pedestrian flows, including shockwave propagation when a queue dissolves (Helbing et al., 2007;Portz and Seyfried, 2010), density dispersion when pedestrians enter a wide space (Hoogendoorn et al., 2014), and self-organized lane formation in bi-directional flows (Hoogendoorn and Daamen, 2004).
The structure of this paper is as follows.Section 2 describes the class of pedestrian flow models that we propose to solve in this paper using the SPH method.Section 3 describes the main principles of Smoothed Particle Hydrodynamics (SPH) as numerical method.Section 4 presents our SPH framework to simulate pedestrian flow at the macroscopic level using the model presented in Section 2. Next, Section 5 presents the setup of a series of simulation experiments that were performed to test mathematical properties and numerical advantages of the model, and to show that the SPH framework can produce meaningful and accurate simulation results.Section 6 presents the results of the simulation experiments.Finally, the key findings of this paper are reviewed in Section 7.

Macroscopic pedestrian flow model
This section presents the type of pedestrian flow model that we propose to solve numerically using the SPH method.First, in Section 2.1, we introduce the generic formulation of macroscopic pedestrian flow models, as described already in Jiang et al. (2010) and Twarogowska et al. (2014).Next, we specify a macroscopic pedestrian flow model based on that generic formulation (Section 2.2).The modeling principles are similar to those presented by Hoogendoorn et al. (2014), but they are formulated in different coordinate systems.Later on, the specified model will be used to evaluate the performance of the SPH method (see Section 5).

Generic macroscopic pedestrian flow model
Consider a large number of pedestrians moving through a two-dimensional continuous walking facility.The mass conservation equation must be satisfied at any point of the pedestrian flow.Here, this equation is written in Lagrangian form, where a Lagrangian coordinate system that moves with the motion of individual pedestrians is applied.This is in contrast to the Eulerian formation, which uses a coordinate system that is fixed in space (van Wageningen-Kessels et al., 2016).
where t and denote time and pedestrian density (i.e., number of pedestrians per unit area), respectively; x is the position vector (which can be decomposed into two dimensional components, x and x ); and v is the velocity vector (which can also be decom- posed into two dimensional components, v and v ).The symbol D Dt / denotes the material time derivative (i.e., the derivative taken with respect to a coordinate system that moves along with the pedestrian flow), which is defined as • .Furthermore, we assume that the acceleration of any point of the pedestrian flow is given by the following acceleration equation: where corresponds to the navigation term (which steers pedestrians towards their destination); and P denotes the pressure (a quantity that will be defined below); and E denotes external forces (e.g., interaction with wall boundaries).
Generally, in order to steer pedestrians towards their desired destination, the navigation term ( ) continually adjusts the velocity of any given point of the pedestrian flow so as to bring it closer to a certain equilibrium velocity (denoted by v e ).Given a fixed destination point ( x d ), the equilibrium velocity at any point x is determined by the local density ( ), which influences the walking speed of pedestrians at that point, and the relative position of point x with respect to x ), which determines the desired walking direction of these pedestrians.Therefore, given the destination, the equilibrium velocity of the pedestrian flow depends on the location ( x ) and the pedestrian density at that location ( ), that is to say, . So the navigation term ( ) is also function of those two variables, plus the current velocity, In physics, static pressure is defined as the amount of force exerted from all directions on a fluid particle per unit area.Here, we adapt this concept and define pressure as the psychological discomfort experienced by pedestrians as a result of their proximity to other pedestrians.In this light, we can define pressure as an increasing function of density ( = P P ( ), with > dP d / 0, ).Note that the pressure term P/ (see Eq. ( 2)) depends on the pressure gradient and, therefore, on the density gradient.Consequently, this term steers pedestrians towards low-density areas and away from high-density areas while they move towards their destination.In other words, whereas the navigation term ( ) accounts for global path choice behavior, the pressure term takes into account local path choice considerations (Hoogendoorn et al., 2015).
The last term ( E ) in Eq. (2) aims to introduce new physics (behavioral assumptions) to the flow, and it only activates for particular applications.This is one of the advantages of the SPH method, where new physics can be easily added in the governing equations without paying too much attention to numerical discretization.Two application examples will be further discussed in Sections 5.2 and 5.3, respectively.
Different assumptions regarding pedestrian behavior can lead to different formulations of P , , and E (and thus, of Eq. ( 2)), which in turn lead to the specification of different macroscopic pedestrian flow models (Twarogowska et al., 2014).Note, however, that the speed of pedestrian flows is limited by the physical characteristics of individual pedestrians.Depending on how Eq. ( 2) is specified, the acceleration given by this equation may not be realistic from a behavioral viewpoint under all circumstances.In particular, the model needs to satisfy the following constraint to ensure numerical stability: where u 0 is the free-flow speed (parameter).

Specific macroscopic pedestrian flow model
Let us develop a macroscopic pedestrian model based on the generic formulation presented in Section 2.1.In order to do that, we need to specify the acceleration equation (Eq.( 2)).We define the navigation term ( ) based on the concept of driving force proposed by Helbing and Molnár (1995): where is the relaxation time (parameter), which accounts for the fact that pedestrians generally do not accelerate or decelerate instantaneously to the equilibrium velocity but they do so over a certain time period.Thus, the navigation term represents a relaxation to the equilibrium velocity v e (see Eq. ( 4)).The equilibrium velocity can be defined as the product of the equilibrium speed (u e ), which is a scalar variable representing the magnitude of vector v e , and a unit vector corresponding to the desired walking direction ( ): (5) We assume that the equilibrium speed is only dependent on the local pedestrian density ( = u u ( ) e e ), and this relation can be described using a triangular fundamental diagram (van Wageningen-Kessels et al., 2016) where c and jam denote critical density and jam density (parameters), respectively.Furthermore, we define the desired walking direction as follows (Hoogendoorn et al., 2014): where v 0 is the desired velocity in homogeneous density conditions, when the desired path is a straight line to the destination point; and is the weight of the local route choice component ( 0).The desired velocity in homogeneous density conditions can be defined as follows: Thus, given the destination point, this variable is only dependent on the position, i.e., . It is important to remark that, as shown in Eqs. ( 7) and ( 8), our definition of assumes that the desired walking direction depends on the position within the pedestrian flow (global path choice) as well as on local traffic conditions (since depends on the density gradient), that is to say, ).Indeed, if > 0, the desired walking direction may deviate from the straight direction to the destination in order to avoid high-density areas.Parameter regulates the influence of local path choice behavior on the desired walking direction.In Hoogendoorn et al. (2014), it is shown that this formulation of can be derived from a microscopic pedestrian model.For comparison purpose regarding model simulation, a case study based on the macroscopic pedestrian model presented in Hoogendoorn et al. (2014) that is solved by a state-of-the-art grid-based numerical scheme will be used to further demonstrate some of the differences between the grid-based and meshfree approaches (Section 5.2).
Recall that, in the acceleration equation of the generic model (Eq.( 2)), the role of the pressure term P/ is to steer pedestrians towards low-density areas and away from high-density areas while they move towards their destination (see Section 2.1).However, the definition that we propose for the navigation term (see Eqs. ( 4)-( 8)) already takes into account this property of local path choice behavior (since = x ( , , )).For this reason, and for simplicity, we leave the pressure term out of the acceleration equation, so Eq. ( 2) becomes . Therefore, the acceleration equation is completely defined by Eqs. ( 4)-( 8), and corresponding expressions for external forces, if necessary.

Basic principles of the SPH method
The SPH method discretizes continuum media (e.g., water or, in our case, pedestrian flow) into a finite number of unconnected discrete elements, generally called particles.Each particle represents a finite mass of the discretized continuum and, as such, possesses material properties (such as density and velocity) and moves in accordance with mass conservation and acceleration laws.In order to determine how the material properties of every given particle evolve over time, the SPH method approximates the mass conservation and acceleration equations (in macroscopic pedestrian flow modeling, Eqs. ( 1) and ( 2)) through interpolation, using information about the material properties of the nearest neighboring particles.More specifically, the SPH method reformulates these equations by means of two approximations of the field functions ( f x ( )) that constitute the variables of these equations: a) kernel approximation; and b) particle approximation.The approximated variables are density ( ), velocity ( v ), pressure (P) and their spatial derivatives.Through these approximations, the partial differential equations that describe the mass conservation and acceleration laws are transformed into a set of ordinary differential equations that can be solved using a standard time-integration algorithm (Liu and Liu, 2009).The concepts of kernel and particle approximation are described in Sections 3.1 and 3.2, respectively.

Kernel approximation
Kernel approximation consists in transforming the relevant field functions and their spatial derivatives into integral equations by means of an interpolation function (or smoothing function), resulting in the so-called kernel estimates of the field functions and their derivatives at a given point x (Liu and Liu, 2009).Kernel estimates are calculated using the field function values at other points ( x ), which are weighted using the smoothing function (W).Therefore, the kernel estimates of a given field function f x ( ) and its spatial derivatives can be defined as: where h is the smoothing length (or radius), which defines the so-called support domain of the kernel estimate (i.e., the area around point x within which points x need to be located in order to be included in the integral).The value of the smoothing function depends on the distance between points x and x relative to h (in general, W decreases with this distance).Note that the SPH method requires that the smoothing function satisfies certain mathematical conditions (Monaghan, 2005).An important one is the so-called compact support condition: where is a constant that defines the effective (non-zero) area of the smoothing function at point x (i.e., the support domain at that point).Several types of smoothing functions are used in SPH applications, such as Gaussian, quadratic and cubic spline functions (Liu and Liu, 2009).The most common, however, is the piecewise cubic spline function proposed by Monaghan and Lattanzio (1985), which is defined as follows (see also Fig. 1): where = h 15/(7 ) d 2 in two-dimensional space problems; and R is the relative distance between points x and x ( = R x x h / ).

Particle approximation
As mentioned above, the SPH method discretizes continuum media into a set of unconnected particles.Each particle represents a finite mass of the discretized continuum; for instance, in macroscopic pedestrian flow modeling, these particles represent a certain mass of pedestrians (and not necessarily individual pedestrians).Particle approximation consists in transforming the continuous integral representations of the field functions and their derivatives (Eqs.( 9) and ( 10)) into discretized representations that can be applied to a computational domain divided into particles.The discretized representations of the kernel estimates of a field function and its spatial derivatives at a particle i (located at point x i ) can be defined as follows: where the smoothing function (W ij ) and its spatial derivative ( W i ij ) are defined as:  Monaghan and Lattanzio (1985) and its first order derivative.
In Eqs. ( 13)-( 15): = … j N 1, , is an index identifying all particles located within the support domain of particle i (which are usually called nearest neighboring particles); m j denotes the mass of particle j (which is constant); j is the density of particle j; and r ij is the distance between particles i and j ( = r x x ij i j , see Fig. 2a).These equations state that the value of (the gradient of) a field function at particle i can be approximated using the average of the values of that function at all nearest neighboring particles, weighted by the value of (the gradient of) the smoothing function, which depends on the distance between particles (see also Fig. 2a).The nearest neighboring particles are generally identified using a search algorithm.For an overview of search algorithms used in SPH applications, we refer to Liu and Liu (2009).
To summarize, the SPH method transforms the partial differential equations that describe the mass conservation and acceleration laws (in macroscopic pedestrian flow modeling, Eqs. ( 1) and ( 2)) by approximating field functions (e.g., density, velocity) and their derivatives using Eqs.( 13) and ( 14).This results in a set of ordinary differential equations that can be solved using a standard timeintegration algorithm (Liu and Liu, 2009).Next section presents the SPH approximations of the macroscopic pedestrian model equations introduced in Section 2.

SPH framework for macroscopic crowd simulation
SPH approximations of the mass conservation and acceleration equations for pedestrian flows (Eqs.( 1) and ( 2)) can be derived using Eqs.( 13)-( 15).The result is a set of ordinary differential equations describing how the density and velocity of a given particle i (which represents a finite mass of the discretized pedestrian flow) evolves over time as a function of the flow properties of the neighboring particles.Sections 4.1 and 4.2 presents the SPH approximations of the mass conservation and acceleration equations, respectively.Various SPH approximations of these equations are possible; here, we use symmetrized approximation equations, which are known to reduce numerical errors (Liu and Liu, 2009).

SPH approximation of the mass conservation equation
Eq. ( 1) can be rewritten as: The following identity can be derived from the product rule for derivatives: An SPH approximation of the expression on the right-hand side of Eq. ( 17) can be derived using Eq. ( 14): which can be rewritten as: where i is the index of the particle at which we are approximating , 1, , is the index identifying the nearest neighboring particles, and W x / ij i is the derivative of the smoothing function (W ij ) with respect to spatial dimension x .Substituting Eq. ( 19) into Eq.( 16), one obtains the SPH approximation of the mass conservation equation for pedestrian flows: (20)

SPH approximation of the acceleration equation
From the product rules for derivatives, it follows that: Furthermore, from the reciprocal rule for derivatives, it follows that: = ( ) Subsituting Eq. ( 22) into Eq.( 21) and rearranging the terms, we obtain the following equation: An SPH approximation of the expression on the right-hand side of Eq. ( 23) can be derived using Eq. ( 14): which can be rewritten as: Substituting Eq. ( 25) into Eq.( 2), one obtains the SPH approximation of the acceleration equation for pedestrian flows: where ij is the artificial viscosity.In SPH applications, it has often been observed during flow (material) compression, that numerical solutions present large unphysical oscillations and are unstable if a dissipative term is not introduced into the governing equations (Bui et al., 2008).Therefore, this term is included in SPH approximations of the acceleration function, and it needs only to be present in order to improve the stability of the simulation.It is defined by four parameters ( c , , , and ) and its formulation can be found in Liu and Liu (2009).However, recall that in our case, the pressure term is not included, so Eq. ( 26) can be rewritten as: The SPH approximation of the navigation term is defined as: The equilibrium speed and desired walking direction of particle i are defined by the following two equations: (29) In Eqs. ( 28) and ( 29), i an denotes the anisotropic density.Here, we define the fundamental diagram using the anisotropic density, so as to take into account that pedestrians adapt their speeds mainly according to the traffic conditions in front of them.We approximate this density by means of Eq. ( 13) and considering only the nearest neighboring particles that are located in the direction of motion of particle i (i.e., the direction of vector v i ), which we identify with index = … g G 1, , (see Fig. 2b): In Eq. ( 31), we multiply by two to account for the fact that we only consider the nearest neighboring particles in half of the support domain.This mainly accounts for the interaction with pedestrians located in front.More advanced formulations can be developed by adopting a constitutive model taking into account a more realistic anisotropic behavior of pedestrian flow.Also, we can implement more advanced smoothing functions where the boundary of the support domain is defined by e.g., a dimpled limacon (Helbing et al., 2002).However, this is beyond the scope of this work.
Note that Eq. ( 26) describes the dynamics of the first component of the velocity vector (v i ), which is denoted by v i ; the total derivative of component v i is calculated using an equivalent formula for the other dimension.Also, W x / ij i and W x / ij i are determined solely by the relative position of particles i and j (more specifically, x x r ( )/ i j ij and x x r ( )/ i ij , respectively), given a choice of smoothing function and value of h.The pressures P i and P j can be calculated using densities i and j based on function = P P ( ).

Simulation setup
We performed simulations using the proposed SPH framework (see Section 4) with the objective of showing that pedestrian flows can indeed be simulated at the macroscopic level using the SPH method and that this approach can yield qualitatively valid solutions.We considered three case studies.The first and second cases simulate uni-directional flow scenarios, whereas the third one simulates a bi-directional flow scenario.The first case evaluates the ability of the proposed SPH method to approximate the exact solution of the model equations (i.e., the numerical accuracy).It also analyzes the sensitivity of the method to changes in key parameters, such as the smoothing length (h) and the relaxation time ( ).The second case shows the ability of an SPH-based simulation model to reproduce density dispersion in wide spaces, a key phenomenon observed in real crowds.The results are compared to the output of a counterpart model which is solved using a state-of-the-art grid-based numerical scheme, in order to highlight key differences between the two methods.The third case aims to test the capability of the simulation framework to reproduce lane formation in bi-directional flows (also a phenomenon observed in real crowds).
Parameters of the triangular fundamental diagram (van Wageningen-Kessels et al., 2016) are given in Table 1.Case-dependent parameter specifications are shown in Table 2 and explained more in detail in Sections 5.1, 5.2, 5.3.In all three cases, we used the linked-list algorithm for nearest neighboring particle searching (i.e., to identify what particles are within the support domain of every particle every simulation time step).Furthermore, we used the Leap Frog (LF) algorithm for time integration (i.e., to update the density, velocity and position of every particle every simulation time step).The characteristics of these two algorithms are explained in Liu and Liu (2009).The smoothing function (W) is the one proposed by Monaghan and Lattanzio (1985) and defined in Eq. ( 12) (see also Fig. 1).

Case 1: Uni-directional flow, moving from standstill
In the first case, we simulate a crowd that is initially standing still at a high density (jam density: = 3.7 0 peds/m ) and has a rectangular shape (50 x 9 m).The left and right sides of the rectangle are the short ones.The pedestrians located at the right side of the rectangle start walking to the right (i.e., away from the crowd, thus no specific destinations are defined for each particle).The other pedestrians start to follow them (thus the process resembles the dissolution of a queue), as shown in Fig. 3a.There are no walls nor any other type of obstacle.

Table 1
Parameters of the fundamental diagram (Eq.( 29)).The main objective of this case is to determine whether the numerical solution obtained using the proposed simulation method approximates well the analytical solution calculated using kinematic wave theory (Lighthill and Whitham, 1955;Richards, 1956).To facilitate this comparison, we use a simplified version of our macroscopic model.Firstly, the relaxation time is set to the same value as the simulation time step length ( = = t 10 4 ), so as to allow particles to accelerate instantaneously to the equilibrium velocity (v e ), which is one of the characteristics of the kinematic wave model.In cases 2 and 3, is set to 0.6 s, a value suggested by Twarogowska et al. (2014).Secondly, the weight of the local route choice component is set to zero ( = 0) so as to avoid density dispersion and make sure that all particles move always in exactly the same direction (note that the kinematic wave model assumes one-dimensional flow movement).Finally, the artificial viscosity is not included in the acceleration equation.
Furthermore, this case is used to test the sensitivity of the SPH-based model to changes in several key parameters: • initial mesh size for particle generation ( x) with a range: 0.5, 1.0, and 1.5 m; • smoothing length (h) with a range: x x x x x x 1.2 , 1.5 , 2 , 2.5 , 3 , 3.5 and x 4 ; • relaxation time ( ) with a range: 10 4 , 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 s.
Based on the results from the sensitivity analysis, we select the initial mesh size x as 0.5 m, the smoothing length h as = x 2.5 1.25 m, for all three cases to ensure numerical accuracy and stability.The reasons behind this choice can be found in Section 6.1.Note that other parameters that are linked to the physical behavior of pedestrian flows should be defined specifically for the three case studies.The mass of each particle is constant throughout of the simulation, defined by x • 2 0 .

Case 2: Uni-directional flow, entering a wide space
In the second case, a crowd flows out of a 10-m wide corridor and enters a much wider space, as shown in Fig. 3b.The desired moving direction does not change when pedestrians exit the corridor (they continue moving towards the right).This case aims to show the influence of the weight of the local route choice component ( ) on the simulation results.Simulations are run with three  different values of (0, 1, and 2 m 4 /peds • s).In the scenarios with > 0, we expect that pedestrians will deviate a bit from their intended path to avoid high-density areas, causing density dispersion.A realistic model should be able to reproduce this phenomenon.
Note that in this test case we need to model the interaction of the crowd with the walls of the corridor.Pedestrians should neither touch nor go through these walls.In most SPH applications, boundaries such as walls are modeled using virtual particles.Interaction between particles that represent the fluid and virtual particles is modeled by including a repulsive boundary force ( B i ) in the fluid particle acceleration function (Liu and Liu, 2009).Here, we use the same approach, considering the boundary force B i as a form of external force E as defined in Eq. ( 2).
We model walls using arrays of virtual particles with constant zero velocity and constant non-zero density.Furthermore, the acceleration equation (Eq.( 27)) is changed to: where the boundary force is defined as a summation of the individual forces produced by all virtual particles = … k K 1, , , located within the support domain of particle i, thus: The pairwise boundary forces (more specifically, their dimensional components) are defined as follows (see Liu and Liu (2009)): where r ik denotes the distance between particle i and virtual particle is the interaction strength; r 0 is the interaction range (generally it is equal to the value of the initial mesh size); 1 and 2 are constants.The parameter values used in this test case are shown in Table 2.
From Eq. ( 34) it follows that if the distance between virtual particle k and particle i is equal or higher than the threshold r 0 ( 1 r rik 0 ), then particle k has no influence on the acceleration behavior of particle i ( = B 0 ik ).Instead, if virtual particle k is closer than the threshold distance r 0 to particle i ( > 1 r rik 0 ), then particle k produces a repulsive force on particle i that increases with decreasing distance between the particles based on a power function.
Note that virtual particles are included as nearest neighboring particles in the mass conservation and acceleration equations of other particles (thus contributing to their dynamics) if they are within their support domain.However, the position, velocity and density of virtual particles is not updated during the simulation (Liu and Liu, 2009).
For the purpose of comparing the performance of the SPH method to that of a grid-base numerical method, the macroscopic pedestrian model presented in Hoogendoorn et al. (2014), with the same definition for walking choice behavior, is used as the counterpart model.This model is formulated in Eulerian coordinate system, and it is solved by a state-of-the-art grid-based numerical scheme: the two-dimensional (2D) approximate Godunov scheme developed by van Wageningen-Kessels et al. (2018).In the counterpart simulation, the spatial area of the same size is divided into cells, time is divided into time steps of selected length.The number of pedestrians in each cell is calculated for each time step.The scheme determines the fluxes at the cell interfaces by taking the minimum of the maximum outflow from the cell and the maximum inflow into an adjacent cell at their interface.The total simulation time in both cases is set to 30 s.The same conditions (including simulation time step, initial density value, and free flow speed) apply to the two models.

Case 3: Bi-directional flow
In the third case, two groups of pedestrians move in opposite directions inside a 14-m wide corridor.The two groups meet inside the corridor and need to pass through each other in order to reach their destinations (a point directly opposite to the particle's initial position), as shown in Fig. 3c.By means of this case, we aim to determine whether the proposed simulation framework can reproduce lane formation in bi-directional flows, a self-organized structure observed in real crowds (Hoogendoorn and Daamen, 2004).As will be shown in Section 6.3, the basic model produces unrealistic result (gridlock).We propose to solve this by reformulating the acceleration equation (Eq.( 32)), adding a so-called seepage force or infiltration force (S i ), as done in Bui et al. (2007) for water-soil interaction.This force can be viewed as an interaction force or drag force between the two directional flows that penetrate each other.In general, this interaction force can be modeled in many forms (including social forces) as long as it can describe the physical behavior of the two opposite pedestrian flows in the process of lane formation.Note that similar to the boundary force (see Eqs. ( 32)-( 34)), the seepage force is treated as an external force, a concept that has already been included in the acceleration equation (Eq.( 2)).To summarize, in this study case, the particle acceleration equation is defined as follows: where an artificial viscosity term ij is included to improve numerical stability in conditions of flow compression (the chosen parameter values in Table 2 are considered reasonable (Liu and Liu, 2009)).The seepage force is defined based on Darcy's law (Biot, 1956): In Eq. ( 36), = … a M 1, , is an index given to all particles located within the support domain of particle i that are moving in the opposite direction of particle i, is a permeability factor (its value used in this test case is shown in Table 2), and µ i is the -di- mension component of unit vector µ i , which indicates the direction of the seepage force.This vector is calculated similarly to the vector representing the desired walking direction ( i ), but in this case the density gradient corresponds to the density gradient of particles moving in the opposite direction.Thus, in the vicinity of particles moving in the opposite direction, the seepage force gives an extra acceleration incentive to particles to go towards their destination but at the same time it steers them away from particles that are moving in the other direction.
Note that this formulation assumes that the seepage force increases if the density of particles a decreases (see Eq. ( 36)).This is intuitively consistent, since the porosity of the crowd moving in the opposite direction increases if its density decreases, which facilitates the infiltration process.

Results
This section presents the results of the three case studies described in Section 5.The main objectives of these case studies are to show that pedestrian flows can be simulated at the macroscopic level using the SPH method and that this approach can yield qualitatively valid solutions, as well as to test the mathematical properties of the method regarding the influence of key parameters, and to demonstrate some of its advantages compared to grid-based numerical methods.

Case 1: Numerical accuracy and sensitivity analysis
The simulation results presented in Fig. 4 show that the proposed numerical method can reproduce well how a queue dissolves according to kinematic wave theory.When the pedestrians located at the head of the queue start moving, the other pedestrians follow them and the queue gradually dissolves.As shown in Fig. 4a-e, a shockwave arises between the part of the crowd where pedestrians are standing still and the part of the crowd that is moving.This shockwave moves upstream and its speed approximates well the shockwave speed calculated analytically using the kinematic wave model (the analytical shockwave is shown in Fig. 4e by a black dashed line).As expected, the traffic states that arise as pedestrians accelerate to their desired speed coincide with the traffic states that define the congested branch of the fundamental diagram (see Fig. 4f).

5.
Examples of particle trajectory plots from the sensitivity analysis on initial mesh size versus smoothing length (the black dashed line represents the theoretical shockwave speed, given the fundamental diagram; color bar: speed in m/s; = ).Fig. 5 shows particle trajectory plots of simulations with different values for two key simulation parameters (initial mesh size x, and smoothing length h).This comparison shows that within reasonable value ranges, a fine resolution with smaller initial mesh size and wider smoothing length leads to more accurate simulation results compared to the analytical solutions.In particular, the smoothing length has a very important influence on accuracy.This is at the cost of more computational effort.Therefore, a trade-off between accuracy and efficiency needs to be considered for distinct applications.Note that for specific applications the selection of the initial mesh size is also constrained by the dimension of simulated cases.For example, considering a bottleneck case with a narrow door, it would be crucial to select a smaller initial mesh size than the door width.
In this work, we select the initial mesh size as 0.5 m for the demonstration of all simulation cases.The smoothing length is set as x 2.5 , since a more detailed analysis shows that with a h larger than this threshold, the gain in convergence (the consistency between the observed wave speed and the analytical wave speed) is marginal, and local information will be smoothing out the approximation domain, which can in turn influence the computational results.This set of simulation parameters to ensure the numerical stability of the SPH method is applied to the three study cases.
Fig. 6 gives the results of testing the sensitivity of the model parameter ( ).Fig. 6b-d indicate that an increase in leads to a decrease in the observed wave speed due to the slower reaction of pedestrians.Fig. 6a also shows a delay in traffic state moving to the equilibrium if a large relaxation time is incorporated.For practical applications, the calibration of this parameter needs to be performed with empirical data.Note that the case with = is only used to compare with the analytical results derived from kinematic wave theory regarding numerical accuracy.

Figs. 7a and 8a show that if
= there is no density dispersion and pedestrians adhere to the planned route described by the desired velocity in homogeneous density conditions ( v 0 ).This happens because pedestrians do not deviate in response to local variations in density (see Eq. ( 30)).Instead, if > the local variations in density affect path choice and cause density dispersion: in Figs.7d, g, 8d, and g, pedestrians divert from their intended path and move towards low-density regions while walking to their destination, as expected.The degree of density dispersion is regulated by parameter .With = dispersion becomes greater than with = (see Figs. 7d, g, 8d, and g).This study case shows that including the local route choice component in the equation used to calculate the desired walking direction (Eq.( 30)) yields qualitatively valid simulation results.For practical applications, the dispersion parameter would need to be calibrated using empirical data.
Graphical comparison of Figs. 7 and 8 clearly demonstrates different characteristics of SPH and the grid-based numerical scheme.As shown in the middle columns of Figs. 7 and 8, grids discretize the continuous space area and, as a result, the state boundaries are not clearly identifiable (there is no obvious transition between flow states).Instead, state transition can be distinguished by individual particle units in the result of the SPH model, as shown in the left columns and the overlaid results in the right due to its meshfree, Lagrangian, and particle-based nature.Also, numerical diffusion (van Wageningen-Kessels et al., 2016) is observed in the simulation output of the grid-based numerical scheme (Figs.7b, e, h, 8b, e, and h), around the flow state boundaries (shocks) and inside of the flows.However, this numerical artifact is constrained in the results of the SPH approach (Figs.7a, d, g, 8a, d, and g), particularly in the non-local-route-choice case (as in the first rows of the plots).This complies with the previous observation that continuum traffic flow models that are formulated and solved in Lagrangian coordinates have many numerical advantages over their Eulerian counterparts (Yuan et al., 2012;van Wageningen-Kessels et al., 2016).
Regarding computation efficiency, we can only draw qualitative conclusion.The same computer (with an i7-2720QM processor and 8 GB RAM) was used to perform the simulations using the two models (with the same simulation period and time step).The gain from the SPH method is marginal.The most important result of the comparison is that the SPH method shows its potential to improve pedestrian simulation regarding accuracy, since it generates less numerical diffusion and produces more clearly identifiable state boundaries, which is of importance for numerical solutions to continuum macroscopic models.

Case 3: Lane formation
The results show that if we use the basic acceleration equation (Eq.( 32)), particles cannot pass through the crowd moving in the opposite direction and a gridlock situation emerges (see Fig. 9).This happens because the equilibrium speed of particles (u e i , ) decreases to almost zero when the two groups of particles meet (because the pedestrian density increases to high values), which leads to a very low value for the navigation term and, thus, low D v Dt / i values, leading to very low particle velocities (see Eq. ( 35)).This gridlock situation is unrealistic assuming initial crowd densities of 1.0 ped/m 2 , which suggests that the macroscopic model needs to be modified in order to reproduce lane formation.
As mentioned in Section 5, we tested the effects of including a seepage force term in the acceleration equation (see Eq. ( 35)).This term incentivizes pedestrian flow particles to continue moving towards their desired destination even if they meet a large group of particles moving in the opposite destination, so they can pass through it instead of stopping because of the high pedestrian density.As shown in Fig. 10, a high enough permeability factor ( = 600 peds/m 4 /s) -thus a high enough seepage force-makes particles pass through the crowd moving in the opposite direction, leading to the formation of various lanes integrated by particles moving in the same direction.The lane formation process is very dynamic: the number of lanes and their spatiotemporal patterns vary considerably during the simulation (see Fig. 10).Note that this permeability factor would need to be tuned using empirical data for realistic model applications.Also, the authors acknowledge that there may be better ways to reformulate the acceleration equation in order to enable the model to reproduce lane formation in bi-directional flows.The seepage force approach proposed in this paper is just a first attempt in that direction.

Conclusions
The grid-based numerical methods generally used to solve macroscopic pedestrian flow models have various shortcomings that limit their applicability (for example, need for grid generation and remeshing, and unclear state boundaries).Smoothed Particle Hydrodynamics (SPH) is a meshfree numerical method that can potentially overcome some of these limitations.This paper presents for the first time a generic SPH framework for solving macroscopic crowd models and showed its potential for macroscopic simulation.In this respect, the paper shows that a macroscopic pedestrian flow model can be solved using the SPH method, and this leads to meaningful simulation results, by means of three simulation case studies.The first case shows that the proposed SPH method can approximate well the exact solution of the equations that define the macroscopic pedestrian model (i.e., it is numerically accurate).In addition, a sensitivity analysis provides insight into the mathematical properties of the method regarding the influence of three key model parameters.The comparison study in the second case with a grid-based counterpart model demonstrates that the meshfree method can potentially be more accurate regarding determining state boundaries.Furthermore, the results from the three case studies suggest that our method can reproduce some of the key phenomena observed in pedestrian flows, including shockwave propagation when a queue dissolves (case 1), density dispersion when pedestrians enter a wide space (case 2), and self-organized lane formation in bi-directional flows (case 3).For the latter, a reformulation of the acceleration equation was necessary.
As future work, the mathematical properties of the simulation method need to be further investigated, including the influence of the choice of smoothing function (to account for more realistic anisotropic pedestrian behavior), simulation time step length and artificial viscosity on model performance.Also, we propose to explore whether defining an adaptive smoothing length (e.g., dependent on density) can improve the accuracy of the simulation results, particularly in scenarios with high particle dispersion.Note that the current paper adopts a very traditional SPH method, further improvement in formulation is possible, such as using kernel normalization to improve the accuracy of SPH approximation.In this paper, we reproduced lane formation by including a seepage force term in the particle acceleration equation, as done for water-soil interaction in Bui et al. (2007).Further research is necessary to investigate whether this can better be achieved by modifying the navigation term or by introducing a form of social force that can describe and reproduce this phenomenon.We also consider it necessary to test the simulation model and develop its formulation further in order to allow reproducing more characteristics of pedestrian flows in different motion base scenarios (e.g., crossing flows, presence of obstacles, rounding a corner) and various density conditions.To achieve this, development of the underlying macroscopic crowd model may be needed (e.g., by inclusion of pedestrian anticipation).Finally, the performance of the proposed SPH method needs to be systematically validated and compared with that of other (grid-based) numerical methods -e.g., Jiang et al. (2010), Hoogendoorn et al. (2015) and van Wageningen-Kessels et al. ( 2016)-with respect to numerical accuracy, computational speed, and ability to reproduce key phenomena observed in real crowds, using empirical observations.

Fig. 2 .
Fig. 2. Illustration of SPH approximation.(a) Support domain of particle i given a smoothing function W with smoothing length h.The support domain is circular with radius equal to κh.Blue particles are inside the support domain of particle i, whereas red ones are not.(b) Calculation of anisotropic density.Only blue particles are included in the calculation.

Fig. 7 .
Fig. 7. Simulation results of case 2 at time t = 12.5 s (color bar: density in peds/m 2 ).Left column: SPH method; Middle column: 2D Godunov scheme (van Wageningen-Kessels et al., 2018); Right column: Overlaid results from the two nu.merical methods.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. Simulation results of case 2 at time t = 25.0 s (color bar: density in peds/m 2 ).Left column: SPH method; Middle column: 2D Godunov scheme (van Wageningen-Kessels et al., 2018); Right column: Overlaid results from the two nu.merical methods.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Simulation results of case 3 with the acceleration equation excluding the seepage force term (color bar: density in peds/m 2 ).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Simulation results of case 3 with the acceleration equation including the seepage force term (color bar: density in peds/m 2 ).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2
Case dependent settings.
*This density is the same for pedestrian flow particles and virtual particles.