High fidelity modeling of aerosol pathogen propagation in built environments with moving pedestrians.

A high fidelity model for the propagation of pathogens via aerosols in the presence of moving pedestrians is proposed. The key idea is the tight coupling of computational fluid dynamics and computational crowd dynamics in order to capture the emission, transport and inhalation of pathogen loads in space and time. An example simulating pathogen propagation in a narrow corridor with moving pedestrians clearly shows the considerable effect that pedestrian motion has on airflow, and hence on pathogen propagation and potential infectivity.

INTRODUCTION

Advances in computational fluid and crowd dynamics (CFD, CCD), as well as computer hardware and software, have enabled fast and reliable simulations in both disciplines. A natural next step is the coupling of both disciplines. This would be of high importance for evacuation studies where fire, smoke, visibility and inhalation of toxic materials influence the motion of people, and where a large crowd can block or influence the flow in turn. The same capability could also be used to simulate with high fidelity the transmission of pathogens in the presence of moving pedestrians, enabling a much needed extension of current simulation technologies [56].

The present work considers a tight, bi‐directional coupling, whereby the flow (and any pathogens in it) and the motion of the crowd are computed concurrently and with mutual influences. Enabling technologies that made this tight coupling feasible include:

Development of immersed boundary methods;

Implementation of fast search techniques for information transfer between codes; and.

Strong scaling to tens of thousands of cores for CFD codes.

Before describing the numerical methodologies, a quick overview of pathogen, and in particular virus transmission is given. This defines the relevant physical phenomena, which in turn define the ordinary and partial differential equations that describe the flow and the particles. Thereafter, the models used for pedestrian motion are outlined. The sections on numerical methodologies conclude with a description of the coupling methodology employed. Several examples illustrate the influence of pedestrian motion on air motion, and hence aerosol and pathogen transport, and show the potential of the proposed methodology.

VIRUS INFECTION

Before describing the proposed model for virus transmission a brief description of virus propagation and lifetime is given. Viruses are usually present in the air or some surface, and make their way into the body either via inhalation (nose, mouth), ingestion (mouth) or attachment (eyes, hands, clothes). In many cases the victim inadvertently touches an infected surface or viruses are deposited on its hands, and then the hands or clothes touch either the nose, the eyes or the mouth, thus allowing the virus to enter the body.

An open question of great importance is how many viruses it takes to overwhelm the body's natural defense mechanism and trigger an infection. This number, the so‐called infectious dose, could depend on numerous factors, among them the state of immune defenses of the individual, the timing of viral entry (all at once, piece by piece), and the amount of hair and mucous in the nasal vessels [7]. In principle, a single organism in a favorable environment may replicate sufficiently to cause disease [65]. Data from research performed on biological warfare agents [20] suggests that both bacteria and viruses can produce disease with as few as 1–100 organisms (e.g. brucellosis 10–100, Q fever 1–10, tularaemia 10–50, smallpox 10–100, viral haemorrhagic fevers 1–10 organisms, tuberculosis 1). Compare these numbers and consider that as many as 3000 organisms can be produced by talking for 5 minutes or a single cough, with sneezing producing many more [43, 57, 58, 69, 75]. Figure 1, reproduced from [69], shows a typical number and size distribution. It should be pointed out that this is a histogram of the relative counts of difference particle sizes. There is often a large variation about the mean, which can impact disease risk estimates.

FIGURE 1

Counts of particles of various diameters in air expelled by 90 coughs [57]

VIRUS LIFETIME OUTSIDE THE BODY

Current evidence for Covid‐19 points to lifetimes outside the body that can range from 1–2 hours in air to several days on particular surfaces (so‐called fomite transmission mode) [38, 72]. There has also been some documentation of lifetime variation depending on humidity.

VIRUS TRANSMISSION

Sneezing and coughing

