Direct numerical simulation of turbulent dispersion of evaporative aerosol clouds produced by an intense expiratory event

Responsible for more than 1.5 × 10⁶ deaths worldwide during 2020, the recent COVID-19 pandemic has made evident the need to better understand how airborne infectious diseases spread among individuals. It is an accepted scientific fact that flow physics plays a central role in the transmission of the SARS-CoV-2 virus (Severe Acute Respiratory Syndrome CoronaVirus), responsible for COVID-19, and other airborne pathogens.¹ On one hand, the turbulent fluid set in motion by exhalation is responsible for the aerosolization of ejected liquid spewed during a respiratory event.^2–4 This turbulent flow also controls the dispersion of the associated virus-laden particle cloud in the environment, especially dispersion close to the source of the aerosols. On the other hand, ambient currents produced by ventilation and air conditioning systems dominate the long-term transport⁵ and the surface deposition of the smallest droplets. Finally, occurrence of new infections is also affected by the process of inhaling suspended aerosols in air and the deposition process in the human upper airway.⁶ Thus, understanding how particles are transported by the action of moving fluids is key in designing and implementing preventative strategies including the use of masks,^7–9 indoor ventilation,^10,11 and social distancing measures.¹² Numerical simulations of specific indoor and public transportation environments such as isolation rooms,¹³ airliner cabin sections,¹⁴ urban buses,¹⁵ elevators,¹⁶ parking areas,¹⁷ restaurants,¹⁸ and offices¹⁹ have been conducted in order to design robust transmission mitigation strategies.

Identified as a major pathway for SARS-CoV-2 transmission, the airborne route occurs when pathogen-laden particles are expelled during an expiratory event including talking, singing, coughing, or sneezing. The breakup of the spewed material due to multiple instabilities²⁰ generates droplets and aerosols of different sizes.^21–24 The diameter of the particle and the characteristics of background flow produced by the exhalation are the two most relevant factors controlling the initial residence time and the range of the particle-laden cloud of infectious particles. Historically, particles larger than 100 μm (Refs. 25 and 26) in diameter have been classified as droplets with gravity-dominated dynamics characterized by ballistic trajectories. While these relatively large particles may contain significant amounts of virus copies, their horizontal range and residence time afloat are very limited.²¹ In contrast, the dynamics of particles of 100 μm in diameter or less is dominated by the hydrodynamic drag exerted by the background flow that allows them to remain suspended in the air longer and travel greater distances.

These distinctive transport properties according to each particle size are complicated when the ambient conditions (temperature and humidity) favor the evaporation of the aqueous phase in droplets and aerosols.^7,27 Most droplets are suspected to evaporate in a few seconds after being released,²⁸ forming droplet nuclei consisting mostly of virions, solid residue,^29,30 and some nonevaporable residual.³¹ With a typical size on the order of 1 μm, each nucleus may potentially remain afloat for hours while carrying multiple virions with an approximate half-life of 1 h.³² When a single virion is suspected to potentially cause a new infection,³³ the importance of this type of droplet nucleus in the transmission of COVID-19 is key. As particles shrink due to evaporative effects, the dominant contribution to their dynamics may transition from gravity to hydrodynamic drag, affecting the particle cloud dispersion.

Using both experiments^34,35 and numerical simulations,^36,37 recent studies analyzed several key aspects of the hydrodynamics produced by different expiratory events, including sneezing,³⁴ coughing,³⁶ talking, and breathing,³⁵ and the dispersion processes of the resulting aerosol cloud.

On the numerical front, most studies have relied on some form of the Reynolds Averaged Navier Stokes (RANS) approach to incorporate the impact of turbulence to the overall transport of momentum and scalars, mainly temperature and relative humidity.

Thus, Dbouk and Drikakis used the RANS approach to investigate the effect of background wind currents³⁶ and face masks³⁸ on the aerosol dispersion of droplets and aerosols produced by expiatory events. Using the RANS methodology for compressible multiphase mixtures and the $κ−ω$ turbulence model for the carrier phase along with a Lagrangian approach for the saliva droplets and aerosols, Dbouk and Drikakis numerically simulated the jet produced by a cough.³⁶ They found that an ambient wind speed of 15 km h⁻¹ can transport airborne particles to distances of up to 6 m from the source.

