# A Cellular Automaton Model for Pedestrians’ Movements Influenced by Gaseous Hazardous Material Spreading

Authors: J. Makmul
Journal: Modelling and Simulation in Engineering
 (2020)

---

## Abstract

A cellular automaton (CA) model is proposed to simulate the egress of pedestrians while gaseous hazardous material is spreading. The advection-diffusion with source term is used to describe the propagation of gaseous hazardous material. It is incorporated into the CA model. The navigation field in our model is determined by the solution of the Eikonal equation. The state transition of a pedestrian relies on the arrival time of cells in the Moore neighborhood. Numerical experiments are investigated in a room with multiple exits, and their results are shown.


---

## Body

## 1. Introduction

During the recent decades, the research on pedestrian flow has become an interesting and important issue to study. A large number of scientists from different research fields have paid attention to studying and modeling the pedestrians’ movements during evacuations. Pedestrian models can help planers and designers to build safe public places and provide important information to understand pedestrian dynamics.Presently, the main methods to study pedestrian evacuation are based on experiment and simulation modeling. Many pedestrian evacuation models have been investigated by researchers on different levels of description, such as on macroscopic models and on microscopic models. The macroscopic models usually apply to the case of large crowd and involve averaged quantities, in particular density, velocity, and energy. Examples of macroscopic models used for the pedestrian flow model can be found in [1, 2] for the first order macroscopic models (or scalar models) and in [3, 4] for second order macroscopic models. The microscopic models describe the time evolution of the position of each individual, addressed as a discrete particle. It mainly contains a social force model [5], an optimal-velocity model [6], a magnetic force model [7], cellular automata models [8], and a discrete choice model [9].The spreading of some gaseous hazardous material (e.g., smoke, a gas cloud) is one of the most important factors that have a great impact on an evacuation. Quite a number of researches focus on pedestrian evacuations under smoke and fire conditions. Zhao and Gao proposed an extended floor-field model to study the pedestrian evacuation under the influence of smoke diffusing. Smoke avoiding and herding behavior of pedestrians are observed in their model [10]. Nguyen et al. integrated the smoke effect and the blind evacuation strategy into the fire evacuation. The results were confirmed by empirical data from metro supermarket [11]. Zheng et al. studied the influence of the fire and the smoke on the movement of pedestrians. An extended Floor-Field model is carried out for study. Their simulation results show that the fire location in the room and the spreading rates of the fire and the smoke highly influence the pedestrian evacuation dynamics [12].In this work, we are interested on the evacuation of pedestrians during the spread of some gaseous hazardous material. A cellular automaton model is used to simulate pedestrian movements. It is combined with the advection-diffusion equation, which is applied to the gaseous hazardous material density. This equation is used in many applications in science and engineering for fluid motion, heat transfer, and flow of gas or pollutant [13]. We solve it numerically by the operator splitting method, which is an efficient approach to solve problems in multidimensions. For the navigation field in our model, the Eikonal equation is applied to attain the arrival time of each cell in the domain. A pedestrian chooses a cell in the Moore neighborhood according to the travelling time in the next time step. The Eikonal equation is also incorporated into the CA model and solved numerically by the fast marching method [14].The main objective of this paper is to include the advection-diffusion equation [13] and the Eikonal equation [15] into a cellular automaton model. The developed model is extended to describe the pedestrian flow, while hazard gas is propagating. The path field is determined by the arrival time of cells, which is obtained from the Eikonal equation. The arrival time of cells depends on the pedestrian density as well as on the hazard density. The influences of the pedestrian density and the hazard source location on the arrival time and evacuation time are investigated and discussed.The framework of this paper is organized as follows. In Section2, we present the cellular automaton model and explain a way to couple it with the advection-diffusion and the eikonal equation. The update rules and the numerical methods to solve the advection-diffusion equation and the Eikonal equation are prescribed in this section. Numerical experiments and results are demonstrated in Section 3. Finally, conclusions are made in Section 4.


## 2. Model