In the sequel, we consider sneezing and coughing as the main conduits of virus transmission. Clearly, breathing and talking will lead to the exhalation of air, and, consequently the exhalation of viruses from infected victims [1, 2, 3, 14]. However, it stands to reason that the size and amount of particles released ‐ and hence the amount of viruses in them—is much higher and much more concentrated when sneezing or coughing [2, 3, 18, 19, 37, 44, 69].

The velocity of air at a person's mouth during sneezing and coughing has been a source of heated debate, particularly in the media. The experimental evidence points to exit velocities of the order of 2–14 m/sec [23, 24, 66, 67]. A typical amount and size of particles can be seen in Figure 1.

Sink velocities

Table 1 lists the terminal sink velocities for water droplets in air based on the diameter [15, 61]. One can see that below diameters of the sink velocity is very low, implying that these particles remain in and move with the air for considerable time (and possibly distances).

TABLE 1

Sink velocities and Reynolds number for water particles in air

Diameter [mm]	Sink velocity [m/s]	Re
1.00E‐01	3.01E‐01	1.99E+00
1.00E‐02	3.01E‐03	1.99E‐03
1.00E‐03	3.01E‐05	1.99E‐06
1.00E‐04	3.01E‐07	1.99E‐09

Evaporation

Depending on the relative humidity and the temperature of the ambient air, the smaller particles can evaporate in milliseconds. However, as the mucous and saliva evaporate, they build a gel‐like structure that surrounds the virus, allowing it to survive. This implies that extremely small particles with possible viruses will remain infectious for extended periods of times ‐ up to an hour according to some studies [38, 72].

An important question is whether a particle/droplet will first reach the ground or evaporate. Figure 2, taken from [76], shows that below 0.12 the particles evaporate before falling 2 m (i.e. reaching the ground).

FIGURE 2

Evaporation time and falling time of droplets of varying diameter ( ). From [76]

PHYSICAL MODELING OF AEROSOL PROPAGATION

When solving the two‐phase equations, the air, as a continuum, is best represented by a set of partial differential equations (the Navier–Stokes equations) that are numerically solved on a mesh. Thus, the gas characteristics are calculated at the mesh points within the flowfield. The droplets/particles, which are relatively sparse in the flowfield, are modeled using a Lagrangian description, where individual particles (or groups of particles) are monitored and tracked in the flow, allowing for an exchange of mass, momentum and energy between the air and the particles.

Equations describing the motion of the air

As seen from the experimental evidence, the velocities of air encountered during coughing and sneezing never exceed a Mach‐number of . Therefore, the air may be assumed as a Newtonian, incompressible liquid, where buoyancy effects are modeled via the Boussinesq approximation. The equations describing the conservation of momentum, mass and energy for incompressible, Newtonian flows may be written as

Here denote the density, velocity vector, pressure, viscosity, gravity vector, coefficient of thermal expansion, temperature, reference temperature, specific heat coefficient and conductivity respectively, and momentum and energy source terms (e.g. due to particles or external forces/heat sources). For turbulent flows both the viscosity and the conductivity are obtained either from additional equations or directly via a large eddy simulation (LES) assumption through monotonicity induced LES (MILES) [6, 11, 21, 22, 30].

Equations describing the motion of particles/droplets

In order to describe the interaction of particles/droplets with the flow, the mass, forces and energy/work exchanged between the flowfield and the particles must be defined. As before, we denote for fluid (air) by and the density, pressure, temperature, conductivity, velocity in direction , viscosity, and the specific heat at constant pressure. For the particles, we denote by and the density, temperature, velocity in direction , equivalent diameter, specific heat coefficient and heat transferred per unit volume. In what follows we will refer to droplet and particles collectively as particles.

Making the classical assumptions that the particles may be represented by an equivalent sphere of diameter , the drag forces acting on the particles will be due to the difference of fluid and particle velocity:

The drag coefficient is obtained empirically from the Reynolds‐number : as (see, e.g. [61]):

The lower bound of is required to obtain the proper limit for the Euler equations, when . The heat transferred between the particles and the fluid is given by where is the film coefficient and the radiation coefficient. For the class of problems considered here, the particle temperature and kinetic energy are such that the radiation coefficient may be ignored. The film coefficient is obtained from the Nusselt‐Number : where is the Prandtl‐number of the gas as