The numerical work by Busco et al.³⁴ used the RANS approach to study the jet flow produced by a sneeze assuming the continuous phase to be a mixture of air and water vapor. Equipped with a $κ−ϵ$ model to account for the turbulent transport contribution and assuming the saliva droplets to be Lagrangian particles, they demonstrated the impact of the exhalation release angle on the airborne droplet cloud dispersion.

The impact of human physiological factors (e.g., illness, stress condition, and anatomy) was addressed by Fontes et al. who used numerical simulations to simulate a sneeze.³⁷ These authors modeled the carrier phase under the Eulerian framework using Detached Eddy Simulations (DES) with a hybrid solver that combined unsteady Reynolds Averaged Navier–Stokes (URANS) for the flow within the boundary layer and Large Eddy Simulations (LES) in the outer region. Again, droplets were assumed to be Lagrangian particles. Their results suggested that changes in the nasal and buccal passages may have a significant impact on the horizontal range traveled by the ejected droplets with nasal obstruction, potentially leading to increases of up to 60%. In addition, physico-chemical characteristics of the saliva were found to affect the overall dispersion properties of the spray generated by a sneeze.

Using LES, Pendar and Pascoa studied the dispersion of saliva droplets produced by violent air exhalation events³⁹ such as coughing and sneezing. Similar to Fontes et al., they used an Eulerian approach for the carrier phase assuming respiratory droplets to obey the Lagrangian dynamics. Assuming a compressible flow, the hydrodynamics was modeled using LES with a subgrid-scale stress model based on a one equation eddy-viscosity model. The authors found that for a strong sneeze, large droplets of 540 μm in diameter could travel up to 4 m with droplets taking up to 3 s to settle down in the absence of background flow. They observed that bending of the head during sneezing could reduce the droplet travel distance by 22% and that wearing a face mask could further reduce the droplet travel distance to 0.6 m. In addition, they also concluded that background flow turbulence must play a key role in the flow physics of jets produced by expiratory events, stressing the need for accurate modeling of the fluctuating transport contribution to the flow dynamics.

Alternative approaches to turbulent modeling based on random walk methodology were used by Wang et al. to determine the effect of several environmental factors on droplet transport and dispersion.⁴⁰ Renzi and Clarke showed that the dynamics of the jet produced by expiratory events could be modeled by extending the theory of buoyant vortex rings.⁴¹ By coupling the integral model for a continuous phase and the Lagrangian particle-tracking approach, these authors explored the effect of initial conditions on the horizontal range of the droplet cloud. The reported key role of the vortex structure on keeping droplets afloat and increasing the residence time of airborne particles provides additional motivation to study in detail the turbulent structure of the flow produced by the rapid air exhalations produced by coughs.

In order to gain insight into the role played by the particle size and the evaporation on the dispersion of particle clouds produced by a mild violent expiratory event, we numerically simulated the trajectories of thousands of particles evolving in the flow produced by an idealized cough solved using Direct Numerical Simulations (DNS). This simulated cough can be divided into two stages: during the exhalation that lasts for 0.4 s, air generates a warm and moist jet that horizontally penetrates into an initially quiescent environment at lower temperature and water vapor content. Once the exhalation ceases, the jet evolves into a thermal puff that continues to grow as it bends in the vertical direction by the action of buoyancy and dissipates due to viscous effects. By considering both evaporable and nonevaporable particles with seven different initial diameters, ranging from 4 to 256 μm, we could numerically investigate with an unprecedented detail the evolution of respiratory particles traveling in the air under typical temperature and relative humidity conditions.

From the point of view of fluid dynamics, respiratory events generate multiphase jets with more or less rich dynamics. The specific length and spatial scales of the problem determine the flow regime and the distribution of exhalation velocity and particle sizes at release.⁴² The flow realization presented here lasts from the onset of the cough at t = 0 until t = 1.65 s. This time span covers the initial dispersion regime over which particle dispersion is dominated by the hydrodynamic drag of the jet-to-plume flow up to approximately $t≈0.75$ s. At that point, the average speed of the particle cloud has dropped to typical indoor air speed levels at which time ambient air currents would take over the dispersion of the aerosols. By assuming the reasonable values of the viral concentration in the spewed material, the aerosol concentration field obtained from the fully resolved simulation allows us to estimate the spatial distribution of the COVID-19 infection risk at the last stages of a mild cough under typical ambient conditions.