We consider the pedestrian evacuation problem of a multiple exits room where there is a source of gaseous hazard material inside the room. We assume that the smoke density has no effect on the pedestrians’ visibility and their health. All pedestrians have a global knowledge of the physical setting of the room. The room is divided into a uniform rectangular grid of cells. The size of a cell is taken as0.4m×0.4m, which is the typical space occupied by a person in a dense crowd [16]. The time domain is discretized into a series of t1,t2,…,tm,…, where m is an integer. Pedestrians can only move to an empty cell of the Moore neighborhood [17]. A pedestrian selects a cell in the Moore neighbours to move in the next time step according to the arrival time of the cells, which is computed from the following Eikonal equation [15]:(1)∇Tx=1Fx,x∈Ω,Tx=0,x∈Γ0⊂Ω,Front=Γt=xTx=t,where Ω is a simulation domain, Tx is the arrival time of a front crossing the point x, and Γ0 stands for an initial front or area in the domain where the pedestrians want to go. Γt is the front at time t. Fx≥0 is a moving speed of the front and depends on the position of x. It is set to(2)Fx=0.001,x∈Ωb,Uρx,x∈Ω/Ωb,where Ωb are areas that are obstructed by obstacles [15] or areas with high hazard gas density. U is the speed-density function. It describes the relationship between the speed and the density of pedestrians. Presently, there are many speed-density functions available. In our simulation, we choose the following [18].(3)Uρx=Umax1−ρxρmax,ρx=1πR2∑j,x−xj<R1,where Umax and ρmax are the maximum speed and the density of the pedestrians, respectively. R is the radius of a ball used to compute the density ρ. 1/πR2 is a circle area with radius R. In the proposed model, the gaseous hazard density is integrated into the CA model which is used to describe the pedestrians’ movements while hazard gas is spreading. The development of hazard material is expressed through the following linear advection-diffusion equation [13, 19, 20]:(4)∂C∂t+w⋅∇C=κd∇2C+Scs,t∈Ω×ℝ+,with the Dirichlet boundary conditions C=0 on ∂Ω, the diffusion constant κd>0, the velocity field w=w1,w2∈ℝ2, the source term Scs,t, and the simulation domain Ω⊂ℝ2. For simplification, we assume that the source of hazard emits gas at a constant rate Qcg/s from a single source point cs=xs,ys. Thus, the source term can be written as(5)Scs=Qcδx−xsδy−ys,where δ is the Dirac delta function given by(6)δx=1,x=0,0,x≠0.


### 2.1. Solving the Advection-Diffusion Equation

There are various numerical techniques in the literature for solving the problems of the advection-diffusion equation, for example, the finite difference method, finite element method, finite volume method, sprectral method, and method of lines. For our purpose, we apply the operator splitting method, which usually solves multidimensional problems efficiently. In two dimensions, it is split into two parts. Thex-direction and y-direction are handled separately over two time steps. First, we generate a grid and use the same size as pedestrian CA cells, i.e.,(7)xi=a+ihx,i=0,1,2,…,Nx,hx=0.4,yj=c+jhy,j=0,1,2,…,Ny,hy=0.4,where Ω=a,b×c,d. Then, we perform the operator splitting on the two-dimensional advection-diffusion equation by writing them as(8)Ci,j∗=Ci,jt−Δtw1∗∂C∗∂x−κd∂2C∗∂x2,(9)Ci,jt+1=Ci,j∗−Δtw2t+1∂Ct+1∂y−κd∂2Ct+1∂y2+ΔtSi,jt,where(10)∂C∗∂x=Ci,j∗−Ci−1,j∗Δx,∂2C∗∂x2=Ci+1,j∗−2Ci,j∗+Ci−1,j∗Δx2,∂Ct+1∂y=Ci,jt+1−Ci,j−1t+1Δy,∂2Ct+1∂y2=Ci,j+1t+1−2Ci,jt+1+Ci,j−1t+1Δy2.We assume that the hazard gas density on the boundaries is zero. Therefore, the Dirichlet boundary conditions are applied:(11)C0,y,t=0,CNx,y,t=0,Cx,0,t=0,Cx,Ny,t=0.Initial conditions are given by(12)Cx,y,0=Qcδx−xsδy−ys,where xs,ys is a hazard source point and Qc is gas concentration at initial time. First, we solve implicitly the equation (8) in x-direction and obtain Ci,j∗ for all grid points i,j. The value Ci,j∗ is used to solve equation (9) in y-direction. Finally, we receive the hazard gas density for the next time step Ci,jt+1 for all grid points i,j. For further details and the convergence of this method, we refer to the reference [21].