Having established the forces and heat flux, the particle motion and temperature are obtained from Newton's law and the first law of thermodynamics. For the particle velocities, we have:

This implies that: where . The particle positions are obtained from:

The temperature change in a particle is given by: which may be expressed as: with . Equations (5.2.9), (5.2.10), (5.2.12) may be formulated as a system of Ordinary Differential Equations (ODEs) of the form: where denote the particle unknowns, the position of the particle and the fluid unknowns at the position of the particle.

Equations describing the motion of diluted viral loads

Viral loads may be obtained directly from the particles in the flowfield. An alternative for small, diluted particles that are floating in air is the use of a transport equation of the form:

Here denote pathogen concentration, the diffusivity and the source terms (due to exhalation or inhalation).

Numerical integration of the motion of the air

The last six decades have seen a large number of schemes that may be used to solve numerically the incompressible Navier–Stokes equations given by Equations (5.1.1), (5.1.2), (5.1.3)). In the present case, the following design criteria were implemented:

Spatial discretization using unstructured grids (in order to allow for arbitrary geometries and adaptive refinement);

Spatial approximation of unknowns with simple linear finite elements (in order to have a simple input/output and code structure);

Edge‐based data structures (for reduced access to memory and indirect addressing);

Temporal approximation using implicit integration of viscous terms and pressure (the interesting scales are the ones associated with advection);

Temporal approximation using explicit, high‐order integration of advective terms;

Low‐storage, iterative solvers for the resulting systems of equations (in order to solve large 3‐D problems); and.

Steady results that are independent from the timestep chosen (in order to have confidence in convergence studies).

The resulting discretization in time is given by the following projection Scheme [47, 48, 49]:

Advective‐Diffusive Prediction:

Pressure Correction:

which results in

Velocity Correction:

denotes the implicitness‐factor for the viscous terms ( : 1st order, fully implicit, : 2nd order, Crank‐Nicholson). are the standard low‐storage Runge–Kutta coefficients . The stages of Equation (5.4.2) may be seen as a predictor (or replacement) of by . The original right‐hand side has not been modified, so that at steady‐state , preserving the requirement that the steady‐state be independent of the timestep Δ. The factor denotes the local ratio of the stability limit for explicit timestepping for the viscous terms versus the timestep chosen. Given that the advective and viscous timestep limits are proportional to: we immediately obtain or, in its final form:

In regions away from boundary layers, this factor is , implying that a high‐order Runge–Kutta scheme is recovered. Conversely, for regions where , the scheme reverts back to the usual 1‐stage Crank‐Nicholson scheme. Besides higher accuracy, an important benefit of explicit multistage advection schemes is the larger timestep one can employ. The increase in allowable timestep is roughly proportional to the number of stages used (and has been exploited extensively for compressible flow simulations [36]). Given that for an incompressible solver of the projection type given by Equations (5.4.1), (5.4.2), (5.4.3), (5.4.4), (5.4.5), (5.4.6), (5.4.7)) most of the CPU time is spent solving the pressure‐Poisson system Equation (5.4.6), the speedup achieved is also roughly proportional to the number of stages used.

At steady state, and the residuals of the pressure correction vanish, implying that the result does not depend on the timestep Δ.

The spatial discretization of these equations is carried out via linear finite elements. The resulting matrix system is re‐written as an edge‐based solver, allowing the use of consistent numerical fluxes to stabilize the advection and divergence operators [49].

The energy (temperature) equation (Equation (5.1.3)) is integrated in a manner similar to the advective‐diffusive prediction (Equation (5.4.2)), i.e. with an explicit, high order Runge–Kutta scheme for the advective parts and an implicit, 2nd order Crank‐Nicholson scheme for the conductivity.

Numerical integration of the motion of particles/droplets

