A main route for SARS-CoV-2 (severe acute respiratory syndrome coronavirus) transmission involves airborne droplets and aerosols generated when a person talks, coughs, or sneezes. The residence time and spatial extent of these virus-laden aerosols are mainly controlled by their size and the ability of the background flow to disperse them. Therefore, a better understanding of the role played by the flow driven by respiratory events is key in estimating the ability of pathogen-laden particles to spread the infection. Here, we numerically investigate the hydrodynamics produced by a violent expiratory event resembling a mild cough. Coughs can be split into an initial jet stage during which air is expelled through mouth and a dissipative phase over which turbulence intensity decays as the puff penetrates the environment. Time-varying exhaled velocity and buoyancy due to temperature differences between the cough and the ambient air affect the overall flow dynamics. The direct numerical simulation (DNS) of an idealized isolated cough is used to characterize the jet/puff dynamics using the trajectory of the leading turbulent vortex ring and extract its topology by fitting an ellipsoid to the exhaled fluid contour. The three-dimensional structure of the simulated cough shows that the assumption of a spheroidal puff front fails to capture the observed ellipsoidal shape. Numerical results suggest that, although analytical models provide reasonable estimates of the distance traveled by the puff, trajectory predictions exhibit larger deviations from the DNS. The fully resolved hydrodynamics presented here can be used to inform new analytical models, leading to improved prediction of cough-induced pathogen-laden aerosol dispersion.

Outbreaks of infectious diseases have profoundly affected human societies on several occasions through history. The ongoing COVID-19 (coronavirus disease) pandemic is the most recent example of the extraordinary impact on health, economics, and sociopolitics of highly contagious pathogens in an increasingly globalized world. By mid-December 2020, SARS-CoV-2 has infected more than 75 × 106 people worldwide killing more than 1.7 × 106 of them.1,2

Airborne droplets and aerosols are released into the air when infected people cough, sneeze, or talk.3,4 These particles constitute the primary transmission route for SARS-CoV-2 (severe acute respiratory syndrome coronavirus), Measles morbillivirus, Mycobacterium tuberculosis, chickenpox, and influenza virus, among others. Understanding the airflow generated by respiratory events is key in predicting how pathogen-laden aerosols are dispersed in the environment and providing valuable information to design measures and interventions aimed at hindering airborne transmission. In the classical approach, dispersion of large aerosols or droplets (>50 μm in diameter), mostly controlled by gravity, is characterized by ballistic trajectories, short residence times, and limited horizontal range. On the other end of the diameter range, small droplets or aerosols (<50 μm) are mostly dispersed by the action of hydrodynamic drag, leading to the formation of aerosol clouds capable of remaining afloat over long periods of time and reaching relatively larger distances. Crucially, enhanced transport of these tiny aerosols by relatively weak background currents as those generated by heating, ventilation, and air conditioning (HVAC) has been linked to transmission between individuals separated over distances much larger than the 2 m-rule suggested by social distancing measures.5 While the effect of HVAC systems on far-field transmission of respiratory viruses is currently being debated and studied,6,7 the general agreement is that during close contacts the short-range airborne route is the dominant mode of respiratory virus transmission.8,9

Recent studies have used experiments10,11 and numerical simulations12,13 to analyze various aspects of the jet flow and the accompanying aerosol transport for different expiratory events including sneezing,10 coughing,12 talking, and breathing.11 Numerical simulations with different turbulence modeling techniques have been used to study aerosol and droplet dispersion during violent expiration events. Two of the first studies that explored the effect of factors like ambient wind12 and masks14 were conducted by Dabouk and Drikakis using the Reynolds averaged Navier–Stokes (RANS) modeling approach. Dabouk and Drikakis used the compressible multiphase mixture RANS equations and the κ ω turbulence model to simulate the carrier expiratory jet of the cough12 and modeled the saliva droplets as Lagrangian particles. One of the main findings of Dabouk and Drikakis was the determination of the wind-speed influence on the airborne droplet traveled distance with droplets traveling up to 6 m for background wind speeds of 15 km h−1. Busco et al. also simulated the jet produced during a sneeze using the compressible RANS equations with the continuous phase being modeled as a mixture of air and water-vapor.10 The turbulence closure used in their model was κ ϵ , and the saliva droplets were modeled as Lagrangian particles. Even though sneezing is not among the primary symptoms of Covid-19, Busco et al. elucidated the importance of capturing the angle of release of the expiration jet in order to accurately model the dispersion of airborne droplets.

Fontes et al. also used numerical simulations to investigate the fluid dynamics of a sneeze.13 In particular, they focused on the impact of human physiological factors (e.g., illness, stress condition, anatomy, etc.) on droplet dispersion. Fontes et al. modeled the expiratory jet as a gas flow under the Eulerian framework and the dispersed droplets as Lagrangian particles. Turbulence in the flow was modeled using detached eddy simulations (DES), which is based on the combination of unsteady Reynolds averaged Navier–Stokes (URANS) for the flow within the boundary layer and large eddy simulations (LESs) in the outer region. Fontes et al. found that the change in the nasal and buccal passages have a significant impact on the distance traveled by the droplets with nasal-obstruction resulting in an increase of up to 60%. Additionally, physical characteristics of saliva were found to change different aspects of the spray generated after the sneeze. Pendar and Pascoa conducted LES to model the dispersion of saliva droplets during violent expiratory events15 like coughing and sneezing. They also used an Eulerian approach for the carrier air jet and the Lagrangian approach for the droplets. The jet was modeled using the compressible Navier–Stokes equations and an LES turbulence closure where the subgrid-scale stresses were calculated using the one equation eddy-viscosity model. Pendar and Pascoa found that for a strong sneeze, large droplets (540 μm in diameter) could travel up to 4 m and particles could take up to 3 s to settle down in the absence of a background flow. They also observed that bending of the head during sneezing could reduce the droplet horizontal range by 22% and face masks could further reduce the aerosol traveled distance to 0.6 m. Additionally, they concluded that droplet transport was enhanced by turbulence. This particular finding shows the importance of accurately modeling and understanding the turbulent characteristics of expiratory jets.