### 2.2. Solving the Eikonal Equation

Several methods exist to solve the Eikonal equation, for example, the fast marching method [22], fast marching level set method [14], fast sweeping method [23], and fast iterative method [24]. In our experiments, we employ the fast marching method, which is well known and efficient to solve the Eikonal equation (1). The solution of the Eikonal equation, which is obtained from the fast marching method, is constructed outwards from the boundary data, ordered from the smallest to the largest arrival time. The fast marching method used the upwind scheme to discretize the Eikonal equation.


### 2.3. Update Rule

The main algorithm for updating the pedestrians’ positions in each time step is as follows:Step 0(a)
Discretize the simulation domain into a rectangular grid with size0.4m×0.4m. Three types of grid points are generated. One grid is for pedestrian CA cells. One is for solving the Eikonal equation to receive the arrival time of each cell in the domain. Another one is for solving the advection-diffusion equation to obtain gas hazard density.(b)
Randomly distribute the pedestrians in the simulation domain. Each cell occupies only one occupant or none.(c)
The state of a cell, occupied by an individual, is assigned to 1. Empty cells receive the value 0.Step 1: each pedestrian stays within one current positionCP at time ts.Step 2: solve the advection-diffusion equation applying the operator splitting method with the same size of lattice spacing as the pedestrian cell in the CA model. The hazard density of each cell on the domain is obtained.Step 3: compute the pedestrian density of each cell in the domain.Step 4: solve the Eikonal equation using the fast marching method with the same size of lattice spacing as the pedestrian cell in the CA model.Tx=0, where x is a grid point on the exits. Fx=Umax1−ρx/ρmax is set for walkable cells and Fx=0.001 for obstacle cells or cells with high hazard density Cx≥0.05. The arrival times of each cell on the domain are received.Step 5: each pedestrian chooses randomly one cell in the Moore neighborhood, which has less or equal remaining travel time than his current cell, for the time stepts+1.Step 6 (parallel update):a conflict arises when two or more pedestrians attempt to move to the same cell. Pedestrians, who have no conflict with other, move to their selected cell. For pedestrians with conflicts, the chosen cell is randomly assigned to one of them with equal probability. The selected pedestrian moves to the interrelated cell, and the unselected pedestrians remain in their originalCP without moving until the next time step.Step 7: all pedestrians update theirCP for the time step ts+1.Step 8: update the state of each cell in the domain. The cell state is assigned to 1, if it is occupied by an occupant. Otherwise it is 0.Step 9: setts = ts+1 and return to Step 1 until ts=tend.Remark 1.
In Step 4, we set the threshold smoke concentration asCx≥0.05. It is assumed to be high enough to influence the movements of the pedestrians. This value reduces the visibility range of a pedestrian to about 6m−7m, and his speed is effected significantly, as described in the references [25, 26].


## 2.1. Solving the Advection-Diffusion Equation

There are various numerical techniques in the literature for solving the problems of the advection-diffusion equation, for example, the finite difference method, finite element method, finite volume method, sprectral method, and method of lines. For our purpose, we apply the operator splitting method, which usually solves multidimensional problems efficiently. In two dimensions, it is split into two parts. Thex-direction and y-direction are handled separately over two time steps. First, we generate a grid and use the same size as pedestrian CA cells, i.e.,(7)xi=a+ihx,i=0,1,2,…,Nx,hx=0.4,yj=c+jhy,j=0,1,2,…,Ny,hy=0.4,where Ω=a,b×c,d. Then, we perform the operator splitting on the two-dimensional advection-diffusion equation by writing them as(8)Ci,j∗=Ci,jt−Δtw1∗∂C∗∂x−κd∂2C∗∂x2,(9)Ci,jt+1=Ci,j∗−Δtw2t+1∂Ct+1∂y−κd∂2Ct+1∂y2+ΔtSi,jt,where(10)∂C∗∂x=Ci,j∗−Ci−1,j∗Δx,∂2C∗∂x2=Ci+1,j∗−2Ci,j∗+Ci−1,j∗Δx2,∂Ct+1∂y=Ci,jt+1−Ci,j−1t+1Δy,∂2Ct+1∂y2=Ci,j+1t+1−2Ci,jt+1+Ci,j−1t+1Δy2.We assume that the hazard gas density on the boundaries is zero. Therefore, the Dirichlet boundary conditions are applied:(11)C0,y,t=0,CNx,y,t=0,Cx,0,t=0,Cx,Ny,t=0.Initial conditions are given by(12)Cx,y,0=Qcδx−xsδy−ys,where xs,ys is a hazard source point and Qc is gas concentration at initial time. First, we solve implicitly the equation (8) in x-direction and obtain Ci,j∗ for all grid points i,j. The value Ci,j∗ is used to solve equation (9) in y-direction. Finally, we receive the hazard gas density for the next time step Ci,jt+1 for all grid points i,j. For further details and the convergence of this method, we refer to the reference [21].


