Computational investigation of prolonged airborne dispersion of novel coronavirus-laden droplets.

We have performed highly accurate numerical simulations to investigate prolonged dispersion of novel coronavirus-laden droplets in classroom air. Approximately 10,900 virus-laden droplets were released into the air by a teacher coughing and tracked for 90 min by numerical simulations. The teacher was standing in front of multiple students in a classroom. To estimate viral transmission to the students, we considered the features of the novel coronavirus, such as the virus half-life. The simulation results revealed that there was a high risk of prolonged airborne transmission of virus-laden droplets when the outlet flow of the classroom ventilation was low (i.e., 4.3 and 8.6 cm/s). The rates of remaining airborne virus-laden droplets produced by the teacher coughing were 40% and 15% after 45 and 90 min, respectively. The results revealed that students can avoid exposure to the virus-laden droplets by keeping a large distance from the teacher (5.5 m), which is more than two times farther than the currently suggested social distancing rules. The results of this study provide guidelines to set a new protection plan in the classroom to prevent airborne transmission of virus-laden droplets to students.

Introduction

At present, life shows no sign of returning to normal from the novel coronavirus (COVID-19) pandemic, and the number of new daily infections is about to reach its peak around the world. The features of infection by COVID-19 are different from those of previous coronaviruses and the influenza virus, so infection routes other than the contract transmission route have attracted the attention of researchers. Bahl et al. (2020) observed the virus-laden droplets expelled during coughing and sneezing with a high-speed camera. Asadi et al. (2019, 2020a, 2020b) measured the droplet emission rate during human speech and discussed the diffusion risk. They demonstrated that the virus-laden droplets expelled during human speech can significantly disperse, and they discussed the possibility of infection spread due to asymptomatic patients. The basic mechanisms of inhalation of airborne particles and particle deposition in the respiratory airway were summarized in a review article by Inthavong (2020).

Regarding the transmission route of COVID-19, the possibility of airborne transmission as a major route of infection has recently been in the limelight. In this route, virus-laden droplets expelled from the mouth and nose during coughing, sneezing, and talking can remain in the air for a long time and may travel a long distance. Numerical simulations of the droplet and airborne transmission routes have been reported in multiple computational modeling studies. Thatiparti, Ghia, and Mead (2016) performed numerical simulations of cough aerosol dispersal in a mock airborne infection isolation room. Yang et al. (2018) investigated the effect of a cough jet on the airflow field and transport characteristics of droplets in an airliner cabin section. Yu et al. (2018) simulated bioaerosol emission in an office with mechanical ventilation and air-conditioning systems. Dbouk and Drikakis (2020) investigated the effect of the wind speed on the travel distance of saliva droplets generated by human coughing. Although these studies have provided useful information about droplet and airborne transmission, the simulation times used in the studies were not long enough to accurately evaluate virus-laden droplet dispersion (e.g., at most 20 s). There was also a lack of consideration of the novel coronavirus characteristics, such as the virus life-time. In March and April 2020, after the COVID-19 pandemic started, our research group suggested the possibility of long-distance travel of the coronavirus in several media, such as The New York Times (2020) and .

The aim of this study is to investigate prolonged airborne dispersion of novel coronavirus-laden droplets in detail using numerical simulations. Approximately 10,900 virus-laden droplets were released into the air in a classroom with students by a teacher's cough. While tracking the droplets for 90 min, we took into account the half-life of the novel coronavirus (Doremalen et al., 2020; Stadnytskyi et al., 2020) to assess the risk of virus transmission with high accuracy by numerical simulation.

Numerical approach

Governing equations

In this study, the airflow in the room was first simulated, and then the translational motion of the virus-laden droplets was tracked using the computed airflow field. The three-dimensional incompressible Navier−Stokes equations are given bywhere represents the non-dimensional variable. u a, v a, and w a are the air velocity components in the x, y, and z directions, respectively. Each axis is explained in Section 2.4. p is the pressure and Re is the Reynolds number. The continuity equation is defined by