In the following sections, first, in Sec. II, the numerical model has been described in detail, covering both governing equations and the pertinent boundary conditions for both the continuous and dispersed phases. In Sec. III, the results have been discussed, illustrating the effect of evaporation on aerosol dispersion. Finally, the paper concludes with summarizing the main findings and listing future research directions.

The flow produced by a violent expiratory event was investigated by Fabregat et al.⁴³ who used DNS to numerically simulate an idealized, isolated mild cough occurring in an isothermal and stagnant environment. In the flow setup, exhaled air, with the temperature and relative humidity of $T0=34$ °C and $ξ0=85$%, was injected through a cylindrical pipe of diameter d representing the mouth into a quiescent environment at $T∞=15$ °C and $ξ∞=65$%. To mimic the transient nature of the violent exhalation, the inlet velocity w₀ was assumed to linearly increase from $w0(t=0)=0$ to the peak value $wm=4.8$ m s⁻¹ at $tm=0.15$ s to then linearly decayed back to zero at $tc=0.4$ s, i.e.,

This velocity sequence reproduces the overall characteristics of a cough in agreement with the measurements reported by Gupta et al.⁴² As the exhaled moist air accelerated/decelerated according to Eq. (1), it produced a jet that penetrated into the unperturbed ambient air. Once the cough ceased for $t>tc$⁠, the jet evolved into a buoyant puff characterized by the deflection of the front trajectory due to buoyancy and the decay of turbulence due to dissipative effects.

Assuming flow incompressibility and using the inlet diameter d, the peak cough velocity w_m and the temperature difference $ΔT=T0−T∞$ as length, velocity, and temperature scales, the nondimensional mass, momentum, and energy conservation equations can be written as

where $t̃=twm/d$ is the time, $p̃=p/ρwm2$ is the pressure, $θ̃=(T−T∞)/ΔT$ is the nondimensional temperature perturbation, δ_ij is the Kronecker delta, and $ũi=(ũ,ṽ,w̃)$ is the velocity field with coordinates $x̃i=(x̃,ỹ,z̃)$ in the spanwise, gravity-aligned, and streamwise directions. Note that the tilde symbol is reserved for nondimensional variables. The air thermal conductivity, kinematic viscosity, thermal diffusion coefficient, and thermal expansion coefficient evaluated at $T*=(T0+T∞)/2$ were assumed to remain constant at $ka=0.026$ W m⁻¹ K⁻¹, $νa=1.6×10−5$ m² s⁻¹, $αa=2.24×10−5$ m² s⁻¹, and $βa=0.00347$ K⁻¹, respectively. Density variations with temperature were considered only in the buoyancy term of the vertical momentum equation according to the Boussinesq approximation [with $ρa=ρ(T*)=1.22$ kg m⁻³]. The Reynolds, Richardson, and Péclet numbers were $Re=wmd/νa=6000$⁠, $Ri=gβaΔTd/wm2=5.61×10−4$⁠, and $Pe=wmd/αa=4200$⁠, respectively. The gravity acceleration was $g δi2=−9.8$ m s⁻². The flow was integrated up to $t≈1.7$ s. The hydrodynamics is illustrated in Fig. 1(a), which shows a slice of the instantaneous w field at t = 0.75 s. This panel also shows the computational domain dimensions and some details of the cylindrical inlet and numerical mesh.

High-Order Spectral Element Methods (HO-SEMs) have been shown to be especially well suited for turbulent flows.⁴⁴ The Nek5000 HO-SEM solver⁴⁵ is used here to solve Eqs. (2)–(4) in the velocity–pressure form using the semi-implicit kth-order backward difference formula (BDFk)/kth-order extrapolation (EXTk) time-stepping in which the time derivative is approximated by a kth-order backward difference formula (BDFk), the nonlinear terms (and other forcing) are treated with a kth-order extrapolation (EXTk), and the viscous and pressure terms are treated implicitly (in this work, we used k = 3). This approach leads to a linear unsteady Stokes problem to be solved at each time step, which comprises a Helmholtz equation for each component of velocity (and temperature/scalar) and a Poisson equation for pressure (see, e.g., Sec. 2.2 of Mittal et al.⁴⁶). The numerical simulation presented here, with approximately 370 × 10⁹ mesh points, used 20 Intel Platinum 8168 nodes with 24 cores each interconnected with a 100 Gb/s Infiniband network. The average CPU time per time step is around 5 s. The simulation took $5.19×105$ CPU hours.