The equations describing the position, velocity and temperature of a particle (Equations 5.2.9, 5.2.10, 5.2.12) may be formulated as a system of nonlinear Ordinary Differential Equations of the form:

They can be integrated numerically in a variety of ways. Due to its speed, low memory requirements and simplicity, we have chosen the following k‐step low‐storage Runge–Kutta procedure to integrate them:

For linear ODEs the choice leads to a scheme that is ‐th order accurate in time. Note that in each step the location of the particle with respect to the fluid mesh needs to be updated in order to obtain the proper values for the fluid unknowns. The default number of stages used is . This would seem unnecessarily high, given that the flow solver is of second‐order accuracy, and that the particles are integrated separately from the flow solver before the next (flow) timestep, i.e. in a staggered manner. However, it was found that the 4‐stage particle integration preserves very well the motion in vortical structures and leads to less 'wall sliding' close to the boundaries of the domain [53]. The stability/ accuracy of the particle integrator should not be a problem as the particle motion will always be slower than the maximum wave speed of the fluid (fluid velocity).

The transfer of forces and heat flux between the fluid and the particles must be accomplished in a conservative way, i.e. whatever is added to the fluid must be subtracted from the particles and vice‐versa. The finite element discretization of the fluid equations will lead to a system of ODE's of the form: where and denote, respectively, the consistent mass matrix, increment of the unknowns vector and right‐hand side vector. Given the 'host element' of each particle, i.e. the fluid mesh element that contains the particle, the forces and heat transferred to are added as follows:

Here denotes the shape‐function values of the host element for the point coordinates , and the sum extends over all elements that surround node . As the sum of all shape‐function values is unity at every point: this procedure is strictly conservative.

From Equations (5.2.9), (5.2.10), (5.2.12) and their equivalent numerical integration via Equation (5.5.2), the change in momentum and energy for one particle is given by:

These quantities are multiplied by the number of particles in a packet in order to obtain the final values transmitted to the fluid. Before going on, we summarize the basic steps required in order to update the particles one timestep:

Initialize Fluid Source‐Terms: .

DO: For Each Particle:

DO: For Each Runge–Kutta Stage:

Find Host Element of Particle: IELEM,

Obtain Fluid Variables Required

Update Particle: Velocities, Position, Temperature, …

ENDDO

Transfer Loads to Element Nodes

ENDDO

Particle parcels

For a large number of very small particles, it becomes impossible to carry every individual particle in a simulation. The solution is to:

a) Agglomerate the particles into so‐called packets of particles;

b) Integrate the governing equations for one individual particle; and.

c) Transfer back to the fluid times the effect of one particle.

Beyond a reasonable number of particles per element (typically >8), this procedure produces accurate results without any deterioration in physical fidelity [53].

Other particle numerics

In order to achieve a robust particle integrator, a number of additional precautions and algorithms need to be implemented. The most important of these are:

Agglomeration/Subdivision of Particle Parcels: As the fluid mesh may be adaptively refined and coarsened in time, or the particle traverses elements of different sizes, it may be important to adapt the parcel concentrations as well. This is necessary to ensure that there is sufficient parcel representation in each element and yet, that there are not too many parcels as to constitute an inefficient use of CPU and memory.

Limiting During Particle Updates: As the particles are integrated independently from the flow solver, it is not difficult to envision situations where for the extreme cases of very light or very heavy particles physically meaningless or unstable results may be obtained. In order to prevent this, the changes in particle velocities and temperatures are limited in order not to exceed the differences in velocities and temperature between the particles and the fluid [53].

Particle Contact/Merging: In some situations, particles may collide or merge in a certain region of space.

Particle Tracking: A common feature of all particle‐grid applications is that the particles do not move far between timesteps. This makes physical sense: if a particle jumped ten gridpoints during one timestep, it would have no chance to exchange information with the points along the way, leading to serious errors. Therefore, the assumption that the new host elements of the particles are in the vicinity of the current ones is a valid one. For this reason, the most efficient way to search for the new host elements is via the vectorized neighbor‐to‐neighbor algorithm described in [45, 49].