The variables in Eq. (1)−(3) are non-dimensionalized as follows: where L 0 is a characteristic length, U 0 is a characteristic velocity, ρ a is the density of the air (1.2 kg/m3), and ν a is the kinematic viscosity of the air. Both the characteristic length and velocity are explained in Section 2.4.

Numerical schemes

Unstructured grids with high shape-reproducibility were used to calculate the airflow field in the body-fitted coordinate system. The governing equations (Eq. (1)) were discretized using the unstructured finite-volume method in which the physical quantities were defined at the center of each cell. Integrating Eq. (1) over the control volume Ω using the Gauss theorem giveswhere l (l = 1, 2, 3, 4) is the number of each boundary surface of the control volume, and n x, n y, and n z are the components of the outward normal vector in the x, y and z directions, respectively. Using the fractional step method, Eq. (6) was divided into two equations for the velocity and pressure, which were iteratively solved using the lower-upper Symmetric-Gauss-Seidel (LU-SGS) method (Yoon & Jameson, 1988) and the biconjugate gradient stabilized method (Bi-CGSTAB) method (Vorst, 1992), respectively. The convergence values for both the velocity and pressure were set to less than 0.1% of their initial residuals. The simulations were performed using an in-house computational fluid dynamics code.

Virus-laden airborne droplet motion

In this study, because we mostly focused on very small virus-laden airborne droplets, we ignored the influence of the droplet motion on the flow field. We simulated the traces of the airborne droplets in the precalculated airflow field, as described in Section 2.1. Considering that droplets expelled from the mouth are composed of saliva, the droplets in this study were water. Additionally, we assumed that all droplets were spherical and the interactive forces the among droplets were ignored.

The droplet translational motion is defined aswhere

ρ d, V, and r are the density (998 kg/m3), volume, and radius of an individual droplet, respectively. u d, v d, and w d are the droplet velocity components in the x, y, and z directions, and u ac, v ac, and w ac are the air velocity components in the x, y, and z directions defined at a center of each droplet. g is gravitational acceleration. Here, the air−water interface has a no-slip condition owing to contamination in saliva, such as protein (Takagi & Matsumoto, 2011). For this reason, we applied the following equation for the solid sphere (Flemmer & Banks, 1986) to the drag coefficient CD:

Red is the droplet Reynolds number.

The change in the droplet radius due to evaporation was estimated by the Clausius−Clapeyron relation (Wallace & Hobbs, 2006):where t is the time, D is the diffusion coefficient of water vapor, RH is the relative humidity, e s is the saturation vapor pressure, R v is the gas constant for water vapor, T is the temperature (20 °C), T 0 is the base temperature (0 °C), and L is the latent heat of vaporization of water. At T 0 = 0 °C, e s (0) = 6.11 × 102 Pa. Here, the droplets are assumed to remain in thermal equilibrium with their surroundings and do not cool owing to evaporation. The ratio of the equilibrium diameter to the initial diameter as a function of the humidity for the droplet was reported by Marr et al. (2019). In this study, the minimum radii of the individual droplets, which are no longer affected by evaporation, were calculated based on this ratio.

The half-life of the COVID-19 virus (i.e., 1.1 h) was introduced into the simulations following Doremalen et al. (2020). This half-life plays a large role in risk assessment because the virus-laden droplets were tracked for 90 min in the numerical simulations. Based on the half-life of the virus, we evaluated whether the virus was active using random numbers at 1 s intervals, and droplets with inactive viruses were accordingly removed to assess the risk of active virus transmission.

In terms of droplet coalescence, the radius and position of the droplet after coalescence of droplets were determined based on the following equations: where the subscripts 1 and 2 represent the two droplets before droplet coalescence, while the subscript 3 represents the droplet after coalescence. h is the position vector of the droplet. Note that under our simulation conditions, droplet coalescence very rarely occurred. For this reason, we estimated the infection risk based on the number of droplets instead of the mass of the droplets (see Section 3).