Wang et al. used alternative turbulent models based on random walk to approximate the effect of turbulence fluctuations on aerosol dispersion and study the effect of environmental factors on droplet transport and dispersion.16 Renzi and Clarke showed that the dynamics of the expiratory jet could be modeled by extending the theory of buoyant vortex rings.17 They coupled their integral model of the continuous phase with a Lagrangian particle-tracking model for the droplets and used it to explore the effect of the initial condition on the distance traveled by the droplet cloud. They also observed that the vortex in their model plays a key role in keeping droplets afloat, providing additional motivation to study in detail the turbulent structure of the expiratory jet.

A survey of the literature shows that most of the existing numerical studies focus primarily on the dispersal of the droplets due to expiration events and do not provide enough insight into the turbulent structure of the puff. Thus, we conduct a direct numerical simulation (DNS) of the flow produced by an idealized, relatively violent respiratory event resembling a mild cough, with the motivation of understanding its turbulent dynamics. The finite injection of exhaled air generates a jet that penetrates into the initially isothermal and stagnant environment. Due to the temperature difference between the injected and the ambient air, once the cough ceases, the flow transitions into a thermal puff that bends in the vertical as turbulence decays due to dissipation.

Scorer18 experimentally investigated the evolution of isolated masses of buoyant fluid or thermals by releasing finite volumes of heavy solutions into a water tank. The results were used to validate an analytical model in which the isolated mass of buoyant fluid kept its linear momentum by decelerating as its volume grew through entrainment of lighter ambient fluid. The mean value of the entrainment coefficient, relating the thermal radius r and the vertically traveled distance form the source y, i.e., r = α y , was found to be α = 0.25 (or n = 4 = 1 / α in the original paper). The mean prefactor η in the relation between radius and volume of the thermal, V = η r 3 , was estimated18 to be η = 3.

Puffs, occurring when the isolated mass of fluid has some initial momentum, were investigated by Richards19 who experimentally produced both axial (point source) and cylindrical (line source) thermals reporting estimations for α between 0.13 and 0.53. Experiments on axial thermal puffs by Richards20 were used to derive a relation between the traveled distance z and time t of the form z 4 = C n 3 M g t 2 / ρ , where C is a constant, M is the mass excess, g is the gravity acceleration, and ρ is the ambient density.

By directing the flow with some angle β0 with respect to gravity aligned with the y direction, Richards21 derived a model for inclined puffs of constant total buoyancy by solving the momentum equations P z = | P | 0 cos ( β 0 ) and P y = | P | 0 sin ( β 0 ) + | M g | t , where P is the puff impulse. Aimed to better understand puffs and thermals in the atmosphere, Scorer22 reported a parameter-free model for predicting the front location, velocity, and radius of a turbulent puff.

Using the analytical solution for inclined thermal puffs derived by Richards,21 Bourouiba, Dehandschoewercker, and Bush23 investigated the flow produced by violent expiratory events by comparing the model predictions with experimental results obtained by horizontally injecting particle laden fresh water payloads into a tank filled with salty water. The puff trajectory and shape were obtained by tracking the frontal region and fitting an spheroid (two identical semi-axes) for successive snapshots. In contrast to ellipsoidal shapes characterized by three independent semi-axis lengths, the semi-axis most aligned with the vertical, used to estimate the puff radius r, was assumed to be equal in length to the one in the spanwise direction.

In this work, a similar approach is used to elucidate the three-dimensional trajectory and shape of a puff produced by a numerical cough. Fully resolved results on the flow hydrodynamics suggest that, along with the self-similar hypothesis in the Scorer22 and Bourouiba, Dehandschoewercker, and Bush23 models and the transient behavior of a cough air injection, the spheroidal assumption may contribute to explain the observed differences between the experimental and analytical attempts to characterize the flow produced by a violent expiratory event and the DNS predictions. Given the major role played by the background fluid hydrodynamics in the aerosol dispersion, the new insight presented here can be used to improve our current predictive capabilities of infection risk by airborne-transmitted pathogens.

Although Direct Numerical Simulations have been carried out in the past,24 to our knowledge, this is the first time that details on the temporal evolution of puff centroid, entrainment coefficient, and front topology have been reported.

A mild cough is modeled as a finite injection of air into an initially isothermal and quiescent environment. Exhaled and ambient air temperatures are T 0 = 34 °C and T = 15 °C, respectively. The transient nature of the cough flow is captured25 by imposing an inlet transient velocity w 0 ( t ) that ramps up linearly from zero to the peak velocity w m = 4.8 m s−1 at t m = 0.15 s and then ramps down linearly to zero at t c = 0.4 s, i.e.,

