Direct numerical simulation of the turbulent flow generated during a violent expiratory event

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 × 10⁶ people worldwide killing more than 1.7 × 10⁶ 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.⁵ 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 experiments^10,11 and numerical simulations^12,13 to analyze various aspects of the jet flow and the accompanying aerosol transport for different expiratory events including sneezing,¹⁰ coughing,¹² talking, and breathing.¹¹ 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 wind¹² and masks¹⁴ 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 cough¹² 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⁻¹. 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.¹⁰ 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.¹³ 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 events¹⁵ 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.¹⁶ Renzi and Clarke showed that the dynamics of the expiratory jet could be modeled by extending the theory of buoyant vortex rings.¹⁷ 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.

Scorer¹⁸ 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 estimated¹⁸ to be η = 3.

Puffs, occurring when the isolated mass of fluid has some initial momentum, were investigated by Richards¹⁹ 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 Richards²⁰ 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 β₀ with respect to gravity aligned with the y direction, Richards²¹ 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, Scorer²² 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,²¹ Bourouiba, Dehandschoewercker, and Bush²³ 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 Scorer²² and Bourouiba, Dehandschoewercker, and Bush²³ 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,²⁴ 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 captured²⁵ 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⁻¹ at $t m = 0.15$ s and then ramps down linearly to zero at $t c = 0.4$ s, i.e.,

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 modelled²⁵ 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.²⁶ 

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

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 %$⁠,²⁷  $ρ = 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$⁠.

The results presented here have been obtained using Nek5000,²⁸ 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.²⁹ 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 Mittal³⁰) 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$ × 10⁶. 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 Coles³⁴ 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 Lumley³⁵ (see, for example, the study by Boersma, Brethouwer, and Nieuwstadt³⁶) 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 Coles³⁴ and Panchapakesan and Lumley,³⁵ 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⁻¹) 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.

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

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 V₀ 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 written²³ as $r ( t ) = α s ( t )$⁠.

Previously reported^18,19 values of α for thermals are 0.25 and between 0.13 and 0.53. In their experiments, Bourouiba, Dehandschoewercker, and Bush²³ 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 z_j, the axial velocity W_j and the radial spread R_j are³⁸ 

For the puff stage (⁠ $t ≥ t j , end$⁠), the axial location z_p can be expressed as²² 

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,³⁸  $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,²²  $θ j ≈ 11.3 °$ and n_j = 5 and for puffs,²²  $θ p ≈ 14.0 °$ and n_p = 4.

The trajectory followed by the thermal puff is obtained following an analogous approach to that used by Bourouiba, Dehandschoewercker, and Bush²³ 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 i_k with a prescribed tolerance of $ε = 10 − 3$⁠, i.e.,

The instantaneous temperature integrated along x, $⟨ θ ̃ ⟩ x$⁠, is thought to mimic the photographic images taken during the experiments by Bourouiba, Dehandschoewercker, and Bush.²³ 

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

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 conditions⁴⁰ of 0.1 m s⁻¹ indicated by the horizontal red dashed line.

In the study by Bourouiba, Dehandschoewercker, and Bush,²³ 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 y – z and η is a shape factor that takes the value of $4 3 π$ for spherical puffs. Equating V_p to the volume of an ellipsoid, $V s = 4 π 3 σ 1 σ 2 σ 3$⁠, η can be estimated as

The trajectory of the puff extracted from the temporal evolution of the ellipsoid centroid and the prediction from the model of Richards²¹ 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 (P₀) 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.

The temporal evolution of the thermal puff traveled distance s is shown in Fig. 5(b) that also includes the models of Richards²¹ given in Eq. (5) and Scorer²² 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.

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.²³ 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 y – z 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 (z – x projection), side (z – y projection), and front (x – y projection) views are included in Figs. 8–10, respectively (Multimedia view).

The growth rate of the puff radius r, defined as the longest semi-axis of the projected ellipsoid on the y – z 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 Bush²³ [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 Richards²⁰ and Bourouiba,²³ 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.²³ 

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.⁴¹ 

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.