Immersed body techniques

The information required from CCD codes consists of the pedestrians in the flowfield, i.e. their position, velocity, temperature, as well inhalation and exhalation. As the CCD codes describe the pedestrians as points, circles or ellipses, a way has to be found to transform this data into 3‐D objects. Two possibilities have been pursued here:

Transform each pedestrian into a set of (overlapping) spheres that approximate the body with maximum fidelity with the minimum amount of spheres;

Transform each pedestrian into a set of tetrahedra that approximate the body with maximum fidelity with the minimum amount of tetrahedra.

The reason for choosing spheres or tetrahedra is that one can perform the required interpolation/ information transfer much faster than with other methods.

In order to 'impose' on the flow the presence of a pedestrian the immersed boundary methodology is used. The key idea is to prescribe at every CFD point covered by a pedestrian the velocity and temperature of the pedestrian. For the CFD code, this translates into an extra set of boundary conditions that vary in time and space as the pedestrians move. This is by now a mature technology (see, e.g. chapter 18 in [49] and the references cited therein). Fast search techniques as well as extensions to higher order boundary conditions may be found in [49, 50]. Nevertheless, as the pedestrians potentially change location every timestep, the search for and the imposition of new boundary conditions can add a considerable amount of CPU as compared to 'flow‐only' runs.

MODELING OF PEDESTRIAN MOTION

The modeling of pedestrian motion has been the focus of research and development for more than two decades. If one is only interested in average quantities (average density, velocity), continuum models [29] are an option. For problems requiring more realism, approaches that model each individual are required [70]. Among these, discrete space models (such as cellular automata [4, 5, 12, 17, 39, 40, 42, 62, 68]), force‐based models (such as the social force model [27, 28, 41, 51, 60]) and agent‐based techniques [13, 25, 26, 59, 64, 71, 73] have been explored extensively. Together with insights from psychology and neuroscience (e.g. [71, 74]) it has become clear that any pedestrian motion algorithm that attempts to model reality should be able to mirror the following empirically known facts and behaviors:

Newton's laws of motion apply to humans as well: from one instant to another, we can only move within certain bounds of acceleration, velocity and space;

Contact between individuals occurs for high densities; these forces have to be taken into account;

Humans have a mental map and plan on how they desire to move globally (e.g. first go here, then there, etc.);

Human motion is therefore governed by strategic (long term, long distance), tactical (medium term, medium distance) and operational (immediate) decisions;

In even moderately crowded situations ( ), humans have a visual horizon of , and a perception range of 120 degrees; thus, the influence of other humans beyond these thresholds is minimal;

Humans have a 'personal comfort zone'; it is dependent on culture and varies from individual to individual, but it cannot be ignored;

Humans walk comfortably at roughly 2 paces per second (frequency: ); they are able to change the frequency for short periods of time, but will return to 2 whenever possible.

We remark that many of the important and groundbreaking work cited previously took place within the gaming/visualization community, where the emphasis is on 'looking right'. Here, the aim is to answer civil engineering or safety questions such as maximum capacity, egress times under emergency, or comfort. Therefore, comparisons with experiments and actual data are seen as essential [31, 32, 51].

The PEDFLOW model

The PEDFLOW model [51] incorporates these requirements as follows: individuals move according to Newton's laws of motion; they follow (via will forces) 'global movement targets'; at the local movement level, the motion also considers the presence of other individuals or obstacles via avoidance forces (also a type of will force) and, if applicable, contact forces. Newton's laws: where denote, respectively, mass, velocity, position, force and time, are integrated in time using a 2nd order explicit timestepping technique. The main modeling effort is centered on . In the present case the forces are separated into internal (or will) forces [I would like to move here or there] and external forces [I have been hit by another pedestrian or an obstacle]. For the sake of completeness, we briefly review the main forces used. For more information, as well as verification and validation studies, see [31, 32, 33, 34, 35, 51, 54, 77].

Will force

Given a desired velocity and the current velocity , this force will be of the form