w 0 ( t ) = { w m t m t , 0 t < t m w m w m t c t m ( t t m ) , t m t t c 0 , t > t c .
(1)

Thus, air flow generates an accelerating ( 0 t < t m ) and a decelerating ( t m t t c ) jet that eventually evolves into a thermal puff when injection ceases for t > t c .

Air thermal conductivity, kinematic viscosity, thermal diffusion coefficient, and thermal expansion coefficient evaluated at T a = ( T 0 + T ) / 2 are assumed to remain constant at k f = 0.026 m 1 K 1 , ν = 1.6 × 10 5 m 2 s 1 , α = 2.24 × 10 5 m 2 s 1 , and β = 0.003 47 K 1 , respectively. Density variations with temperature are considered only in the buoyancy term of the vertical momentum equation according to the Boussinesq approximation [with ρ a ( T a ) = 1.22 kg m 3 ]. The air exit is modelled25 as a cylindrical pipe of diameter d  =  0.02 m and length H p = 0.04 m. The subject is assumed to remain steady with a pitch angle perpendicular to gravity at all times. Under the above stated conditions, the airflow during the cough can be modeled as an incompressible flow.26 

The conservation equations for mass, momentum, and energy for incompressible flows can be written as

u i x i = 0 ,
(2)
u i t + u j u i x j = 1 ρ a p x i + ν 2 u i x j x j + g β ( T T ) δ i 2 ,
(3)
T t + u j T x j = α 2 T x j x j ,
(4)

where t is the time, p is the pressure, T is the temperature, δij is the Kronecker delta, and u i = ( u , v , w ) is the velocity field with coordinates x i = ( x , y , z ) in the spanwise, gravity-aligned, and streamwise directions. Using the exit diameter d, the inlet peak velocity w m , and the maximum temperature difference Δ T = T 0 T = 19 ° C as space, velocity, and temperature scales, the Reynolds, Richardson, and Péclet numbers are R e = w m d / ν = 6000 , R i = g β Δ T d / w m 2 = 5.61 × 10 4 , and P e = w m d / α = 4200 , respectively. The gravity acceleration is g δ i 2 = 9.8 m s 2 , and the temperature perturbation is defined as θ ̃ = ( T T ) / Δ T . Note that the tilde symbol is reserved for non-dimensional variables. The variation of the physical properties with the water vapor concentration is neglected. This approximation is reasonable under the current conditions considered. For example, for exhaled air at T = 34 °C and relative humidity R H = 85 % ,27  ρ = 1.130 kg m 3 , and for ambient conditions at T = 15 °C and R H = 65 % , ρ = 1.220 kg m 3 . However, for air at T = 34 °C and R H = 65 % , ρ = 1.134 kg m 3 .

The computational domain dimensions, the coordinate system, and the mesh details are illustrated in Fig. 1 using non-dimensional variables. The inlet boundary conditions in the cylindrical injection section of diameter ( d ̃ = 1 ) and length ( H ̃ p = 2 ) are ( u ̃ , v ̃ , w ̃ ) | x ̃ , y ̃ , z ̃ = 2 , t ̃ = ( 0 , 0 , w ̃ 0 ( t ̃ ) ) and θ ̃ | x ̃ , y ̃ , z ̃ = 2 , t ̃ = 1 , with boundary conditions for the outer walls set to no-slip and adiabatic. An annular indent of the Gaussian profile with the center at z ̃ d = 1 / 2 , a depth of h ̃ d = 0.05 , and a width of σ ̃ d 2 = 0.01 has been used to mimic the complicated passage of exhaled air across the human mouth, cause boundary layer separation, and facilitate the transition to turbulence. The main domain consists of a cylinder of diameter D ̃ = D / d = 50 and length H ̃ = H / d = 80 . Outflow boundary conditions have been imposed at z = H and at the outer cylindrical wall of radius R = D / 2 , which have also been considered adiabatic. These dimensions ensure that the cough flow does not interact with the boundaries over the temporal extent of the simulation. Initial conditions over the entire computational domain are u ̃ = v ̃ = w ̃ = θ ̃ = 0 .

FIG. 1.

Section of the computational domain showing an instantaneous w velocity component field at t = 0.75 s. The inset shows the spectral element mesh (in red) and the LGL nodes (in black).

FIG. 1.

Section of the computational domain showing an instantaneous w velocity component field at t = 0.75 s. The inset shows the spectral element mesh (in red) and the LGL nodes (in black).

Close modal

The results presented here have been obtained using Nek5000,28 an open-source high-order spectral element method (HO-SEM)-based solver for the incompressible Navier–Stokes equations (INSEs). The SEM is a high-order weighted residual method that combines the geometric flexibility of finite elements (Ω is decomposed into K smaller elements) with the rapid convergence of spectral methods. The basis functions in SEM are Nth-order tensor-product Lagrange polynomials on the Gauss–Lobatto–Legendre (GLL) quadrature points inside each element, which lead to fast operator evaluation and low operator storage cost.29 In the SEM-based solver Nek5000, the unsteady INSEs are solved in velocity-pressure form using semi-implicit BDFk/EXTk timestepping 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. 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 in the study by Mittal30) Based on this approach, SEM has proven to be well suited for turbulent flows.31–33 

For the current DNS, we model the domain using K = 215 820 elements with N = 11th order polynomials for the solution. The total number of mesh nodes is K × ( N + 1 ) 3 370 × 106. The simulation has been run on 20 nodes of a Central Processing Unit (CPU) cluster, interconnected with a 100 Gb/s Infiniband network. Each node contains two Intel Platinum 8168 CPUs with 24 cores each. The average CPU time per time step is around 5 s. The simulation took 5.19 × 10 5 CPU h to reach t = 1.68 s.

