Modeling and simulation of the infection zone from a cough.

The pandemic of 2020 has led to a huge interest of modeling and simulation of infectious diseases. One of the central questions is the potential infection zone produced by a cough. In this paper, mathematical models are developed to simulate the progressive time-evolution of the distribution of locations of particles produced by a cough. Analytical and numerical studies are undertaken. The models ascertain the range, distribution and settling time of the particles under the influence of gravity and drag from the surrounding air. Beyond qualitative trends that illustrate that large particles travel far and settle quickly, while small particles do not travel far and settle slowly, the models provide quantitative results for distances travelled and settling times, which are needed for constructing social distancing policies and workplace protocols. © Springer-Verlag GmbH Germany, part of Springer Nature 2020.

Introduction

The pandemic of 2020, due to SARS-CoV-2, named COVID-19 and referred to as coronavirus, has been responsible for hundreds of thousands of deaths in 2020 alone. It is well-established that this virus primarily spreads from person-to-person contact by respiratory droplets produced when an infected person coughs or sneezes. Subsequently, the droplets come into contact with the eyes, nose or mouth of a nearby person or when a person touches an infected surface, then makes contact with their eyes, nose or mouth. Since the virus is small, 0.06–0.14 microns in diameter, it can be contained in or attached to such emitted droplets. Droplets as small as one micron can carry enough viral load to cause an infection. A particular concern is the interaction of droplets with ventilation systems, which potentially could enhance the propagation of pathogens. This has implications on situation-specific safe distancing and the design of building filtration systems, air distribution, heating, air-conditioning and decontamination systems, for example using UV-c and related technologies. In order to facilitate such system redesigns, fundamental analysis tools are needed that are easy to use. Accordingly, this paper develops one type of such needed tools, namely a simulator for the analysis of cough particle tracking, in order to ascertain how large is the potential infection zone and the airborne setting time of cough particles.

In its most basic form, a cough can be considered as a high-velocity release of a random distribution of particles of various sizes, into an ambient atmosphere. We refer the reader to Wei and Li [53], Duguid [11], Papineni and Rosenthal [37], Wei and Li [54], Zhu et al. [59], Chao et al. [9], Morawska et al. [31], VanSciver et al. [49], Kwon et al. [21], Tang et al. [46], Xie et al. [57], Gupta et al. [13], Wan et al. [52], Villafruela et al. [50], Nielson [33], Zhang and Li [58], Chao et al. [8] and Lindsley et al. [26] for extensive reviews of coughs and other respiratory emissions. Following formulations for physically similar problems associated with particulate dynamics from the fields of blasts, explosions and fire embers (Zohdi [64-67]), we make the following assumptions:In order to gain insight, initially, we will discuss the closely related, analytically tractable, Stokesian model next.

We assume the same initial velocity magnitude for all particles under consideration, with a random distribution of outward directions away from the source of the cough. This implies that a particle non-interaction approximation is appropriate. Thus, the inter-particle collisions are negligible. This has been repeatedly verified by “brute-force” collision calculations using formulations found in Zohdi [60-63].

We assume that the particles are spherical with a random distribution of radii , . The masses are given by , where is the density of the particles.

We assume that the cough particles are quite small and that the amount of rotation, if any, contributes negligibly to the overall trajectory of the particles. The equation of motion for the ith particle in the system is with initial velocity and initial position . The gravitational force is , where .

For the drag, we will employ a general phenomenological model where is the drag coefficient, is the reference area, which for a sphere is , is the density of the ambient fluid environment and is the velocity of the surrounding medium which, in the case of interest, is air. We will assume that the velocity of the surrounding fluid medium () is given, implicitly assuming that the dynamics of the surrounding medium are unaffected by the particles.1

Remark

As mentioned, there are a large number of physically similar phenomena to a cough, such as the particulate dynamics associated with blasts, explosions and fire embers. We refer the interested reader to the wide array of literature on this topic; see Plimpton [39], Martin-Alberca and Garcia-Ruiz [29], Brock [6], Russell [40], Shimanzu [43], Werrett [55], Kazuma [19, 20], Wingerden et al. [56] and Fernandez-Pello [12], Pleasance and Hart [38], Stokes [45] and Rowntree and Stokes [42], Hadden et al. [14], Urban et al. [48] and Zohdi [67] .

Analytical characterization: simplified Stokesian model

Analysis of particle velocities

For a (low Reynolds number) Stokesian model, the differential equation for each particle is (Fig. 1)where , where is the viscosity of the surrounding fluid (air) and the local Reynolds number for a particle is and is the fluid viscosity. This can be written in normalized form asThis can be solved analytically to yield, for example in the z directionwherewhere the same holds for the y and z directions. The trends areSome special cases:Again, we note that the equations are virtually the same for the x and y directions, with the direction of gravity and fluid flow being the main differentiators.

Fig. 1

Model problem for release of cough particles

,

,

,

,

As

As

The decay rate is controlled by , indicating that small particles attain ambient velocities extremely quickly.

With no gravity:

With no damping:

Analysis of particle positions