The modeling aspect is included in the function , which, in the non‐linear case, may itself be a function of . Suppose is constant, and that only the will force is acting. Furthermore, consider a pedestrian at rest. In this case, we have: which implies: and

One can see that the crucial parameter here is the 'relaxation time' which governs the initial acceleration and 'time to desired velocity' Typical values are and . The 'relaxation time' is clearly dependent on the fitness of the individual, the current state of stress, desire to reach a goal, climate, signals, noise, etc. Slim, strong individuals will have low values for , whereas fat or weak individuals will have high values for . Furthermore, dividing by the mass of the individual allows all other forces (obstacle and pedestrian collision avoidance, contact, etc.) to be scaled by the 'relaxation time' as well, simplifying the modeling effort considerably.

The direction of the desired velocity will depend on the type of pedestrian and the cases considered. A single individual will have as its goal a desired position that he would like to reach at a certain time . If there are no time constraints, is simply set to a large number. Given the current position , the direction of the velocity is given by where denotes the desired position (location, goal) of the pedestrian at the desired time of arrival . For members of groups, the goal is always to stay close to the leader. Thus, becomes the position of the leader. In the case of an evacuation simulation, the direction is given by the gradient of the perceived time to exit to the closest perceived exit:

The magnitude of the desired velocity ∣ depends on the fitness of the individual, and the motivation/urgency to reach a certain place at a certain time. Pedestrians typically stroll leisurely at 0.6 − 0.8 m/s, walk at 0.8 − 1.0 m/s, jog at 1.0 − 3.0 m/s, and run at 3.0 − 10.0 m/s.

Pedestrian avoidance forces

The desire to avoid collisions with other individuals is modeled by first checking if a collision will occur. If so, forces are applied in the direction normal and tangential to the intended motion. The forces are of the form: where , denote the positions of individuals , the radius of individual , and . Note that the forces weaken with increasing non‐dimensional distance . For years we have used , but, obviously, this can (and probably will) be a matter of debate and speculation (perhaps a future experimental campaign will settle this issue). In the far range, the forces are mainly orthogonal to the direction of intended motion: humans tend to move slightly sideways without decelerating. In the close range, the forces are also in the direction of intended motion, in order to model the slowdown required to avoid a collision.

Wall avoidance forces

Any pedestrian modeling software requires a way to input geographical information such as walls, entrances, stairs, escalators, etc. In the present case, this is accomplished via a triangulation (the so‐called background mesh). A distance to walls map (i.e. a function is constructed using fast marching techniques on unstructured grids), and this allows to define a wall avoidance force as follows:

Note that ∣. The default for the maximum wall avoidance force is . The desire to be far/close to a wall also depends on cultural background.

Contact forces

When contact occurs, the forces can increase markedly. Unlike will forces, contact forces are symmetric. Defining these forces are modeled as follows: and .

Motion inhibition

A key requirement for humans to move is the ability to put one foot in front of the other. This requires space. Given the comfortable walking frequency of , one is able to limit the comfortable walking velocity by computing the distance to nearest neighbors and seeing which one of these is the most 'inhibiting'.

Psychological factors

The present pedestrian motion model also incorporates a number of psychological factors that, among the many tried over the years, have emerged as important for realistic simulations. Among these, we mention:

Determination/Pushiness: it is an everyday experience that in crowds, some people exhibit a more polite behavior than others. This is modeled in PEDFLOW by reducing the collision avoidance forces of more determined or 'pushier' individuals. Defining a determination or pushiness parameter , the avoidance forces are reduced by (1 − .

Comfort zone: in some cultures (northern Europeans are a good example) pedestrians want to remain at some minimum distance from contacting others. This comfort zone is an input parameter in PEDFLOW, and is added to the radii of the pedestrians when computing collisions avoidance and pre‐contact forces.

Right/Left Avoidance and Overtaking: in many western countries pedestrians tend to avoid incoming pedestrians by stepping towards their right, and overtake others on the left. However, this is not the norm everywhere, and one has to account for it.

Numerical integration of the motion of pedestrians

The equations describing the position and velocity of a pedestrian may be formulated as a system of nonlinear Ordinary Differential Equations of the form:

These ODEs are integrated with explicit Runge–Kutta schemes, typically of order 2.

The geographic information required, such as terrain data (inclination, soil/water, escalators, obstacles, etc.), climate data (temperature, humidity, sun/rain, visibility), signs, the location and accessibility of guidance personnel, as well as doors, entrances and emergency exits is stored in a so‐called background grid consisting of triangular elements. This background grid is used to define the geometry of the problem. At every instance, a pedestrian will be located in one of the elements of the background grid. Given this 'host element' the geographic data, stored at the nodes of the background grid, is interpolated linearly to the pedestrian. The closest distance to a wall or exit(s) for any given point of the background grid evaluated via a fast ( ) nearest neighbor/heap list technique ([49, 51]). For cases with visual or smoke impediments, the closest distance to exit(s) is recomputed every few seconds of simulation time.

Linkage to CFD Codes

The information required from CFD codes such as FEFLO consists of the spatial distribution of temperature, smoke, other toxic or movement impairing substances in space, as well as pathogen distribution. This information is interpolated to the (topologically 2‐D) background mesh at every timestep in order to calculate properly the visibility/ reachability of exits, routing possibilities, smoke, toxic substance or pathogen inhalation, and any other flowfield variable required by the pedestrians. As the tetrahedral grid used for the CFD code and the triangular background grid of the CCD code do not change in time, the interpolation coefficients need to be computed just once at the beginning of the coupled run. While the transfer of information from CFD to CCD is voluminous, it is very fast, adding an insignificant amount to the total run‐times.

COUPLING METHODOLOGY

The coupling methodology used is shown in Figure 3. The CFD code computes the flowfield, providing such information as temperature, smoke, toxic substance and pathogen concentration, and any other flow quantity that may affect the movement of pedestrians. These variables are then interpolated to the position where the pedestrians are, and are used with all other pertinent information (e.g. will‐forces, targets, exits, signs, etc.) to update the position, velocity, inhalation of smoke, toxic substances or pathogens, state of exhaustion or intoxication, and any other pertinent quantity that is evaluated for the pedestrians. The position, velocity and temperature of the pedestrians, together with information such as sneezing or exhaling air, is then transferred to the CFD code and used to modify and update the boundary conditions of the flowfield in the regions where pedestrians are present.

FIGURE 3

Coupling CFD and CCD Codes

Of the many possible coupling options (see e.g. [8, 9, 63]), we have implemented the simplest one: loose coupling with sequential timestepping ([46, 52, 55, 60, 63]). This is justified, as the timesteps of both the flow and pedestrian solvers are very small, so that possible coupling errors are negligible. PEDFLOW typically runs with fixed timesteps of Δ, while the timestep chosen by FEFLO depends on mesh size and velocity. Should the timestep of FEFLO be less than the default value for PEDFLOW, then PEDFLOW automatically reduces its timestep to be the same as FEFLO. During the course of many cases run, we have never encountered any stability problems with this loose coupling and timestepping strategy.

EXAMPLES

Corridor with pedestrians

This example considers the corridor of 10.0 m x 2.0 m x 2.5 m shown in Figure 4. Both entry and exit sides have two doors each of size 0.8 m x 2.0 m. For climatisation, 4 entry vanes and 1 exit vane are placed in the ceiling. The vertical air velocity for the entry vanes as set to , while the horizontal velocity was set as increasing proportional to the distance of the center of the vane to a maximum of . The CFD mesh had approximately 3.5 Mels, and as can be seen from Figures 5 and 6 had a uniform element size. Two streams of pedestrians enter and exit through the doors over time. As stated before, the pedestrian dynamics code, which only 'sees' the floorplan (Figure 7) of the problem at hand, computes position, velocity and orientation of the pedestrians, and then produces a tetrahedral mesh for each pedestrian (Figure 8) and then sends this information to the transfer library. This information is then passed on to the flow solver, which treats the pedestrians via the immersed body approach in the flowfield. Should there be smoke, pollutants or pathogens in the flowfield, this information is passed back to the pedestrian dynamics code, which interpolates it at the height of pedestrians in order to update inhalation, intoxication and infection information.

FIGURE 4

Simple Channel: Geometry and boundary conditions

FIGURE 5

Simple Channel: Y‐Plane‐cut of CFD Mesh

FIGURE 6

Simple Channel: Z‐Plane‐cut of CFD Mesh

FIGURE 7

Simple Channel: CCD background Mesh

FIGURE 8

Simple Channel: Surface of tetrahedral elements for immersed BC of pedestrians in CFD Mesh

This simulation was run in 3 phases:

Phase 1: Every code is run independently until 'things settle down', i.e until the flow reaches a quasi‐steady state and the pedestrian streams have formed; for the present case this took 20 sec of physical time;

Phase 2: The restart files from Phase 1 are taken, and the run continues in fully coupled mode, until 'things settle down'; for the present case this took 20 sec of physical time;

Phase 3: The restart files from Phase 2 are taken, and the run continues in fully coupled mode imposing the boundary conditions for a sneezing event.

In order to see the effects of pedestrians, three simulations were carried out: a) Two pedestrians streams in counterflow mode; b) Two pedestrians streams in parallel flow mode; and c) No pedestrians.