Glezer and Coles34 estimated the production of turbulent kinetic energy of a self-similar momentum puff using laboratory measurements. They reported typical values of production of turbulent kinetic energy around Π = 100 ( I / ρ f t 5 ) 1 / 2 within the vortex ring, where I is the puff impulse ( I = 6 × 10 4 N s in the present DNS study). Assuming that production equals dissipation, the ratio between the Kolmogorov length scale ηK and the exit diameter is η K / d 0.01 t 5 / 8 . This result is compatible with the estimation usually reported for DNSs of jets based on the measurements of Panchapakesan and Lumley35 (see, for example, the study by Boersma, Brethouwer, and Nieuwstadt36) According to these authors, for the present Reynolds number based on the maximum velocity during the cough, η K / d 6 × 10 4 x / d . The 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 at x = 15 ( x ̃ = 15 d ) and the estimations of the non-dimensional Kolmogorov length scale are 5 × 10 3 and 10 2 according to the criteria based on the measurements of Glezer and Coles34 and Panchapakesan and Lumley,35 respectively. The grid sizes at the jet axis at this position are Δ x ̃ = Δ y ̃ 0.009 and Δ z ̃ 0.04 , which are of the same order of magnitude as the estimations.

The flow hydrodynamics is illustrated in Figs. 2 and 3 that show a detail of the x = 0 plane for the instantaneous velocity magnitude (in m s−1) and temperature (in °C) at six different times, including the instants of peak velocity ( t = t m = 0.15 s) and the end of the cough ( t = t c = 0.40 s), respectively.

FIG. 2.

Detail of the slice at x = 0 of the velocity magnitude field (in m s−1) at t = 0.15 , 0.25 , 0.30 , 0.40 , 0.75 , 1.5 s. Note that t = t m = 0.15 s and t = t c = 0.40 s correspond to the peak and cough ending times, respectively.

FIG. 2.

Detail of the slice at x = 0 of the velocity magnitude field (in m s−1) at t = 0.15 , 0.25 , 0.30 , 0.40 , 0.75 , 1.5 s. Note that t = t m = 0.15 s and t = t c = 0.40 s correspond to the peak and cough ending times, respectively.

Close modal
FIG. 3.

Detail of the slice at x = 0 of the temperature field (in °C) at t = 0.15 , 0.25 , 0.30 , 0.40 , 0.75 , 1.5 s. Note that t = t m = 0.15 s and t = t c = 0.40 s correspond to the peak and cough ending times, respectively.

FIG. 3.

Detail of the slice at x = 0 of the temperature field (in °C) at t = 0.15 , 0.25 , 0.30 , 0.40 , 0.75 , 1.5 s. Note that t = t m = 0.15 s and t = t c = 0.40 s correspond to the peak and cough ending times, respectively.

Close modal

The distance traveled by a horizontal buoyant puff released in quiescent, neutral ambience can be written as21 

s ( t ) = { 4 η α 3 ρ f P 0 2 2 B 0 [ B 0 t P 0 1 + ( B 0 t P 0 ) 2 + ln ( B 0 t P 0 + 1 + ( B 0 t P 0 ) 2 ) ] } 1 / 4 ,
(5)

where P 0 = ρ f V 0 w ¯ and B 0 = Δ ρ V 0 g are the initial linear momentum and the initial buoyancy of the puff with initial ejected volume V0 and average velocity w ¯ , respectively, α is the entrainment coefficient, and η is a constant that depends on the shape of the puff. Equation (5) assumes that the velocity distributions inside the puff are similar and that the momentum of the puff remains constant during its dispersion. The linear momentum is conserved by the entrainment of the quiescent ambient fluid, resulting in a puff that increases in mass and decelerates as it penetrates into the ambient. In inclined puffs and thermals, where the puff may bend due to buoyancy effects, the radial spread of the puff can be written23 as r ( t ) = α s ( t ) .

Previously reported18,19 values of α for thermals are 0.25 and between 0.13 and 0.53. In their experiments, Bourouiba, Dehandschoewercker, and Bush23 found smaller values ranging between 0.09 α 0.18 for the jet stage and 0.015 α 0.037 for the horizontal buoyant puff phase.

The behavior of a horizontal buoyant puff released with an initial momentum during a short initial period of time ( 0 t t j , end ) can be understood as an initial turbulent jet that evolves to a puff.23,37 During the jet stage ( 0 t t j , end ) the axial location zj, the axial velocity Wj and the radial spread Rj are38 

z j ( t ) = ( 12 K ) 1 / 2 ( W j 0 R j 0 ) 1 / 2 t 1 / 2 ,
(6)
W j ( t ) = 6 W j 0 R j 0 K z j ( t ) ,
(7)
R j ( t ) = R j 0 + ( z j ( t ) z j 0 ) n j .
(8)

For the puff stage ( t t j , end ), the axial location zp can be expressed as22 

z p ( t ) = z j , end + ( 4 W j , end R j , end 3 n p 3 ) 1 / 4 ( t 1 / 4 t j , end 1 / 4 ) .
(9)

The subscripts 0, j, and p indicate initial state, jet, and puff, respectively. In Eqs. (6) and (7), K = 0.457 is a constant for jets,38 R j 0 is the radius of the orifice of the jet, and W j 0 is the exit velocity of the jet. z j 0 in Eq. (8) is the position of the virtual origin of the jet, which can be neglected under the conditions of the simulation. The parameter n in the jet and puff equations defines the angle of spread, θ, n = 1 / tan θ = x / r . For jets,22 θ j 11.3 ° and nj = 5 and for puffs,22 θ p 14.0 ° and np = 4.

The trajectory followed by the thermal puff is obtained following an analogous approach to that used by Bourouiba, Dehandschoewercker, and Bush23 who determined the ellipse that properly enclosed the front region of the dyed puff in consecutive photographic snapshots. The longest of the semi-axes, mostly perpendicular to z along the entire experiment, was used as a measure of the puff radius r.