Laboratory measurements used to estimate the production of turbulent kinetic energy in self-similar momentum puffs by Glezer and Coles.⁴⁷ suggested values around $Π=100(I/ρf t5)1/2$ within the vortex ring, where I is the puff impulse (⁠$I=6×10−4$ N s in the present DNS study). Under the assumption of equal production and dissipation, the ratio between the Kolmogorov length scale η_K and the exit diameter can be estimated as $ηK/d≈0.01 t5/8$⁠. This result is compatible with estimations from DNSs of jets based on the measurements of Panchapakesan and Lumley⁴⁸ (see, for example, Boersma et al.⁴⁹). For the present Reynolds number based on the maximum velocity during the cough, $ηK/d≈6×10−4x/d$⁠. The present simulation shows that the initially laminar jet becomes completely turbulent, with fine scale activity, at t = 0.3 s (⁠$t̃=72$⁠). At this time, the puff is approximately located at z = 15 (⁠$z̃=15d$⁠) with values of the nondimensional Kolmogorov length scale of approximately $5×10−3$ and $10−2$ according to the criteria based on the measurements of Glezer and Coles⁴⁷ and Panchapakesan and Lumley,⁴⁸ respectively. The corresponding grid sizes at the jet axis at this position are $Δx̃=Δỹ≈0.009$ and $Δz̃≈0.04$⁠, which are of the same order of magnitude as the estimations. A detailed discussion on the cough hydrodynamics can be found in Fabregat et al.⁴³ 

Given the relatively initial small particle volume fraction (approximately $10−5$⁠) in real coughs and sneezes,⁵⁰ it is assumed that the particles do not affect the flow (one-way coupling). Under this hypothesis, the solver for the carrier phase hydrodynamics was coupled with a particle advection model to obtain the trajectories $sj(x→j,t)$ of $j=1,…,Np$ particles that represent the ejecta fluid. Although, immediately after being spewed, the ejected material typically undergoes filamentation and a breakup process due to capillary instabilities, the dispersed phase rapidly evolves into droplets and aerosols of spherical shape with sizes, typically ranging between 1 and several hundred μm in diameter.⁵⁰ Mostly composed of water,⁵¹ the dispersed phase density and thermal conductivity are set to $ρp=1000$ kg m⁻³ and $kp=0.606$ W m⁻¹ K⁻¹, respectively.

Particles are released using n_b = 200 batches equally spaced in time over approximately the entire duration of the cough (i.e., one batch released every 2 ms). To investigate the effect of droplet size and evaporation, particles of n_s = 7 different initial diameters (4, 8, 16, 32, 64, 128, and 256 μm), both evaporating and nonevaporating (n_e = 2), are released from n_l = 69 fixed seeding positions, resulting in a total number of particles $Np=nbnsnenl=193 200$ by the end of the expiratory event. The sizes of the particles were selected from the typical ranges reported for coughs.⁵⁰ According to the estimations of the Kolmogorov length scale indicated above (⁠$ηK≈2×10−4 m$⁠), only the largest particles of 256 μm, which are dominated by the gravitational settling, are larger than this length scale. The nondimensional location of the n_l particle seeding positions is shown in Fig. 1(b). The particles are seeded over a circular area centered at $(x̃, ỹ, z̃)=(0, 0, 0)$⁠. Depositions on the inlet vicinity boundaries are prevented by using an outer diameter for the seeding area equal to half the inlet diameter [red dashed line in Fig. 1(b)].

Due to the presence of the nonevaporable material in the ejected fluid mostly composed of proteins and salts, a lower bound for the particle diameter representing the size of the nuclei⁵² has been set to $0.3dp0$⁠, where $dp0$ represents the initial particle size. Thus, the fraction of remaining evaporable water ς in a particle can be defined as

where ς ranges between 1 at release time and 0 when no evaporable water is left. Due to the modest change in particle density as water evaporates, ρ_p has been assumed constant over the duration of the numerical experiment.