From the fundamental equation, relating the position to the velocitywe can write for the z directionwith the same being written for the x and y directions. Integrating and applying the initial conditions yieldsIf and , thenAs As , the travel distance is dramatically shorter. The converse is true, larger particles travel farther.

Settling (airborne) time

The settling, steady-state velocity can be obtained directly fromby setting , one can immediately solve for the steady-state velocityThe trends areIn summary

As , then ,

As , then .

Large particles travel far and settle quickly and

Small particles do not travel far and settle slowly.

Remark

The ratio of the Stokesian drag force to gravity iswhich indicates that for very small particles, drag will dominate the settling process and for larger particles, gravity will dominate.

Computational approaches for more complex models

More detailed characterization of the drag

In order to more accurately model the effects of drag, one can take into account that the empirical drag coefficient varies with Reynolds number. For example, consider the following piecewise relation (Chow [9]):where, as in the previous section, the local Reynolds number for a particle is and is the fluid viscosity.2We note that in the zero Reynolds number limit, the drag is Stokesian. In order to solve the governing equation,we integrate the velocity numericallyThe position is the obtained by integrating again:This approach has been used repeatedly for a variety of physically similar drift-type problems in Zohdi [64-67].

For , ,

For , ,

For , ,

For , ,

For , ,

The piecewise drag law of Chow [9] is a mathematical description for the Reynolds number over a wide range and is a curve-fit of extensive data from Schlichting [44].

Simulation parameters

In order to illustrate the model, the following simulation parameters were chosen:

Starting height of 2 m,

Total simulation duration, 4 s,

The time step size, s,

The cough velocity, (taken from the literature which indicates ),

Density of particles, = 1000 kg/m,

Density of air, , kg/m and

Total mass, kg.

Particle generation

A mean particle radius was chosen to be  m with variations according towhere and a random variable . The algorithm used for particle generation was:

M=0

Start loop:

,

If then stop (determines = particles)

End loop

Initial trajectories

The initial trajectories we determined from the following algorithm

Specify relative direction ‘cone’ parameters: ,

For each particle, , construct a (perturbed) trajectory vector: where , and .

For each particle, normalize the trajectory vector:

For each particle, the velocity vector is constructed by a projection onto the normal vector:

Numerical results

An extremely small (relative to the total simulation time) time-step size of s was used. Further reductions of the time-step size produced no noticeable changes in the results, thus the solutions generated can be considered to have negligible numerical error. The simulations took under 10 s on a standard laptop. The algorithm generated 59,941 particles ranging from m (i.e. ). We used a trajectory cone of and in the example given. Figures 2, 3 illustrate the results for the parameters above (for ). If particles contacted the floor, they were immobilized. The maximum distance travelled from the source located at (0, 0, 2) was 2.72 m (achieved by large particles). Table 1 shows variation in the headwind. For strong tailwind, the larger particles land further away from the cough source. As the analytical theory asserts, successive frames indicate that: (a) Large particles travel far and settle quickly and (b) Small particles do not travel far and settle slowly (when there are no ambient velocities). As observed in the simulations, the settling of the small particles is still not achieved by the end of the simulation time (here 4 s). Accordingly, the simulations were also run for extremely long periods to ascertain that the “mist” of small particles remained airborne for several minutes (as predicted by the theory). For strong opposing headwind, small particles move backwards, and still remain airborne for extended periods of time. This is by far the most dangerous case, since this will encounter other persons at the torso level. We also note that ratio of the general drag to gravity indicates:which indicates that at high velocities, the dynamics are dominated by drag.

Fig. 2

Cough simulation (from a starting height of 2 m, for ): successive frames indicating the spread of particles. a Large particles travel far and settle quickly and b small particles do not travel far and settle slowly

Fig. 3

Zoom on cough simulation (from a starting height of 2 m, for ): successive frames indicating the spread of particles. a Large particles travel far and settle quickly and b small particles do not travel far and settle slowly

Table 1

Maximum distance from the source at the end of T = 4 s

vyf(m/s)	Max-distance (m)	Comments
− 2.0	8.002	Due to small particles moving backwards
− 1.0	4.210	Due to small particles moving backwards
0.0	2.721	Due to large particles moving forwards
1.0	4.937	Due to large particles moving forwards
2.0	8.736	Due to large particles moving forwards

Summary and extensions