For the present numerical results, each numerical snapshot ( k = 1 224 ) of the instantaneous θ ̃ field is first interpolated into a Cartesian grid with similar average resolution to that used in the computational mesh, then integrated along x, and finally binarized using an indicator function ik with a prescribed tolerance of ε = 10 3 , i.e.,

i k = { 1 , θ ̃ k x ε 0 , θ ̃ k x < ε .
(10)

The instantaneous temperature integrated along x, θ ̃ x , is thought to mimic the photographic images taken during the experiments by Bourouiba, Dehandschoewercker, and Bush.23 

Using the contour of ik, the puff front is fitted by a 3D ellipsoid with the centroid ( c x , c y , c z ) and semi-axes ( σ 1 , σ 2 , σ 3 ) .39 To mimic the analysis in Ref. 23, we define r as the longest projected semi-axis on the plane yz.

The puff front temperature and vertical velocity are shown in the top and bottom panels of Fig. 4, respectively. Results suggest that, after the end of the cough event (marked by the vertical dashed blue line), temperature within the puff front decays exponentially as turbulent mixing entrains fresh fluid. The horizontal velocity in the puff front peaks at the same time the injection does (vertical dashed orange line) to later decay to values close to that typical of indoors conditions40 of 0.1 m s−1 indicated by the horizontal red dashed line.

FIG. 4.

Puff front temperature and vertical velocity.

FIG. 4.

Puff front temperature and vertical velocity.

Close modal

In the study by Bourouiba, Dehandschoewercker, and Bush,23 the volume of the puff front is defined as V p = η r 3 , where r is the longest semi-axis of the fitted ellipse projected in plane yz and η is a shape factor that takes the value of 4 3 π for spherical puffs. Equating Vp to the volume of an ellipsoid, V s = 4 π 3 σ 1 σ 2 σ 3 , η can be estimated as

η = 4 π 3 σ 1 σ 2 σ 3 r 3 .
(11)

The trajectory of the puff extracted from the temporal evolution of the ellipsoid centroid and the prediction from the model of Richards21 expressed in Eq. (5) are shown in Fig. 5(a) for two values of the entrainment coefficient. While the model prediction exhibits marginal puff deflection over the duration of the jet (end of cough indicated by the blue vertical dashed line), the numerical results suggest that deflection starts soon after the peak velocity has been reached (orange vertical dashed line). Note that the initial momentum (P0) in Eq. (5) has been calculated with the time averaged velocity w ¯ = w m / 2 . Contrary to the models that assume a constant velocity during the injection, the decrease in the inlet momentum, as injection velocity ramps down, allows buoyancy forces to start deflecting the jet after reaching the peak velocity. Once the cough has finished, the thermal puff continues to rise, while temperature mixes and turbulence decays due to dissipation. The deflection rate of the plume eventually decreases for c z 62 cm. This behavior cannot be reproduced by the model that predicts a continuous increase in the amount of entrained fluid into the puff.

FIG. 5.

Panel (a): Puff front trajectory ellipsoid centroid from the DNS (thin dotted black) and the Richards21 model for α = 0.20 (red) and α = 0.25 (blue). Panel (b): Temporal evolution of the centroid traveled distance from the DNS (thin dotted black), the Richards21 model for α = 0.20 (red) and α = 0.25 (blue), and the Scorer22 model (solid black). Vertical dashed orange and blue lines indicate the peak velocity and cough end times, respectively. Curves for t 1 / 2 and t 1 / 4 added for reference.

FIG. 5.

Panel (a): Puff front trajectory ellipsoid centroid from the DNS (thin dotted black) and the Richards21 model for α = 0.20 (red) and α = 0.25 (blue). Panel (b): Temporal evolution of the centroid traveled distance from the DNS (thin dotted black), the Richards21 model for α = 0.20 (red) and α = 0.25 (blue), and the Scorer22 model (solid black). Vertical dashed orange and blue lines indicate the peak velocity and cough end times, respectively. Curves for t 1 / 2 and t 1 / 4 added for reference.

Close modal

The temporal evolution of the thermal puff traveled distance s is shown in Fig. 5(b) that also includes the models of Richards21 given in Eq. (5) and Scorer22 given in Eqs. (6)–(9). While models predict a monotonic growth in s, the numerical results exhibit an acceleration stage during the laminar jet increase over 0 < t t m (orange vertical dashed line). Once the inlet velocity starts to ramp down, the s curve follows a t 1 / 2 [Eq. (6)] trend up to the end of the cough (blue vertical dashed line). The differences in the initial stages of the cough, explained by the initial laminar regime of the jet and the transient nature of the injection, decrease once the thermal puff stage is reached and s exhibits a t 1 / 4 growth rate [Eq. (9)]. While larger α values result in better predictions of the centroid trajectory, the faster rate of entrainment leads to underpredicted traveled distance s.

The three projection views of the best-fitting ellipsoid to the puff front at s = 5, 30, and 65 cm are shown in Fig. 6. Results suggest that, as the puff penetrates in the environment, the fastest growing semi-axis in this flow realization is the one quasi-parallel to the x direction.

FIG. 6.

The three projections on zy, zx, and xy planes of the best-fitting ellipsoid at three different values of traveled distance s = 5, 30, and 65 cm corresponding to times t 0.1, 0.24, and 1.0 s. All axes in cm.

FIG. 6.

The three projections on zy, zx, and xy planes of the best-fitting ellipsoid at three different values of traveled distance s = 5, 30, and 65 cm corresponding to times t 0.1, 0.24, and 1.0 s. All axes in cm.

Close modal

Details of the temporal evolution of the puff front topology are provided in Fig. 7. Results for the three semi-axes of the ellipsoid suggest that, after reaching the peak velocity, the puff topology starts to depart from the spheroidal shape assumed by Bourouiba, Dehandschoewercker, and Bush.23 The opposite thermal vertical structure across the puff leads to enhanced turbulent vertical mixing in the upper shear layer and diminished mixing in the lower one. As a result, the horizontal location of the top part of the front is delayed and the yz projected ellipsoid appears rotated in the counterclockwise direction for t > t c . Animations showing the temporal evolution of the θ ̃ = 0.025 isosurface in the top (zx projection), side (zy projection), and front (xy projection) views are included in Figs. 8–10, respectively (Multimedia view).

FIG. 7.

Temporal evolution of the 3D ellipsoid semi-axes σi (red, blue, and green) and the largest projected semi-axis on the zy plane σl (black) with local values of the entrainment coefficient α for the present DNS realization. Peak and end cough times are shown as vertical orange and blue dashed lines, respectively.

FIG. 7.

Temporal evolution of the 3D ellipsoid semi-axes σi (red, blue, and green) and the largest projected semi-axis on the zy plane σl (black) with local values of the entrainment coefficient α for the present DNS realization. Peak and end cough times are shown as vertical orange and blue dashed lines, respectively.

Close modal
FIG. 8.

Video 1: Top xz view showing the temporal evolution of the θ = 0.025 isosurface for the DNS of a mild cough. Time is in seconds. Multimedia view: http://dx.doi.org/10.1063/5.0042086.1

FIG. 8.

Video 1: Top xz view showing the temporal evolution of the θ = 0.025 isosurface for the DNS of a mild cough. Time is in seconds. Multimedia view: http://dx.doi.org/10.1063/5.0042086.1

Close modal
FIG. 9.

Video 2: Side xy view showing the temporal evolution of the θ = 0.025 isosurface for the DNS of a mild cough. Time is in seconds. Multimedia view: http://dx.doi.org/10.1063/5.0042086.2

FIG. 9.

Video 2: Side xy view showing the temporal evolution of the θ = 0.025 isosurface for the DNS of a mild cough. Time is in seconds. Multimedia view: http://dx.doi.org/10.1063/5.0042086.2

Close modal
FIG. 10.

Video 3: Front yz view showing the temporal evolution of the θ = 0.025 isosurface for the DNS of a mild cough. Time is in seconds. Multimedia view: http://dx.doi.org/10.1063/5.0042086.3

FIG. 10.

Video 3: Front yz view showing the temporal evolution of the θ = 0.025 isosurface for the DNS of a mild cough. Time is in seconds. Multimedia view: http://dx.doi.org/10.1063/5.0042086.3

Close modal

The growth rate of the puff radius r, defined as the longest semi-axis of the projected ellipsoid on the yz plane (black line), clearly exhibits different regimes over the duration of the cough. The results indicate that the decelerating jet exhibits a relatively small entrainment coefficient of α = 0.06 . Once the cough ceases, the puff seems to be characterized by intermittent episodes of rapid entrainment with coefficient values ranging over 0.18 < α < 0.25 . These events can also be seen in the measurements of Bourouiba, Dehandschoewercker, and Bush23 [see, for example, their Fig. 9(b)]. A value of α = 0.2 would, therefore, overpredict the cough growth during its early stage and underpredict it, at least locally, once the turbulent puff is totally developed.

The temporal evolution of the shape factor η defined in Eq. (11) is shown in Fig. 11. The pronounced disparity in the length of the three ellipsoid semi-axes suggests that the puff front shape notably differs from that previously assumed by both Richards20 and Bourouiba,23 which lead to underestimation of the puff volume growth rate. Numerical results suggest that after the peak velocity time, the initially slender puff rapidly grows until reaching its largest volume soon after the end of the cough. Near the end of the simulation at t 1.65 s, the puff volume shrinks to values of η approximately twice the value used in the study by Bourouiba, Dehandschoewercker, and Bush.23 

FIG. 11.

Temporal evolution of η as defined in Eq. (11) for the present DNS realization. Orange and blue vertical dashed lines indicate the peak and end of cough times, respectively. The value η = 3 suggested by Bourouiba, Dehandschoewercker, and Bush23 and η = 4 π / 3 are shown as a horizontal black and blue dashed lines, respectively.

FIG. 11.

Temporal evolution of η as defined in Eq. (11) for the present DNS realization. Orange and blue vertical dashed lines indicate the peak and end of cough times, respectively. The value η = 3 suggested by Bourouiba, Dehandschoewercker, and Bush23 and η = 4 π / 3 are shown as a horizontal black and blue dashed lines, respectively.

Close modal

Turbulence is a chaotic process in nature, and therefore, the sole DNS realization of the flow produced by a mild cough reported here does not provide information regarding turbulence properties. In addition, the flow configuration assumed a simplified respiratory track geometry, a stagnant and homogeneous environment, and an ambient air temperature typical of outdoors conditions at mid latitudes. Although changes in the “mouth” diameter and exit velocity would affect the Reynolds number, the variability of this parameter over the typical range observed for violent expiratory events including coughs and sneezes is expected to be relatively limited and, therefore, does not significantly change the integral quantities of interest, namely, the puff horizontal penetration and lateral growth rate. Regarding the ambient air temperature, a larger value of T closer to typical indoors conditions should result in weakened vertical puff deflection. With a vertical displacement of approximately 0.05 cm per horizontally traveled cm, an increase in the overall horizontal penetration in warmer ambient conditions is expected to be modest.

Current effort is directed to explore the effects of potentially key variables capable of significantly impacting the flow hydrodynamics of expiratory events. On the one hand, the rapid kinetic energy dissipation observed once the exhalation ceases in the present results suggests that, for times as short as t = 1 s after the cough onset, the puff average velocity decays to values typical of indoor conditions. This suggests that background currents produced by heating, ventilation, and air conditioning may play a significant role in the late stages of the thermal puff evolution. Similarly, ambient thermal stratification may notably change the turbulent dissipation by suppressing the vertical momentum and heat transport, which, in turn, could affect the horizontal, drag-dominated dispersion of pathogen-laden aerosols. Finally, milder exhalation events including talking or signing are suspected to exhibit significantly different dispersion patterns due to their sustained and transient nature.

The resolved numerical simulation of an idealized violent expiratory event provides unprecedented details on the temporal and spatial evolution of the jet and thermal puff produced by a mild cough. New insight into the flow hydrodynamics generated by a mild cough could help in developing new models for dispersion of pathogen-laden aerosols that are responsible for airborne transmission of infectious diseases. The DNS results presented here could be used to validate/improve reduced order models that can then be used in physics-based transmission-models.41 

While the numerical prediction of horizontal range agrees with theory, the trajectory of the puff front and the entrainment rate are found to differ significantly. Contrary to theoretical models and experiments, the DNS presented here accounts for the discontinuous injection that characterizes respiratory events where exhaled air accelerates until peaking and then decreases to the end of the cough. This laminar to turbulent transition over the jet-to-puff evolution, combined with the assumption of spherical puff topology and the turbulence role neglect under the linear momentum conservation premise, explain the shortcomings of the theoretical model in capturing the dynamics of thermal puffs produced by mildly violent respiratory events.

The authors report no conflict of interest.

This work was funded by Spanish Ministerio de Ciencia, Innovación y Universidades through Grants DPI2016-75791-C2-1-P and RTI2018-100907-A-I00 (MCIU/AEI/FEDER, UE) and also by the Generalitat de Catalunya through Grant 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 for testing and smoothening the computational mesh. We thank Professor Josep Anton Ferré comments and suggestions that undoubtedly helped to improve the manuscript.

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

1.
CSSE
, https://coronavirus.jhu.edu/map.html
for
COVID-19 Dashboard by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University
, The Johns Hopkins Coronavirus Resource Center,
2020
.
2.
H.
De-Leon
and
F.
Pederiva
, “P
article modeling of the spreading of coronavirus disease (COVID-19
,”
Phys. Fluids
32
,
087113
(
2020
).
3.
J. M.
Creeth
, “
Constituents of mucus and their separation
,”
Br. Med. Bull.
34
,
17
24
(
1978
).
4.
P. S.
Padmanabha
and
M. R.
Prasnna Simha
, “
Universal trends in human cough airflows at large distances
,”
Phys. Fluids
32
, 0
81905
(
2020
).
5.
K.
Nissen
,
J.
Krambrich
,
D.
Akaberi
,
T.
Hoffman
,
J.
Ling
,
A.
Lundkvist
,
L.
Svensson
, and
E.
Salaneck
, “
Long-distance airborne dispersal of SARS-CoV-2 in COVID-19 wards
,”
Sci. Rep.
10
,
19589
(
2020
).
6.
S.
Dutta
,
S.
Cheng
,
A.
Parikh
,
T.
Truscott
,
P.
Fischer
, and
L.
Chamorro
, “
On turbulence and particle transport in closed rooms
,” in Bulletin of the American Physical Society (
2020
).
7.
M.
Abuhegazy
,
K.
Talaat
,
O.
Anderoglu
, and
S. V.
Poroseva
, “
Numerical investigation of aerosol transport in a classroom with relevance to covid-19
,”
Phys. Fluids
32
,
103311
(
2020
).
8.
W.
Chen
,
N.
Zhang
,
J.
Wei
,
H.-L.
Yen
, and
Y.
Li
, “
Short-range airborne route dominates exposure of respiratory infection during close contact
,”
Build. Environ.
176
,
106859
(
2020
).
9.
Z.
Zhang
,
T.
Han
,
K. H.
Yoo
,
J.
Capecelatro
,
A. L.
Boehman
, and
K.
Maki
, “
Disease transmission through expiratory aerosols on an urban bus
,”
Phys. Fluids
33
,
015116
(
2021
).
10.
G.
Busco
,
S. R.
Yang
,
J.
Seo
, and
Y. A.
Hassan
, “
Sneezing and asymptomatic virus transmission
,”
Phys. Fluids
32
,
073309
(
2020
).
11.
M.
Abkarian
,
S.
Mendez
,
N.
Xue
,
F.
Yang
, and
H. A.
Stone
, “
Speech can produce jet-like transport relevant to asymptomatic spreading of virus
,”
Proc. Natl. Acad. Sci. U. S. A.
117
,
25237
25245
(
2020
).
12.
T.
Dbouk
and
D.
Drikakis
, “
On coughing and airborne droplet transmission to humans
,”
Phys. Fluids
32
,
053310
(
2020
).
13.
D.
Fontes
,
J.
Reyes
,
K.
Ahmed
, and
M.
Kinzel
, “
A study of fluid dynamics and human physiology factors driving droplet dispersion from a human sneeze
,”
Phys. Fluids
32
,
111904
(
2020
).
14.
T.
Dbouk
and
D.
Drikakis
, “
On respiratory droplets and face masks
,”
Phys. Fluids
32
,
063303
(
2020
).
15.
M.-R.
Pendar
and
J. C.
Páscoa
, “
Numerical modeling of the distribution of virus carrying saliva droplets during sneeze and cough
,”
Phys. Fluids
32
,
083305
(
2020
).
16.
B.
Wang
,
H.
Wu
, and
X.-F.
Wan
, “
Transport and fate of human expiratory droplets—A modeling approach
,”
Phys. Fluids
32
,
083307
(
2020
).
17.
E.
Renzi
and
A.
Clarke
, “
Life of a droplet: Buoyant vortex dynamics drives the fate of micro-particle expiratory ejecta
,”
Phys. Fluids
32
,
123301
(
2020
).
18.
R. S.
Scorer
, “
Experiments on convection of isolated masses of buoyant fluid
,”
J. Fluid Mech.
2
,
583
594
(
1957
).
19.
J. M.
Richards
, “
Experiments on the penetration of an interface by buoyant thermals
,”
J. Fluid Mech.
11
,
369
384
(
1961
).
20.
J. M.
Richards
, “
Puff motions in unstratified surroundings
,”
J. Fluid Mech.
21
,
97
(
1965
).
21.
J. M.
Richards
, “
Inclined buoyant puffs
,”
J. Fluid Mech.
32
,
681
692
(
1968
).
22.
R. S.
Scorer
,
Dynamics of Meteorology and Climate
, 1st ed. (
Wiley
,
1997
).
23.
L.
Bourouiba
,
E.
Dehandschoewercker
, and
J. W. M.
Bush
, “
Violent expiratory events: On coughing and sneezing
,”
J. Fluid Mech.
745
,
537
563
(
2014
).
24.
S.
Diwan
,
S.
Ravichandran
,
R.
Govindarajan
, and
R.
Narasimha
, “
Understanding transmission dynamics of COVID-19-type infections by direct numerical simulations of cough/sneeze flows
,”
Trans. Indian Natl. Acad. Eng.
5
,
255
261
(
2020
).
25.
J.
Gupta
,
C.-H.
Lin
, and
Q.
Chen
, “
Flow dynamics and characterization of a cough
,”
Indoor Air
19
,
517
525
(
2009
).
26.
D. D.
Gray
and
A.
Giorgini
, “
The validity of the Boussinesq approximation for liquids and gases
,”
Int. J. Heat Mass Transfer
19
,
545
551
(
1976
).
27.
E.
Mansour
,
R.
Vishinkin
,
S.
Rihet
,
W.
Saliba
,
F.
Fish
,
P.
Sarfati
, and
H.
Haick
, “
Measurement of temperature and relative humidity in exhaled breath
,”
Sens. Actuators, B
304
,
127371
(
2020
).
28.
P.
Fischer
,
J. W.
Lottes
, and
S. G.
Kerkemeier
, http://nek5000.mcs.anl.gov for nek5000 web page,
2008
.
29.
M.
Deville
,
P.
Fischer
, and
E.
Mund
, “
High-order methods for incompressible fluid flow
,”
Appl. Mech. Rev.
56
,
B43
(
2003
).
30.
K.
Mittal
,
S.
Dutta
, and
P.
Fischer
, “
Multirate timestepping for the incompressible Navier–Stokes equations in overlapping grids
,” arXiv:2003.00347 [physics.flu-dyn] (
2020
).
31.
A.
Fabregat
,
A.
Poje
,
T.
Özgökmen
, and
W.
Dewar
, “
Dynamics of multiphase turbulent plumes with hybrid buoyancy sources in stratified environments
,”
Phys. Fluids
28
,
095109
(
2016
).
32.
E.
Merzari
,
A.
Obabko
, and
P.
Fischer
, “
Spectral element methods for liquid metal reactors applications
,” arXiv:1711.09307 [cs.CE] (
2017
).
33.
R.
Vinuesa
,
P.
Negi
,
M.
Atzori
,
A.
Hanifi
,
D.
Henningson
, and
P.
Schlatter
, “
Turbulent boundary layers around wing sections up to rec = 1,000,000
,”
Int. J. Heat Fluid Flow
72
,
86
99
(
2018
).
34.
A.
Glezer
and
D.
Coles
, “
An experimental study of a turbulent vortex ring
,”
J. Fluid Mech.
211
,
243
283
(
1990
).
35.
N. R.
Panchapakesan
and
J. L.
Lumley
, “
Turbulence measurements in axisymmetric jets of air and helium. Part 1. Air jet
,”
J. Fluid Mech.
246
,
197
223
(
1993
).
36.
B. J.
Boersma
,
G.
Brethouwer
, and
F. T. M.
Nieuwstadt
, “
A numerical investigation on the effect of the inflow conditions on the self-similar region of a round jet
,”
Phys. Fluids
10
,
899
909
(
1998
).
37.
S.
Chaudhuri
,
S.
Basu
,
P.
Kabi
,
V. R.
Unni
, and
A.
Saha
, “
Modeling the role of respiratory droplets in Covid-19 type pandemics
,”
Phys. Fluids
32
,
063309
(
2020
).
38.
N.
Abani
and
R.
Reitz
, “
Unsteady turbulent round jets and vortex motion
,”
Phys. Fluids
19
,
125102
(
2007
).
39.
Yury
(
2020
). “Ellipsoid fit,” Matlab Central File Exchange. https://www.mathworks.com/matlabcentral/fileexchange/24693-ellipsoid-fit
40.
International Organization for Standardization
, ISO 7730 2005-11-15 Ergonomics of the Thermal Environment: Analytical Determination and Interpretation of Thermal Comfort Using Calculation of the PMV and PPD Indices and Local Thermal Comfort Criteria,
2005
.
41.
S.
Chaudhuri
,
S.
Basu
, and
A.
Saha
, “
Analyzing the dominant SARS-CoV-2 transmission routes toward an ab initio disease spread model
,”
Phys. Fluids
32
,
123306
(
2020
).