The position of an idealized spherical and smooth particle X_i can be written as

where the particle velocity U_i can be obtained by solving the balance between hydrodynamic drag, buoyancy, and thermophoresis forces,

The lift force, the pressure gradient force, and the virtual mass and Basset forces are neglected according to the relatively small density ratio between the fluid and the particle. The hydrodynamic drag accounts for the fluid resistance or friction a particle experiences as it moves with respect to the carrier phase and can be characterized by the particle Stokes number τ_p defined as the ratio of the particle and flow characteristic times. Small Stokes numbers are associated with particles that rapidly react to changes in accelerations in the underlying flow, while large Stokes are found in inertia-dominated particles. The Stokes number for a particle of nondimensional diameter $d̃p$ can be written as

where the last term is a correction of the Stokes drag when the particle Reynolds number $Rep=dp||ui−Ui||/νa$ does not remain small.⁵³ In this equation, C_c is a correction factor that depends on the Knudsen number Kn defined as

where $λ=6.8×10−8$ m is the free mean path for air at $T*$⁠.⁵⁴ 

The magnitude of the buoyancy force experienced by a particle floating in a fluid with different densities under a gravitational field with acceleration magnitude g can be written as

The thermophoretic force arises because of gradients in the temperature field across the particle characteristic length, d_p, for a spherical particle.⁵³ The prefactor of the nondimensional temperature gradient, n_th, is defined as

where the factor K_tp, which is defined as⁵⁵ 

contains three constants set to $Cs=1.17, Cm=1.14$⁠, and $Ct=2.18$⁠.

Note that the negative sign in n_th indicates that the force points in the direction of maximal temperature decay.

After solving for the carrier phase velocity $ũi$ by integrating in time the hydrodynamics [Eqs. (2)–(4)], the dispersed phase [Eqs. (6) and (7)] is temporally integrated using the Exponential-Lagrangian Tracking Scheme (ELTS) derived by Barton.⁵⁶ While classical explicit methods exhibit prohibitive computational costs due to numerical instability constraints for relatively small values of τ_p (e.g., $Δt<2τp$ for fourth-order Runge–Kutta methods), the ELTS is inherently stable and allows us to use the same time step for both the carrier and disperse phases.

The nondimensional equation for the rate of change in particle size due to evaporation of the aqueous fraction, which can be obtained from the unsteady mass balance of the particle,⁵⁷ can be written as

where ρ_s is the saturation water vapor concentration at the droplet surface, ρ_f is the local ambient water vapor concentration, D_v is the water vapor diffusivity in air, and $Scv=νaDv$ is the Schmidt number for water vapor.

The evolution of the particle temperature, derived from the thermal energy balance of the particle,⁵⁷ can be written as

where $νp=1.015×10−6$ m² s⁻¹ is the water kinematic viscosity at $T*$⁠, θ_p and θ_f are the particle and ambient fluid nondimensional temperatures, $Cpp=4179.6$ J kg⁻¹ K⁻¹ is the liquid water heat capacity, $kv=0.026$ W m⁻¹ K⁻¹ and $kp=0.606$ W m⁻¹ K⁻¹ are the water vapor and liquid water thermal conductivities, $αp=kp/(ρpCpp)=1.45×10−7$ m² s⁻¹ is the liquid water thermal diffusivity, $ΔHv=2.257×106$ J kg⁻¹ is the water enthalpy of vaporization, and $Prp=νp/αp=7.0$ is the dispersed phase Prandtl number. Equation (14) assumes that the Sherwood and Nusselt numbers are constant (⁠$Sh=Nu=2$⁠)⁵³ given the small particle Reynolds numbers involved. We found that the use of a higher Reynolds correlation does not influence significantly the results. Specifically, for the largest particles considered, of 256 μm, falling freely from a height of 1.5 m, the differences of the diameter using $Sh=Nu=2$ or the correlation proposed by Ranz and Marshall⁵⁸ when they reach the floor are less than 1%.

The saturation water vapor concentration ρ_s is determined, using the ideal gas law, as

where the water vapor pressure p_s in Pa is computed from the correlation⁵⁹ 

where T is in K. Recall that $T=T∞+ΔTθ̃$⁠.

The local ambient water vapor concentration ρ_f is determined as