For general cough conditions, there can be cases where the change in the surrounding fluid’s behavior, due to the motion of the particles and cough, may be important. The result is a system of coupled equations between the particles and the fluid, requiring spatio-temporal discretization (high-fidelity Finite Elements or Finite Differences) of the classical equations governing the surrounding fluid mechanics (Navier Stokes)where is the density field of the fluid, is the fluid velocity field, is the fluid stress field, is the fluid velocity gradient field, is the body force field, is the fluid pressure field, and are fluid material property fields.3 It is important to emphasize that physically compatible boundary data must be applied, and this is not a trivial matter for compressible flow. Additionally, the first law of thermodynamics should be included (along with equations for various chemical reactions), which reads aswhere w(x) is the stored energy in the fluid, is the heat flux field, z is the heat source field per unit mass. Generally such models are ineffective for rapid real-time use, but are quite useful for detailed offline background analyses, where a rapid response is a nonissue. The continuum discretization is usually combined with a Discrete Element Method for the particle dynamics. There are a variety of such approaches, for example, see Avci and Wriggers [2], Onate et al. [34, 35], Leonardi et al. [23], Onate et al. [36], Bolintineanu et al. [4] and Zohdi [60, 63]. Such models are significantly more complex than the models used in the current paper. More detailed analyses of fluid-particle interaction can be achieved in a direct, brute-force, numerical schemes, treating the particles as part of the fluid continuum (as another fluid or solid phase), and thus meshing them in a detailed manner. In such an approach (for example see Avci and Wriggers [2])Along these lines, in Zohdi [60, 63], more detailed, computationally intensive models were developed to characterize the motion of small-scale particles embedded in a flowing fluid where the dynamics of the particles affects the dynamics of the fluid and vice-versa. In such a framework, a fully implicit Finite-Difference discretization of the Navier-Stokes equations was used for the fluid and a direct particle-dynamics discretization is performed for the particles. Because of the large computational difficulty and expense of a conforming spatial discretization needed for large numbers of embedded particles, simplifying assumptions are made for the coupling, based on semi-analytical computation of drag-coefficients, which allows for the use of coarser meshes. Even after these simplifications, the particle-fluid system was strongly-coupled. The approach taken in that work was to construct a sub-model for each primary physical process. In order to resolve the coupling, a recursive staggering scheme was constructed, which was built on works found in Zohdi [60-63]. The procedure was as follows (at a given time increment): (1) each submodel equation (fluid or particle-system) is solved individually, “freezing” the other (coupled) fields in the system, allowing only the primary field to be active, (2) after the solution of each submodel, the associated field variable was updated, and the next submodel was solved and (3) the process is then repeated, until convergence. The time-steps were adjusted to control the rates of convergence, which is dictated by changes in the overall physics. Specifically, the approach was a staggered implicit time-stepping scheme, with an internal recursion that automatically adapted the time-step sizes to control the rates of convergence within a time-step. If the process did not converge (below an error tolerance) within a preset number of iterations, the time-step was adapted (reduced) by utilizing an estimate of the spectral radius of the coupled system. The developed approach can be incorporated within any standard computational fluid mechanics code based on finite difference, finite element, finite volume or discrete/particle element discretization (see Labra and Onate [22], Onate et al. [34, 35], Rojek et al. [41] and Avci and Wriggers [2]). However, while useful in many industrial applications where high precision is required, the use of such a model for the coarser applications of interest in this work is probably unwarranted.

A fluid-only problem is solved, with (instantaneous) boundary conditions of at each point on the fluid-particle boundaries, where the velocity of the points on the boundary are given by where is the center of mass and where is the particle angular velocity for each of the individual particles and is a vector from the mass center to the surface.

For each particle, one would solve: and where the forces and moments would have a contribution from the fluid drag (with particle occupying domain and outward surface normal ) is defined as and

At a time-step, the process is iteratively driven by solving the fluid-only problem first, then the particles-only problem, and repeated until convergence in an appropriate norm.

In closing, we remark on a closely related theme to the one described in this paper, namely decontamination. For example, decontamination based on UV technology has become ubiquitous, with many variants now being proposed. UV light varies in wavelength from 10 to 400 nm, thus making it shorter that visible wavelengths and larger than x-rays. Short wave UV light (UV-c) can damage DNA and sterilize surfaces making in useful in the medical industry. This was first noted in 1878 (Downes and Blunt [10]) when the effect of short-wavelength light killing bacteria was discovered. By 1903 it was known the most effective wavelengths were around 250 nm (UV-c), for which Niels Finsen won a Nobel Prize (for skin-based tuberculosis eradication using UV light). Contaminants in the indoor environment are almost entirely organic carbon-based compounds, which break down when exposed to high-intensity UV at 240–280 nm. Despite the attractiveness of using UV-c light, the literature has shown that it is difficult to ensure that all surfaces are completely decontaminated due to shadowing effects. Thus, the use of ultraviolet germicidal irradiation (UVGI) is effective only as a component in a multistage process-it alone carries the risk of residual contamination. Thus, purely UV-c protocols should be adopted if there is no other choice. However, they can be an integral part of a multistage process involving a combination of (a) gas vapors and (b) heat and humidity. The topic of decontamination technologies is of paramount interest (see references Anderson et al. [1], Battelle [3], Boyce et al. [5], Card et al. [7], Heimbuch and Harish [15], Heimbuch et al. [16], Ito and Ito [17], Lin et al. [24], Kanemitsu [18], Lindsley et al. [25], Lore et al. [27], Marra et al. [28], Mills et al. [30], Tseng and Li [47], Viscusi et al. [51] and Nerandzic et al. [32]), and the corresponding simulation of such processes has recently been undertaken in Zohdi [68] and is a topic of ongoing research.