Computational domain and modeling conditions

We simulated coronavirus-laden airborne droplet transmission in a school classroom. Regarding previous studies on airborne transmission in a school environment, Abuhegazy et al. (2020) carried out fluid-particle dynamics simulations with the steady-state Reynolds average Navier−Stokes incompressible solver to investigate aerosol transport and surface deposition in a realistic classroom environment. Shao et al. (2021) performed quantitative risk assessment of airborne transmission of COVID-19 under a small classroom setting by combining novel in situ measurements and computational fluid dynamics simulations with the Eulerian−Lagrangian framework. In contrast to previous studies, we incorporated the half-life of the novel coronavirus into the simulations and considered longer simulation time than the above studies, as described in Section 2.3.

The classroom model with a teacher and nine students in the classroom is shown in Fig. 1 . The classroom dimensions in the x, y, and z directions (i.e., L x × L y × L z) are 7.5 m × 8.0 m × 3.0 m, and there are four air inlet vents (0.55 m × 0.25 m) and an outlet vent (0.35 m × 0.35 m). The classroom model geometry with the student numbers is shown in Fig. 2 . The teacher is standing in front of the class and facing the students. We assumed that the teacher is infected by coronavirus. The computational grids of the model are shown in Fig. 3 . The modeled teacher and students are configured using unstructured grids. The total number of elements in the model is 5,984,579. Here, a grid independence test was performed using three meshes with 2,645,201 (coarse), 5,984,579 (medium; used in the final simulations), and 10,778,504 (fine) cells. The average displacement of the airborne droplets (sizes of 0.9, 1.5, 3.0, 5.0, and 7.0 μm) between the final and fine meshes was less than 0.01%, indicating that the results were grid independent.

Fig. 1

Illustration of the classroom model with a teacher and students.

Fig. 2

Illustration of the classroom model geometry with the student numbers.

Fig. 3

Computational model grids. (a) Boundary surfaces and (b) cross-sectional view around the teacher and a student.

As an initial condition, the air velocity and pressure (u a, v a, w a, and p at all of the cells) were set to zero. As a boundary condition, the Neumann condition was used for the flow at the inlet vents, and p was set to zero. Two flow speeds at the outlet vent were used (4.3 and 8.6 cm/s), while the Neumann condition was used for the pressure. The non-slip condition for the velocity and the Neumann condition for the pressure were imposed on all of the walls. The characteristic velocity U 0 and characteristic length L 0 in Eqs. (4), (5) were set to the flow speed at the outlet vent and its height, respectively, so that the Reynolds numbers based on U 0 and L 0 were Re = 1000 and 2000. Additionally, we used three values for the relative humidity (RH = 30%, 60%, and 90%). We performed simulations for four cases (Re = 1000 and RH = 30%, Re = 1000 and RH = 60%, Re = 1000 and RH = 90%, and Re = 2000 and RH = 60%). Considering the low Reynolds number conditions, the use of fine meshes, and the relationship between the time required for a cough (i.e., 0.5 s) and the whole simulation time (i.e., 90 min), we applied the laminar flow model (Kim et al., 2018; Zhang & Kleinstreuer, 2002) to the simulations. Additionally, we assumed that the flow field was isothermal.

We applied different breathing flow rates for the teacher and students (Fig. 4 ) (Fenn & Rahn, 1965). In the middle of breathing, the teacher coughed once for approximately 0.5 s. The coughing flow rate of the teacher used in the model is shown in Fig. 5 (Gupta et al., 2009). The non-dimensional time step was 0.0001. At Re = 2000, a cough and a breathing cycle took 1119 time steps and 9874 time steps, respectively. Before the cough, the indoor airflow, whose behavior is only determined by the vents, was computed for 200 s. During this period, the teacher and students did not breathe. After this initial 200 s period, the students began breathing at different phases of the breathing cycle, and the teacher coughed. Following the cough, the teacher resumed breathing. The distribution of the diameters of the droplets (d) expelled from the teacher's mouth during a cough is shown in Fig. 6 . This was determined by considering several experimental results (Duguid, 1946; Yang et al., 2007). The total number of droplets produced by the teacher's cough (N 0) was 10,908.