where the relative humidity ξ can, then, be calculated using⁵⁹ 

with the total pressure taken to be $p0=101 325$ Pa. Since both heat and water vapor transport equations are only distinguished by their Peclet number and given the similitude between the Prandtl and Schmidt numbers (Pr = 0.7 and Sc = 0.6, respectively), differences in the transport between both scalars due to diffusive effects have been neglected and the temperature field has been used as a proxy for the water content in the air. Thus, the mass fraction of water in the air x_v can be computed assuming a linear relationship between temperature and the mass fraction of water vapor (i.e., for $T=T∞=15$ °C or $θ̃=0, xv=0.007$ and for $T=T0=34$ °C or $θ̃=1, xv=0.0273$⁠),

Given the relatively small particle volume fraction, this model neglects the heat transfer from the dispersed to the carrier phase.

At each time step, the disperse phase transport equations [(6) and (7)] are integrated in time with the updated values of particle diameter and particle temperature obtained by integrating in time the ordinary differential equations [(14) and (15)] with the same BDF3 scheme used to integrate the hydrodynamics.

The simulated expiratory event allows us to investigate the role of initial particle size and evaporation on the resulting aerosol cloud dispersion under typical temperature and relative humidity conditions. Inspection of the relevance of each force in Eq. (7) demonstrates that the contribution due to thermophoresis is negligible in comparison to hydrodynamic drag and buoyancy terms.

Temporal evolution of instantaneous averages across cloud particles, represented by $⟨·⟩$⁠, for the fraction of remaining evaporable water ς defined in Eq. (5) and the streamwise velocity W are shown in the top and bottom panels of Fig. 2, respectively. The time to complete evaporation is found to strongly depend on the initial particle size. By the end of the experiment, particles of initial diameter 16 μm and smaller have completely evaporated, while droplets of 64 μm and larger still retain more than 75% of the free water at release.

As shown in the bottom panel of Fig. 2, the horizontal speed of the cloud peaks at approximately t_m (red dashed vertical line) when exhalation reaches its maximum velocity. By the end of the numerical experiment, the cloud has lost most of its momentum with the values of the averaged vertical velocity component (⁠$⟨W⟩$⁠) around 0.05 m s⁻¹. This result suggests that while particle transport is initially controlled by the flow induced by the expiratory event, once the momentum dissipates and the cloud velocity declines to typical indoor conditions⁶⁰ with an average velocity of about 0.1 ms⁻¹, the dispersion processes, thereafter, will most likely be dominated by the ambient air currents.

The temporal evolution of the cloud centroid $⟨Xi⟩$ for each particle size and evaporation mode is shown in Fig. 3. The top and bottom panels show the results for the drag-dominated (4–32 μm) and gravity-dominated (64–256 μm) particles, respectively. Each symbol represents the cough end time t_c for each nonevaporative case (for the sake of clarity, the very similar cough end time for the evaporative case is not shown).

The small differences between the cloud centroid trajectory for evaporating (dashed lines) and nonevaporating (solid lines) particles under 32 μm in diameter suggest that evaporation has a negligible impact on the aerosol dispersion. The drag-dominated transport of these aerosols exerted by the background flow can only intensify as the particle diameter further decreases due to evaporation. Notably, over the duration of the numerical experiment, both the 128 and 256 μm particles have mostly reached the boundaries of the domain located at $r=R̃=25$ (see Fig. 1). The largest ones, with larger inertia, are the particles that travel the farthest before reaching to bottom domain limits.

Conversely of the particle distribution, particles with a diameter above 32 μm also exhibit negligible differences in their cloud centroid trajectory. In this case, and as shown in Fig. 2, the shrinking in size due to evaporation over the duration of the experiment is too small to significantly change the particle dynamics dominated by the gravitational action.

Notably, evaporation is found to significantly modify the path followed by the 32 μm cloud. While nonevaporative particles follow a quasi-parabolic trajectory similar to that described by larger particles, the evaporative cloud counterpart reverses its downward trajectory. This behavior suggests that for this specific size, evaporation leads to a transition from gravity to drag dominated dispersion.

To compare both phases, the trajectory followed by the frontal region of the thermal puff produced by the cough⁴³ is included in Fig. 3. As expected, as the particle size decreases, the transport contribution due to hydrodynamic drag increases, and both the puff front and particle cloud paths tend to overlap. As gravitational effects increase for particles above 16 μm, the trajectory departure from the rising thermal puff strengthens.