## 2.2. Solving the Eikonal Equation

Several methods exist to solve the Eikonal equation, for example, the fast marching method [22], fast marching level set method [14], fast sweeping method [23], and fast iterative method [24]. In our experiments, we employ the fast marching method, which is well known and efficient to solve the Eikonal equation (1). The solution of the Eikonal equation, which is obtained from the fast marching method, is constructed outwards from the boundary data, ordered from the smallest to the largest arrival time. The fast marching method used the upwind scheme to discretize the Eikonal equation.


## 2.3. Update Rule

The main algorithm for updating the pedestrians’ positions in each time step is as follows:Step 0(a)
Discretize the simulation domain into a rectangular grid with size0.4m×0.4m. Three types of grid points are generated. One grid is for pedestrian CA cells. One is for solving the Eikonal equation to receive the arrival time of each cell in the domain. Another one is for solving the advection-diffusion equation to obtain gas hazard density.(b)
Randomly distribute the pedestrians in the simulation domain. Each cell occupies only one occupant or none.(c)
The state of a cell, occupied by an individual, is assigned to 1. Empty cells receive the value 0.Step 1: each pedestrian stays within one current positionCP at time ts.Step 2: solve the advection-diffusion equation applying the operator splitting method with the same size of lattice spacing as the pedestrian cell in the CA model. The hazard density of each cell on the domain is obtained.Step 3: compute the pedestrian density of each cell in the domain.Step 4: solve the Eikonal equation using the fast marching method with the same size of lattice spacing as the pedestrian cell in the CA model.Tx=0, where x is a grid point on the exits. Fx=Umax1−ρx/ρmax is set for walkable cells and Fx=0.001 for obstacle cells or cells with high hazard density Cx≥0.05. The arrival times of each cell on the domain are received.Step 5: each pedestrian chooses randomly one cell in the Moore neighborhood, which has less or equal remaining travel time than his current cell, for the time stepts+1.Step 6 (parallel update):a conflict arises when two or more pedestrians attempt to move to the same cell. Pedestrians, who have no conflict with other, move to their selected cell. For pedestrians with conflicts, the chosen cell is randomly assigned to one of them with equal probability. The selected pedestrian moves to the interrelated cell, and the unselected pedestrians remain in their originalCP without moving until the next time step.Step 7: all pedestrians update theirCP for the time step ts+1.Step 8: update the state of each cell in the domain. The cell state is assigned to 1, if it is occupied by an occupant. Otherwise it is 0.Step 9: setts = ts+1 and return to Step 1 until ts=tend.Remark 1.
In Step 4, we set the threshold smoke concentration asCx≥0.05. It is assumed to be high enough to influence the movements of the pedestrians. This value reduces the visibility range of a pedestrian to about 6m−7m, and his speed is effected significantly, as described in the references [25, 26].


## 3. Numerical Experiments and Results