Figures 9A‐G, 10A‐G and 11A‐G show the solutions obtained for these different modes.

FIGURE 9

A, Counterflow Movement: Solution at . B, Counterflow Movement: Solution at . C, Counterflow Movement: Solution at . D, Counterflow Movement: Solution at . E, Counterflow Movement: Solution at . F, Counterflow Movement: Solution at . G, Counterflow Movement: Solution at

FIGURE 10

A, Parallel Movement: Solution at . B, Parallel Movement: Solution at . C, Parallel Movement: Solution at . D, Parallel Movement: Solution at . E, Parallel Movement: Solution at . F, Parallel Movement: Solution at . G, Parallel Movement: Solution at

FIGURE 11

A, No Pedestrians: Solution at . B, No Pedestrians: Solution at . C, No Pedestrians: Solution at . D, No Pedestrians: Solution at . E, No Pedestrians: Solution at . F, No Pedestrians: Solution at . G, No Pedestrians: Solution at

In the cases shown different temporal scales appear:

The fast, ballistic drop of the larger ( ) particles, occurring in the range of ;

The slower drop of particles of diameter , occurring in the range of ; and

The transport of the even smaller particles through the air, occurring in the range of .

We have attempted to show these phases in the results, and for this reason the results are not displayed at equal time intervals. Unless otherwise noted, the particles have been colored according to the logarithm of the diameter, with red colors representing the largest and blue the smallest particles.

Note the very large differences in the flowfield with and without pedestrians. The main reason for this difference is the discrepancy in velocities: humans walk at approximately , while the perception of discomfort due to air motion being at around , implying that in most of the volume of any built environment where humans reside, these lower velocities will be encountered. The different velocities between walking pedestrians, and in particular counterflows as the one shown, lead to large‐scale turbulent mixing, enhancing the spread of pathogens emanating from infected victims.

Note the very large difference in mixing and viral transmission due to the presence of pedestrians.

CONCLUSIONS AND OUTLOOK

A high fidelity model for the propagation of pathogens via aerosols in the presence of moving pedestrians has been implemented. The key idea is the tight coupling of computational fluid dynamics and computational crowd dynamics in order to capture the emission, transport and inhalation of pathogen loads in space and time in the presence of moving pedestrians.

The example of a narrow corridor with moving pedestrians clearly shows the considerable effect that pedestrian motion has on airflow, and hence on pathogen propagation and potential infectivity.

At present, the 'pedestrians' appear in the flow code as rigid bodies. The incorporation of leg and arm movement while walking would be a possible improvement to the model.