Figure 4 shows the cloud size, estimated from the standard deviation of the particle positions with respect to the centroid σ_i,

for particles of 4–32 μm (left panels) and 64–256 μm (right panels) in diameter for both evaporating (dashed) and nonevaporating (solid) modes. The landing of most of the 128 and 256 μm particles on the computational domain boundary, where they remain deposited eventually, leads to constant values of σ_i in all directions. The start and end times of this deposition process roughly span from the peak to the nearly zero values of σ_y for these two particle sizes.

Estimated from σ_x, the spanwise width inversely grows with the particle size for the 4–32 μm range with very similar values in σ_y and σ_z.

To gain insight into the effects of evaporation and estimate the impact of the droplet size decrease on the particle dispersion, Fig. 5 shows the probability density function of the particle position difference between the evaporative and nonevaporative cases at the latest simulation time. The results for the distribution of vertical positions (top row of Fig. 5) suggest that evaporation leads to significant rises in y only for release diameters above 16 μm with values above 64 μm (not shown here for simplicity) being negligible. The results indicate that there is a significant skewness value in the vertical position difference distribution for 32 and 64 μm particles with a mean vertical rise of approximately 4 cm for the evaporative case.

The distribution of the z position difference (bottom row of Fig. 5) suggests that the impact on the horizontal range for particles smaller than 32 μm and larger than 64 μm (not shown) in diameter is negligible. The distribution for 32 and 64 μm, although slightly skewed toward positive values, still exhibits zero average.

Particle position distribution in each spatial direction for the evaporative case at the end of the numerical experiment is shown in Fig. 6 for all particle sizes (left panels) and only the five smaller diameters (right panels). As indicated by the standard deviation results in Fig. 4, lateral cloud spreading increases as the particle diameter decreases. Differences in the vertical position distributions are negligible between particles from 4 to 16 μm. As the particle diameter increases, the cloud falls due to gravity intensify with 64 μm droplets still afloat and with 128 and 256 μm droplets completely landed on the domain boundaries. The largest particles, capable of retaining more injection momentum, exhibit the largest horizontal reach up to approximately 125 cm from the source before landing. Interestingly, the 4–16 μm particles, capable of staying afloat within the thermal puff, travel farther than the rest of the diameters. The intermediate size of 64 μm, large enough for the fall outside the thermal puff relatively early but too small to retain the injection momentum, exhibits the smallest horizontal spread.

The 4–32 μm and 64–256 μm particle instantaneous positions for the nonevaporative case are shown in Figs. 7(a) and 7(b), respectively. Markers are scaled according to particle size and colored by the batch number with the first batch released at t = 0 and the last one n_b = 200 at approximately t = 0.4 s. For reference, the inlet section has been projected as a cylinder spanning over the entire domain. The frontal puff region in Fig. 7(a) is characterized by the presence of particles released at the peak injection velocity (batch 70 in black). While particles released at early times during the initial stages of the cough (for instance, in white) tend to fall due to gravitational effects, those injected with similar velocity at the end of the expiratory event (in red) remain afloat dragged by the growing puff.

In the case of large particles shown in Fig. 7(b), both whitish and reddish particles capable of rapidly falling outside the puff due to gravity exhibit similar vertical position distributions. Both 128 and 256 μm appear at the bottom of the computational domain after describing quasi-parabolic trajectories.

Figures 8, 9, and 10 (Multimedia view) show the of nondimensional particle position at the latest time $t≈1.68$ s for both evaporative (right panel) and nonevaporative (left panel) particles of initial diameters 8, 32, and 64 μm, respectively. The marker color corresponds to the value of ς defined in Eq. (5).