Fig. 4

Breathing flow rates for the teacher and students.

Fig. 5

Coughing flow rate for the teacher.

Fig. 6

Distribution of droplet diameters expelled from the teacher's mouth for a cough.

Results and discussion

Distribution of the virus-laden airborne droplets

The distributions of the virus-laden airborne droplets in the classroom at t = 5, 15, 30, 60 and 90 min after the teacher's cough are shown in Fig. 7 (Re = 1000 and RH = 60%). The red dots in the figure indicate the droplets expelled from the teacher's mouth. Large droplets settle close to the teacher (around the teacher's legs and on the top of their table). Smaller droplets travel farther towards the area where the students are sitting (hereafter called the student area). At t = 5 min, a large number of airborne droplets are around the center of the first row in the student area. The airflow in this region is dominated by both the cough and the subsequent normal breathing of the teacher and students because the region is far from the outlet vent and the outlet flow is slow (i.e., 4.3 cm/s). From t = 15–60 min, the droplets disperse over a relatively wide range of the student area. At t = 90 min, the droplets finally tend to travel toward the outlet vent because a number of the droplets seem to distribute at the lower right side of the figure indicated by the dashed blue circle. As described in Section 2.3, this simulation takes the half-life of the novel coronavirus (i.e., 1.1 h, Doremalen et al., 2020) into account, and, hence, inactive viruses are removed. Owing to the effect of the inactivation and expansion of the spatial dispersion of the airborne droplets, the local concentration of droplets decreases with time.

Fig. 7

Distributions of virus-laden airborne droplets in the classroom at (a) t = 5 min, (b) t = 15 min, (c) t = 30 min, (d) t = 60 min and (e) t = 90 min after a cough (Re = 1000 and RH = 60%).

Local airborne droplet concentration at different heights

To evaluate the temporal change in the spatial distribution of the virus-laden airborne droplets in detail, the local airborne droplet concentrations (α) at different heights (locations in the z direction) for Re = 1000 and RH = 60% are shown in Fig. 8 . Here, α is defined as the ratio of the number of airborne droplets in the unit volume of 0.21 m × 0.22 m × 0.5 m (in the x, y, and z directions, respectively) to the total number of droplets (N 0). At t = 5 and 15 min, the droplet concentrations are relatively high at z/L z = 0.25 and 0.42. This trend is maintained even at t = 30 min. However, at t = 60 min, regions with relatively high droplet concentrations appear at z/L z = 0.75. Note that the red point at z/L z = 0.42 can be regarded as a singular point because the droplets accidently accumulate at a stagnation point in the flow field. Finally, at t = 90 min, most of droplets disappear at z/L z = 0.25, 0.42, and 0.58, except for the singular point, while droplets still remain at z/L z = 0.75. However, this region is sufficiently far from the students, so the droplets in the region would barely infect the students.

Fig. 8

Local airborne droplet concentrations α at (a) t = 5 min, (b) t = 15 min, (c) t = 30 min, (d) t = 60 min and (e) t = 90 min when Re = 1000 and RH = 60%.

Next, we evaluate the effect of the RH and Reynolds number on the local airborne droplet concentration (Fig. 9 ). At several z/L z locations, the airborne droplet concentration averaged over the x-y plane is included in the figure. In terms of the RH at the constant Reynolds number (i.e., Re = 1000), the dependency of the RH on the droplet concentration is not very clear. However, a change in the droplet concentration with Re is observed at constant humidity (RH = 60%). When t = 5 min, at z/L z = 0.25, the droplet concentration is higher at Re = 1000 than at Re = 2000, while it is the opposite at z/L z = 0.42. A similar trend is observed at t = 15 min. At t = 30 min, there is a large difference in the concentrations between Re = 1000 and 2000 at z/L z = 0.58, and the concentration at Re = 2000 is higher than that at Re = 1000. That is, at Re = 2000, more virus-laden droplets remain in the air 30 min after the cough compared with at Re = 1000.

Fig. 9

Local airborne droplet concentrations α at different Reynolds numbers and humidity levels. (a) t = 5 min, (b) t = 15 min and (c) t = 30 min.

Time evolution of the number and size distribution of the airborne droplets

The time evolution of the airborne droplet number (N) normalized by the total number of droplets (N 0) at different Reynolds numbers and humidity levels is shown in Fig. 10 . In all cases, around 40% of the initial droplets are airborne 45 min after the cough, and approximately 15% of the droplets remain in the air 90 min after the cough. Under the constant Reynolds number condition (i.e., Re = 1000), there is little difference in the ratios of airborne droplets (N/N 0) between RH = 30% and 60% throughout the 90-min duration. In addition, at RH = 90%, N/N 0 is lower than the N/N 0 values at RH = 30% and 60%. Because of less evaporation, the droplets at RH = 90% are relatively heavy compared with those at RH = 30% and 60%. Such droplets easily settle, so the number of airborne droplets decreases. Moreover, under the constant humidity condition (i.e., RH = 60%), the difference in the ratios of airborne droplets between Re = 1000 and 2000 is negligible until around t = 30 min. After this time, the difference gradually increases, and at t = 90 min, the airborne droplet ratio is much higher at Re = 2000 than at Re = 1000. This is contrary to our initial expectation because the flow speed at the outlet is higher at Re = 2000 than at Re = 1000. However, this can be explained as follows. In the computation domain, there is an upward flow from the inlets towards the outlet. This flow prevents very small airborne droplets from settling, and it becomes slightly stronger at Re = 2000 (Fig. 11 ). As a consequence, more droplets can remain in the air at Re = 2000 than at Re = 1000. Note that this trend would change at higher Reynolds numbers (e.g., Re > 5000), because under such Reynolds numbers, an induced strong flow can sweep more airborne droplets towards the outside of the room. We therefore speculate that the trend described here is restricted to the case of low Reynolds numbers.

Fig. 10

Time evolution of the airborne droplet number N normalized by initial number of droplets N0 at different Reynolds numbers and relative humidity levels.

Fig. 11

Distributions of air velocity in the z direction wa at t = 90 min. (a) Re = 1000, RH = 60% and (b) Re = 2000, RH = 60%.

The size distributions of the airborne droplets at t = 5, 30 and 90 min are shown in Fig. 12 . For reference, the size distribution of the droplets produced by the cough is included in the figure. However, the result at RH = 30% is not included in this figure because there is little difference in the N/N 0 values between 30% and 60%, as shown in Fig. 10. At t = 5 min, in all cases (excluding the initial condition), the distribution shifts towards smaller diameter and there are no droplets with diameters larger than 22.5 μm. This is mainly due to evaporation from the droplet interface of smaller droplets, as well as larger droplets falling to the floor. In addition, the distribution is almost independent of Re and RH. At t = 30 min, the droplets with diameters larger than 10 μm settle and no airborne droplets are observed. Moreover, the difference in the distribution between Re = 1000 and 2000 is very small, while there is only a significant difference between RH = 60% and 90% in the range 1.5 μm ≤ d ≤ 3 μm. In all cases, at t = 90 min, the distribution further shifts towards smaller diameter and there are no airborne droplets larger than 5 μm.

Fig. 12

Size distributions of airborne droplets at (a) t = 5 min, (b) t = 30 min and (c) t = 90 min.

Contact of virus-laden droplets with the students

To clarify contact of virus-laden droplets with the students in the room, the time evolution of the cumulative number of droplets that reach each student is shown in Fig. 13 (Re = 1000 and RH = 60%). When an airborne droplet exists in a cell within the student's body surface or mouth, it is considered “in contact with” or “attached to” a student. Illustration of the classroom model geometry with the student numbers is shown in Fig. 3. Among the students, student 4, who sits in front of the teacher, is exposed to the most virus-laden droplets, and exposure lasts for approximately 60 min. In particular, approximately 30 droplets come into contact with the student within the first 3 min after the cough, which results not from airborne transmission, but from direct droplet transmission from the cough. Virus-laden airborne droplets that remained above the students gradually start to attach to the other students at around 20 min after the cough. This is induced by dispersion of virus-laden airborne droplets with diameters of less than 20 μm. Students 7 and 1, who sit in the first row in the student area, are exposed to virus-laden droplets until around t = 50 and 80 min, respectively. Students 2, 5, and 8, who sit in the second row in the student area, are still exposed to droplets even at t = 90 min. However, there are no virus-laden droplets attached to students 3, 6 and 9, who sit in the third row in the student area. Student 7 stops being exposed to new virus-laden droplets at around t = 50 min, which is the shortest time among the students sitting in the first and second rows in the student area. This is mainly owing to the upward flow, because student 7 is the closest to the outlet.

Fig. 13

Time evolution of the cumulative number of droplets that reach each student n at Re = 1000 and RH = 60%.

The time evolution of the cumulative number of droplets that reach each student at different Reynolds numbers and humidity values is shown in Fig. 14 . Regardless of the Re and RH levels, most of the students are exposed to virus-laden droplets for more than 30 min, except for the students sitting in the last row from the teacher. In addition, under the constant Reynolds number condition (i.e., Re = 1000), more exposure to droplets occurs at RH = 90% than at RH = 60%. Moreover, under the constant humidity condition (i.e., RH = 60%), the students are exposed to fewer virus-laden droplets at higher Re. In particular, a significant decrease in droplet attachment occurs for students 1 and 2. This is related to the stronger upward flow from the inlets towards the outlet induced at Re = 2000, and droplet attachment to students 1 and 2 is more easily affected by the induced flow because they sit around the edge of the airborne droplet propagation region, as shown in Fig. 7. By contrast, in all cases, students 3, 6, and 9, who sit in the third row, have no exposure to the virus-laden droplets for the whole 90-min duration. The distance between these students and the teacher is at least 5.5 m, which is more than twice the recommended social distance (i.e., 2 m) used for prevention of droplet transmission. This result will be useful as a preventative measure when considering airborne transmission in a classroom where the airflow speed at the outlet vent is low, like the simulations.

Fig. 14

Time evolution of the cumulative droplet number that reach each student n at different Reynolds numbers and relative humidity levels. (a) Re = 1000, RH = 30%, (b) Re = 1000, RH = 60%, (c) Re = 1000, RH = 90% and (d) Re = 2000, RH = 60%.

Summary and conclusions

We have investigated airborne dispersion of novel coronavirus-laden droplets by numerical simulations. In the simulations, the infected teacher stood facing nine students sitting in a classroom with four inlet vents and one outlet vent. Where the airflow speed at the outlet vent was low, the virus-laden droplets expelled by the teacher's cough dispersed over a relatively wide range of the student area. The rates of remaining airborne virus-laden droplets produced by the teacher's cough were 40% and 15% after 45 and 90 min, respectively. This indicates that in situations without appropriate air ventilation, there is a risk of prolonged airborne transmission of virus-laden droplets in the room. However, students can avoid exposure to the virus-laden droplets by staying a large distance from the teacher (5.5 m), which is more than two times farther than the suggested social distancing rules. The results of this study provide guidelines to set a new protection plan in the classroom to prevent airborne transmission of novel coronavirus-laden droplets to students.

Declaration of competing interest

There is no conflict of interest in this paper.