We perform the following numerical experiments. The pedestrians are randomly distributed in the room at the initial time of the simulation and have to be evacuated from the room if there is a source of some gaseous hazardous material (e.g., smoke or a gas cloud). The modeling area is a rectangle16m×20m with two exits which are located at the bottom and on the right side of the room. The exits are labeled Exit 1 and Exit 2, respectively. The width of each exit is 2m, see Figure 1. The simulation domain is meshed into grid cells. Each cell has a size of 0.4m×0.4m which is the typical space occupied by a person in a dense area [27]. Thus, the average movement of an individual in each time step is 0.48m (parallel movement 0.4m or diagonal movement 0.42m). As the average velocity of pedestrians in a nervous state is about 1.65m/s [28, 29], one time step in the CA model is 0.29s. In order to obtain the general results, the numerical simulations are operated ten trial runs for the same example and their average is recorded. No re-entering mechanism was allowed. The evacuation process runs for 65 s. All simulation programs are implemented in MATLAB.Figure 1
The modeling area for the numerical experiments.To give an impression of the gaseous hazardous material that is spreading, we assume that the hazard density at the source point is rather high at the initial time and restrict our study to only one source. Furthermore, the hazard gas particles are supposed to generate from the source with a constant rate, i.e.,(13)Qc=10g/s,t=0,0.1g/s,t>0.The velocity fieldw1,w2 in equation (4) is chosen randomly out of the interval −0.5,0.5 at each time step. The maximum velocity Umax is set to 3 m/s, and the maximum density ρmax is 10 ped/m2. Figure 2 shows the snapshots of the spreading of the hazard where the source is located in the middle of the room and in front of Exit 1 at time intervals t = 5.22 s, 15.08 s, and 39.05 s. One can see that the gas particles, which are produced by the source of the hazard, spread around in the course of the experiment.Figure 2
The spread of gaseous hazardous material in scenarioi and scenario ii at time 5.22 s, 15.08 s, and 39.15 s. The velocity field w1,w2 in equation (4) is chosen randomly out of the interval −0.5,0.5. (a) t = 5 : 22 s, (b) t = 15 : 08 s, (c) t = 39 : 15 s, (d) t = 5 : 22 s, (e) t = 15 : 08 s, and (f) t = 39 : 15 s.
(a)
(b)
(c)
(d)
(e)
(f)First, we start with exploring the impact of a ball radiusR, which is used to compute the pedestrian density in equation (3) on the evacuation time. A ball radius is set to 1,2,…,10m for analyzing. The evacuation process begins with 300 pedestrians placed in a homogeneously distributed arrangement throughout the room. The pedestrians escape from the room, while the source of the gaseous hazardous material is located in front of Exit 1. We measure the mean evacuation time for ten trial runs using each radius. A plot of the mean evacuation time versus the ball radius R is shown in Figure 3. The average evacuation time is rather high when the ball radius is set to 1. As the ball radius increases from 1m to 4m, the evacuation time is reduced. The mean evacuation times are slightly different when the ball radius is more than 4. Therefore, the ball radius R influences the evacuation time of evacuees. Setting a very small radius R means pedestrians consider less pedestrian density for moving, leading to longer evacuation time. Thus, pedestrian density influences the arrival time of the cells and therefore the evacuation time.Figure 3
A plot of the ball radiusR versus the average evacuation time.Then, we examine the effect of the hazard source’s location on the evacuation process. The hazard source positions are considered in two cases: in the middle of the room in scenarioi, and it is located in front of Exit 1 in scenario ii. The number of evacuees in the experiments is 100, 300, and 500. The ball radius R in equation (3) is fixed to 4m for all trials. Table 1 displays the average evacuation time and the percentage of pedestrians who are moving out through Exit 1 considering the hazard spreading effect in scenario i, where the source is in the middle of the room and in scenario ii, where the source is in front of Exit 1. The hazard source’s place has a great impact on the evacuation time. When the source of danger occurs near the exit, it takes more time to evacuate all people compared with the situation of a source in the middle of a room. This can be interpreted as the fact that pedestrians reject to move out through the exit where the hazard source is located, since it refers to danger or insecurity for them to move out. Most of the pedestrians then use Exit 2 to escape. Therefore, the usage of Exit 2 is absolutely high and a large jamming is observed around this exit in scenario ii. The evacuation in scenario ii takes longer duration compared with scenario i. The usage of both exits is slightly different in scenario i, where the hazard source is in the centre of the room. The average number of outside pedestrians of each exit versus the time in Figure 4 is plotted for each scenario. The results show that more pedestrians evacuate through Exit 2 in scenario ii and the usage of the exits is quite balanced in scenario i.Table 1
The average evacuation time and the percentage of pedestrians who move out through Exit 1. The source is in the middle of the room in scenarioi and in front of Exit 1 in scenario ii.
Number of pedestrians
Evacuation times
Exit 1%
Scenarioi
Scenarioii
Scenarioi
Scenarioii
100
15.6020
s
21.1990
s
45.70
2.6
300
21.0250
s
34.5390
s
47.60
1.7
500
30.8850
s
56.2020
s
47.14
2.32Figure 4
Number of outside pedestrians of scenario (i) and scenario (ii). Total number of pedestrians is 300. (a) Source is in the middle and (b) source is near Exit 1.
(a)
(b)Table2 shows the computation time of ten trial runs in scenario i and in scenario ii with 100, 300, and 500 evacuees. The results in Table 2 reveal that the computation times for both scenarios increase as the number of pedestrians in simulations grows up. Scenario ii takes longer computation time than scenario i since the evacuation times in scenario ii are larger than those in scenario i, see Table 1. Large crowds are formed around Exit 2, and pedestrians take long time to move out of the room in scenario ii.Table 2
The computation time of ten trial runs in scenarioi and in scenario ii with 100, 300, and 500 evacuees.
Number of pedestrians
CPU time (10 rounds)
Scenarioi
Scenarioii
100
0.4116
h
0.7349
h
300
0.5132
h
0.8940
h
500
0.7130
h
1.2547
hFigures5 and 6 illustrate evacuees’ movements during the propagation of hazard, as the source is in the middle of the room and in front of Exit 1, respectively.Figure 5
The movements of 300 pedestrians during smoke spreading, where the source is in the middle of the room. The velocity fieldw1,w2 of the advection-diffusion equation (4) is chosen randomly out of the interval −0.5,0.5.Figure 6
The movements of pedestrians during smoke spreading, where the source is at Exit 1. The velocity fieldw1,w2 of the advection-diffusion equation (4) is chosen randomly out of the interval −0.5,0.5.From the results of scenarioii in Figure 6, we see that, at the beginning of the simulation, the pedestrians who are near Exit 1 move out through this exit due to less arrival time and the little gas density effect. As time increases, the hazard density around Exit 1 becomes higher, the pedestrians move away from it, and large crowds are formed around Exit 2. The clogging and the arching effects, which are collective phenomena naturally observed in pedestrian crowds, are remarkable in the model.The movements of 300 pedestrians during smoke spreading are demonstrated in Figure7. The velocity field of the advection-diffusion is set to w1,w2=0.5,0 in equation (4), and the source is in the middle of the room.Figure 7
The movements of 300 pedestrians during smoke spreading, as the velocity field of the advection-diffusion isw1,w2=0.5,0 in equation (4) and the source is in the middle of the room.The contour plots of the estimated remaining travel time (arrival time) of 300 pedestrians for the example in scenarioii and scenario i at t=15.08s are shown in Figure 8. Cells that belong to obstacle cells or high hazard gas density cells receive very large arrival time. In order to plot and see the difference of arrival times between scenario (i) and scenario (ii), we set the arrival times of these cells to 45 when the arrival time is larger than 45. The estimated arrival time around Exit 1 of scenario ii, where the hazard source is placed in front of Exit 1, is quite large compared with scenario i. This is due to fact that the hazard density is combined with the speed of the front Fx in equation (1). When the hazard density around Exit 1 is high, it leads to a large remaining time of cells around Exit 1. Therefore, pedestrians will reject to move to these cells, which surround Exit 1.Figure 8
Contour plots of the dynamic potential (estimated remaining travel time) of 300 pedestrians in scenarioii and scenario i at t=15.08s. (a) scenario (i) and (b) scenario (ii).
(a)
(b)


## 4. Discussion and Conclusions

In this paper, we incorporate the spread of hazard material into the cellular automaton model for the pedestrian evacuation in a room with multiple exits. The Eikonal equation is applied to determine the arrival time of each cell in the Moore neighborhood, and the evacuees select a cell to move in the next time according to this arrival time. Through this simple attempt, realistic pedestrians’ movements, while the hazard material is spreading, are demonstrated. The arching and the clogging effects, which are pedestrian collective phenomena occurring in pedestrian crowds around exits, were observed in the model [28,30]. Different numerical experiments are performed, and their results are shown.In this research, we limited our study with the assumptions that all pedestrians know the physical setting of the room and that the hazard gas density has no effect on the pedestrians’ visibility and their health. For future work, one might include these effects in modeling.



---

*Source: 3402198-2020-01-25.xml*