We computed the potential viral concentration in the air using the positions of the particles at the maximum time reached by the simulation (⁠$t=1.68 s$⁠). As indicated above, at this time, the momentum of the puff has decreased to velocity values typical of indoor conditions. We have counted the particles in bins of $2×2×2$ cm³, and the potential viral concentration has been calculated as the initial volume of the expelled droplet (nl) in a given bin multiplied by the part per thousand of the droplets of a given diameter in the bin, according to a typical droplet size distribution of expiratory events reported by Duguid,⁵⁰ and divided by the volume of the bin (8 cm³). The potential viral concentration nl/1000 cm³ in the air resulting from the injection of 1000 ejected particles is shown in Fig. 11 for the evaporative (left panels) and nonevaporative (right panels) cases for all particles (top row), 64–256 μm (central row) and 4–32 μm (bottom row). The total viral load per air volume can be computed multiplying the potential viral concentration by the number of droplets expelled divided by 1000 and by the concentration of virus per volume of expelled material. The oral fluid average virus RNA load is 7 copies per nl with maximum up to $2.35×103$ copies per nl.⁵² For COVID-19, the minimum infectious dose is unknown, but it is estimated to be around a hundred virus particles.⁶¹ Considering that one cough releases around 5000 droplets,⁵⁰ a value of the potential viral concentration of 0.01 corresponds to 350 copies per liter of air. The volume of inhaled air during normal breathing is about 0.5 l, and consequently, the inhalation of the air in the regions marked in red in Fig. 11 will produce at least the intake of 175 virus copies per inhalation.

As suggested by the results in Fig. 5, the overall effect due to evaporation is limited. Despite the larger horizontal range of the particle cloud, the low particle concentration in the puff front results in a reduced potential viral concentration in that region. Nonzero potential viral concentration levels exclusively due to large particles are mostly found at the location lower than the source (dashed horizontal black lines). On the other hand, small particles are responsible for vertically wider risk maps due to their capacity to rise above the source height level.

Overall, for the conditions considered in this simulation, the horizontal range of the viral concentration cloud before typical ambient air currents will likely dominate the aerosol dispersion is around 80 cm from the emitter.

Characterized by its transition from a laminar and then a turbulent jet into a buoyant puff, the numerical simulations of an idealized mild cough stressed the complex dynamics of the flow produced by expiratory events. Equipped with the fully resolved hydrodynamics and using a transport model for inertial particles, the results presented here demonstrate the impact of the background flow on the dispersion and fate of airborne droplets and aerosols produced by the breakup of ejected mucosa in rapid air exhalations.

Due to the significant impact of the particle size on the ability of the background flow on transporting it, the disperse phase model accounts for diameter reduction due to evaporative effects. Thus, droplets and aerosols produced by a cough are simulated by continuously seeding idealized spherical particles of different initial diameters and evaporation rates over the duration of the expiratory event. Individual trajectories are combined to obtain cloud statistics that allow us to elucidate the overall dispersion properties of the dispersed phase.

On the one hand, evaporation is found to be relevant only for a certain range of particle diameters for which the size reduction due to water content loss produces a transition of the dominant transport contribution from gravity to drag-dominated regimes over the span of the particle excursion.

We show that the cloud dynamics of aerosols under 64 μm is characterized by an upward vertical drift produced by the buoyant puff. In contrast, larger particles tend to settle due to gravitational effects. Combining this detailed particle dispersion information with several reasonable assumptions on the viral load of typical ejecta fluid, we derived a viral concentration map for idealized environments that can be used to assess the infection risk under typical ambient conditions. In the future, the techniques used in the current study will be extended to other expiration events, e.g., talking, which has been postulated to be the primary mechanism for asymptomatic transmission. In the future, the importance of evaporation on long-term Heating, Ventilation, and Air Conditioning (HVAC) induced dispersion of aerosols is also a topic that will be explored.

This work was funded by Spanish Ministerio de Ciencia, Innovación y Universidades through Grant Nos. DPI2016–75791-C2–1-P and RTI2018–100907-A-I00 (MCIU/AEI/FEDER, UE) and also by the Generalitat de Catalunya through Grant No. 2017-SGR-1234. The authors thank the University of Catalonia Services Consortium (CSUC) and the Red Española de Supercomputación (RES) for providing the computational resources used in this work. Computational resources provided by the Covid-19 HPC Consortium through time on Frontera (TACC) and Blue Waters (NCSA) were used to test and smooth the computational mesh. S.D.'s participation was supported by the DOE Office of Science through the National Virtual Biotechnology Laboratory (NVBL), a consortium of DOE national laboratories focused on response to COVID-19, with funding provided by the Coronavirus CARES Act.

The authors report no conflict of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.