Towards realistic simulations of human cough: effect of droplet emission duration and spread angle

Human respiratory events, such as coughing and sneezing, play an important role in the host-to-host airborne transmission of diseases. Thus, there has been a substantial effort in understanding these processes: various analytical or numerical models have been developed to describe them, but their validity has not been fully assessed due to the difficulty of a direct comparison with real human exhalations. In this study, we report a unique comparison between datasets that have both detailed measurements of a real human cough using spirometer and particle tracking velocimetry, and direct numerical simulation at similar conditions. By examining the experimental data, we find that the injection velocity at the mouth is not uni-directional. Instead, the droplets are injected into various directions, with their trajectories forming a cone shape in space. Furthermore, we find that the period of droplet emissions is much shorter than that of the cough: experimental results indicate that the droplets with an initial diameter $\gtrsim 10\mu$m are emitted within the first 0.05 s, whereas the cough duration is closer to 1 s. These two features (the spread in the direction of injection velocity and the short duration of droplet emission) are incorporated into our direct numerical simulation, leading to an improved agreement with the experimental measurements. Thus, to have accurate representations of human expulsions in respiratory models, it is imperative to include parametrisation of these two features.


Introduction
Since the outbreaks of SARS-CoV in 2003 and SARS-CoV-2 in 2019, the role played by the turbulent multiphase flow in the transmission of infectious diseases via the airborne route has received increasing attention. The host-to-host transmission of respiratory disease is a complicated process with multiple stages including the exhalation, dispersion and inhalation of the pathogen (Bourouiba, 2020;Zhou and Zou, 2021;Bourouiba, 2021;Chong et al., 2021;Smith et al., 2020).In particular, a central piece to the puzzle is the dispersion of pathogencarrying droplets in turbulent flows.
A classic approach to mathematically predict the transmission of an infectious disease among populations is the compartmental models (Tolles and Luong, 2020). Developed in 1927, the most basic compartmental model is an SIR model (Kermack and McKendrick, 1927), where the population is separated into three so-called compartments, namely 'Susceptible', 'Infected' and 'Removed'. Individuals can transfer between compartments based on ordinary differential equations while the total population remains constant. The infection rate of susceptible individuals can be improved with more physical insight via dose-response models (Brouwer et al., 2017), or the Wells-Riley model in particular for dis-2 eases transmitted via the airborne route (Noakes et al., 2006). The application of these two groups of models and a detailed comparison are presented in the review paper by Sze To and Chao (2010). In essence, the Wells-Riley model centers on the concept of 'quantum', which is the dose required to start an infection, while the dose-response model uses the amount of pathogen taken in by the susceptible individual. In both approaches, the complicated transmission process that involves a series of stages (exhalation, dispersion, ventilation, inhalation, etc) is parametrised based on some physical assumptions or empirical observations (Lelieveld et al., 2020;Bazant and Bush, 2021;Jones et al., 2021;Nordsiek et al., 2021). As a typical example, the room where the infectious disease is spread from is assumed to be 'well-mixed', i.e. the pathogen distributes uniformly in the room. More recently, some efforts have been made to incorporate the temporal and spatial inhomogeneity in the model (Mittal et al., 2020a;Yang et al., 2020). A sophisticated software tool has been developed to estimate the infection risk under a number of scenarios, including various modes of ventilation and different levels of activities (Mikszewski et al., 2020). These models are relatively easy to use and hence highly appealing when developing social guidelines based on the risk of transmission at a wide range of real-life scenarios. However, the derivation of these models, especially the selection of the parameters hinges on the detailed understanding of the governing physics in each step of the transmission. Alongside the effort in the epidemiology community, there has also been an ongoing endeavour in understanding and modelling the flow physics in the transmission of respiratory diseases. In the pioneering work of Wells (1934), the lifespan of a droplet is illustrated through a simple 'evaporation-falling curve', where the fate of a small droplet is considered to be full evaporation, while that of a large droplet is ground deposition. This simplistic model was then improved by Xie et al. (2007), where the salinity of the droplet, the ambient temperature and relative humidity conditions and the buoyant respiratory jet were taken into consideration. The effect of ambient conditions on the droplet evaporation is then incorporated into one parameter, the effective evaporation diffusivity proposed by Balachandar et al. (2020). The coupling between the carrying fluid and the disperse phase is further incorporated into the puff model by Bourouiba et al. (2014), and a good agreement with the laboratory measurement of a multiphase puff is observed, but a direct comparison with clinical data, such as the speed and the trajectory of the puff is still lacking.
There is now a consensus that the flow physics in the transmission of respiratory diseases involve a multiphase flow (Bourouiba, 2020). Indeed, various human respiratory events, such as breathing, speaking, coughing and sneezing, can be simplified as transient turbulent jets with liquid droplets as the second phase (Balachandar et al., 2020;Mittal et al., 2020b;Bourouiba et al., 2014;Bourouiba, 2020;Ng et al., 2021;Chong et al., 2021;Abkarian et al., 2020;Stadnytskyi et al., 2020). The flow field associated with a single respiratory event have been investigated through both experiments with human subjects (VanSciver et al., 2011;Bahl et al., 2020) and laboratory or numerical simulations (Bourouiba et al., 2014;Wei and Li, 2017;Chang et al., 2020;Abkarian et al., 2020;Ng et al., 2021;Chong et al., 2021;Fabregat et al., 2021). The general picture of a cough or sneeze is that the droplet-laden, warm and moist air is injected at a high speed over a short duration ( 1 s), which propagates under the effect of both the initial momentum and buoyancy (Bourouiba et al., 2014;Ng et al., 2021), while in continuous events such as speaking, a 'puff train' assimilates to a turbulent jet like structure in the far field . The fate of the 4 droplets varies largely depending on their size : large droplets with a diameter ∼ O(1000µm) deposit on the ground in a near ballistic motion, while small droplets with a diameter ∼ O(10µm) are trapped in the humid puff and have prolonged lifetimes compared to the prediction of Wells (1934). Condensation is observed at low ambient temperature and high relative humidity settings . Similar to the need for judicious assessment of the accuracy of theoretical models with clinical data, a direct comparison of numerical simulations with real respiratory events is imperative to test and elevate the accuracy of these models.
Accordingly, in this study, we report a unique pair of datasets of human coughs with both experimental and numerical components with similar conditions. The numerical dataset is obtained using direct numerical simulations (DNS) based on the methods employed by Chong et al. (2021) and Ng et al. (2021), and the experimental component is the extension of the method used by Bahl et al. (2020) to coughs. We reveal that both the droplet emission duration and the initial spread angle have strong influence on the far-field behaviour of the flow field. These two refinements to the flow inlet conditions are validated through the experimental data collected from a human subject in §2. Two numerical simulation cases are formulated based on these refined conditions in §3, and their detailed effect on the simulated flow is discussed in §4.

Experimental observations
To capture the flow dynamics of both the airflow and the droplets expelled during a cough from human exhalations, we utilise two sets of experiments.
Namely, a spirometer is employed in the subject's mouth to capture the volumetric flow rate expelled during a cough ( Fig. 1(b)), and Particle Tracking Velocimetry (PTV) experiments quantify the motion of the expelled droplets.
For the latter, a volume illumination configuration is employed using a high In the present study and Gupta et al. (2009), the flow rate is measured experimentally by a spirometer, whereas the inlet velocity profile used by Chong et al. (2021) is converted to flow rate by assuming a circular mouth shape with a constant diameter D = 23 mm. The scatter of the peak location and magnitude from 25 subjects of Gupta et al. (2009) is shown by the horizontal and vertical error bars at the peak flow rate, and the variation of the cough duration is shown by the horizontal error bar centred around t = 0.5 s. The three black curves (rep1, rep2 and rep3) are three repeating measurements of the same subject.
light from the droplets expelled during a cough. The high-speed camera was This arrangement provided a resolution of 312.5 µm/pixel and a depth of field of 40 mm. Respiratory droplets expelled by the subject were used as tracer particles, and no additional seeding was introduced. The basic setup is depicted in Fig. 1(a). We note that all the experiments were performed at an ambient temperature of 22 • C.
To perform precise particle tracking, the high-speed image sequences are first pre-processed to minimise background/sensor noise and isolate head movement. After that, PTV was performed on the image sequence using Lavision® Davis 8.4. Here an initial Particle Image Velocimetry (PIV) pass on the image sequence was employed to obtain a velocity field estimate and then individual droplets were tracked using PTV algorithm (Cowen and Monismith, 1997). Further details on the experimental setup and processing can be found in Bahl et al. (2020).
In order to draw a direct comparison of the flow rates in our experiments,   Bourouiba et al. (2014) and the inline holography in Shao et al. (2021). In the atomisation process of a high-speed gas jet shearing a liquid film, which is considered as a major mechanism of droplet generation in respiratory tracts (Johnson and Morawska, 2009), there is also experimental evidence that droplets are ejected at an angle of 20 • − 40 • relative to the jet flow direction (Descamps et al., 2008). Note that the spreading angle of the initial velocity field discussed here is a near field quantity (O(D) where D is the mouth diameter), which is different from the far field (> O(5D)) spreading angle of the respiratory jet typically reported to be 10 • − 15 • with respect to the centre-line (Gupta et al., 2009;Abkarian et al., 2020;Dudalski et al., 2020).
By visually inspecting the raw PTV images, the droplet emission rate is found to be the highest at the beginning of a cough, then it quickly diminishes long before the flow injection velocity at the mouth reaches zero (de Silva et al., 2021). In other words, the droplet emission is concentrated at the onset of the cough for a much shorter duration than that of the flow injection. To provide further quantitative evidence regarding this observation, we utilise the numerical dataset in Chong et al. (2021), and compare the number of droplets in a region Note that the PTV method in this study can confidently detect droplets with d 10µm although smaller droplets are also likely to be detected. Therefore, the emission cutoff time discussed above is only strictly applicable to droplets in To summarise, two major features of a real human cough can be concluded from the experiments: first, the initial velocities of droplets are not uni-directional.
Instead, droplets with a relatively high momentum spread into various directions after exiting the mouth, as evidenced by their trajectories which form a cone shape with an apex around 90 • . Second, the majority of droplets with d 10µm are emitted between 0-0.05 s, which is less than 10% of the duration of a cough based on a nonzero injection flow rate. Incorporating these features in numerical simulations are therefore mandatory in order to achieve a better representation of realistic human respiratory events.

Numerical investigations
In this section, we provide details of the numerical cases designed to study the effect of the spread angle and droplet emission cutoff time.
We consider an incompressible fluid of gas phase, with both temperature and vapour concentrations coupled to the velocity field by employing the Boussinesq approximation. The gas phase is solved using DNS by a staggered second-order accurate finite difference scheme and marched in time using a fractional-step third-order Runge-Kutta approach (Verzicco and Orlandi, 1996;. Buoyancy effect of the puff as a result of the temperature and vapour mass fraction difference between the exhalation and the ambient environment is considered in the simulation. For droplets, we apply the spherical point-particle The mouth shape in both cases is assumed to be a circle with diameter D = 23 mm, which remains constant in time. By integrating the streamwise velocity across the mouth area, the flow rate of the exhalation can then be expressed asṁ andṁ where ρ is the air density, G(r) is the Gaussian function and r is the variable of integration with r = 0 at the centre-line of the mouth. φ(r) is the angle between the injection velocity vector and the horizontal axis in the θ = 45 • case, which can be related to r as The injected flow rate and forward momentum govern the propagation of the respiratory jet. For a given flow rate, the injection velocity magnitudes, U 0 • and U 45 • , can then be determined using (1). Here we assume that both cases  Since the droplet distribution in the flow field is relatively sparse, the dropletto-droplet interaction is negligible. For the same reason the modification to evaporation rates due to droplet clustering is also insignificant. Therefore, it is possible to remove droplets emitted after a certain time in post processing to study the effect of the droplet emission cutoff time.

Results and discussion
The spatial distribution of the droplets in the DNS cases at a range of time instances from 0.01 s to 0.2 s are shown in Fig. 6. in the experiment (Fig. 7a), with droplets scattered both above and below z m .
In the θ = 0 • case (Fig. 7b), however, very few droplets can be found above the mouth location z m . The build-up of droplets near the mouth is absent, and the droplet distribution is almost uniform in the streamwise direction. wise deceleration is complete before they reach the floor. Comparing Fig. 8(c) and (d ), one would notice that droplets in the θ = 45 • case settle closer to the mouth. Instead of carried downstream by the jet with a similar velocity, droplets emitted with a large spread angle penetrate deeper sideways into the still ambient fluid, so they experience larger relative velocity (therefore higher aerodynamic drag and stronger deceleration) compared to the θ = 0 • case. The reduction in the settling distance is not limited to droplets emitted with a high wall-normal velocity component: it is also seen in those with small or zero wallnormal component, as a result of the lower streamwise velocity in the widened jet.
The settling behaviour of droplets is further explored by examining the distribution of deposition locations in the ground plane of z = 0 in Fig. 9. The location where a droplet reaches the ground is represented by a dot in the x − y plane, with the colour determined by its initial diameter d 0 . The contaminated region on the ground for each size class is shown by the convex hull. Although the mouth is at z ≈ 0.18 m in the simulation domain, for smaller droplets (d 0 100 µm), the streamwise velocity has already decreased to 0 before reaching the z = 0 plane. Therefore, the deposition pattern reported in Fig. 9 would still be applicable when the mouth is at z ≈ 1.5 m above the ground, which is the case of a standing person.
Comparing Fig. 9(a) and (b), it is apparent that the deposition location scatters over a wider range in the spanwise direction for θ = 45 • , while it is mostly concentrated along the centre-line in the θ = 0 • case. This is a progress towards a more realistic cough, as evidenced by the improved agreement with the data of human subjects   (Fig. 10). In Fig. 9(b), the spanwise width of the convex hull increases with increasing d 0 , presumably because larger droplets have a higher initial momentum and they can continue for a longer distance along the direction of their initial velocity. The streamwise extent of the convex hull is shorter in the θ = 45 • case compared to the θ = 0 • case, indicating that droplets (especially those with d 0 100 µm) have a strong tendency to deposit closer to the mouth.
We note that the initial size distribution of droplets exhaled during respiratory events is presently actively investigated. The simulation results can be affected by the size distribution: for instance, if we use the distribution of , which has more small (O(10) µm) droplets, a higher number of trajectories can be expected in Fig. 8(e-f ). However, here our focus is on capturing the dominant effects of an actual cough.
The scalar fields of relative humidity (RH) are shown in Fig. 11. At t 1 = 0.15 s, despite the higher injected streamwise momentum, the θ = 45 • case has a shorter and wider turbulent jet compared to the θ = 0 • case. The flow also appears more turbulent close to the mouth (x 0.2 m) in the former. The cone angle α of the jet, which can be determined from tan α = R jet (x)/x, where R jet (x) is the jet radius at the streamwise location x, is very different from the spread angle θ. Comparing the white reference lines in the figure with the scalar field of the jet, it is shown that the jet cone angle is approximately 10 • in the θ = 45 • case, while it is only slightly smaller in the θ = 0 • case. Therefore, a 45 • change in the spread angle θ only results in a less than 5 • change in the cone angle α.
At t 2 = 0.6 s ( Fig. 11c and d ), a vortex ring separated from the main jet is observed in both cases. These vortices are formed at the front of the jet, which can be seen attached to the jet at t = t 1 . As the injection velocity decreases and the jet front decelerates, the self-induced velocity of the vortex ring eventually exceeds the velocity of the jet front, then the vortex ring detaches from the jet.
Such vortex rings have also been observed in the large-eddy simulation results by Liu et al. (2021), and their formation, strength and velocity are shown to be very sensitive to the initial conditions of the simulation.

Conclusions and outlook
In this study, we present a unique comparison of both experimental and numerical datasets on a human cough with matched conditions. Specifically, we implemented two important features of the human cough with insight gained from the comparison with the experimental data: the injection of the coughing flow which is shown to be better described with a spread angle of 45 • , and the emission of large droplets only occurs for a short period of time (≈ 0.05 s) at the beginning of a cough, whose length is much longer (≈ 1 s). The contribution of these features to the simulation result is examined by comparing two DNS cases with the conventional and with the modified inlet conditions based on the aforementioned experimental observations from a human cough. Overall, a significant improvement in the agreement with the experimental data is observed after implementing the modifications. We further demonstrate that by considering the 45 • initial spread angle, the simulated jet appears thicker and shorter.
Droplets deposit on the ground closer to the mouth in the streamwise direction, but scatter over a wider range in the spanwise direction, which is consistent with the experimental observations. The aim of future work must be to extend the present comparative approach between experiments and numerical simulations to other respiratory events such as singing and in particular speaking, which accumulatively plays an even more important role in the aerosol release Yang et al., 2020). Larger files such as the scalar fields will be available upon request. The experimental data cannot be directly shared due to ethics limitations.

Acknowledgements
by employing the Boussinesq approximation. The motion of the gas phase is assumed to be incompressible ∂u i /∂x i = 0, and governed by the momentum, energy, and vapour mass fraction equations: The last term on the right-hand side of equation (1) represents the coupling of temperature and vapour mass fraction to the momentum equation.
For droplets, the spherical point-particle model is applied, and we consider the conservation of momentum (Maxey-Riley equation [2]), energy, and mass as follows: du i,n dt = (β + 1) Du i,g,n Dt + (β + 1) 3ν air (u i,g,n − u i,n ) r 2 n + gβê y , The notations we used here are as follows: u i , u i,n , and u i,g,n are the velocity of gas, droplets, and gas at the location of droplets, respectively. Similarly θ g , θ n , and θ g,n are in Kelvin and used to represent the temperature of gas, droplets, and gas at the location of droplets, respectively.
c is the vapour mass fraction, r n the droplet radius, A n the surface area of the droplets, V n the volume of the droplets, m n the mass of the droplets. Also p denotes the reduced pressure. h m is the heat transfer coefficient. L is the latent heat of vaporisation of the liquid. ρ l and c p,l are the density and specific heat capacity of the droplets assumed to consist of water. ν air is the kinematic viscosity of air. k g is the thermal conductivity of gas which relates to the thermal diffusivity of gas D g by k g = D g ρ g c p,g . ρ g and c p,g are the density and specific heat capacity of gas, Note that ρ g c p,g = (ρ a c p,a + ρ v c p,v ), with ρ a and ρ v being the densities of air and vapour and c p,a and c p,v being the specific heat capacities of air and vapour. c drop and c f luid denote the vapour mass fractions of the droplet and the fluid at the location of the droplet. β is the density ratio term defined as β = 3ρ g /(ρ g + 2ρ l ) − 1.
Equation (5) and equation (6) are closed using the Ranz-Marshall correlations [3], which are reasonable assumptions for moving point droplets in our study since the droplet Reynolds number are at most O(100), in accordance with the applicable limits [4]. The correlations give us the estimation of h m and Sh drop for a single spherical droplet.
h m r/(D g c p,g ρ g ) = 2 + 0.6Re where we have droplet Reynolds number Note that, by definition, Sh drop and Re drop are also functions of r. The realistic parameters needed are listed in Table (1). The relative humidity at the location of the droplets is defined as P vap,drop /P sat = c/c sat,vap , where the vapour mass fraction c is solved from equation (3). Therefore, to calculate the local relative humidity, the saturated vapour mass fraction c sat,vap is first determined by the ideal gas law: where P sat and R are the saturated vapour pressure, and specific gas constant of water vapour.

Coupling methodology
To interpolate gas quantities to the droplet locations, we employ the tri-cubic Hermitian interpolation scheme, which are sufficiently accurate for turbulent flows and comparable to B-spline interpolations [5,6]. The backwards forcing of the droplets onto the gas phase uses the tri-linear projection onto the eight nearest nodes to the droplets location [7].

Numerical setup
To numerically solve the equations, we used our finite difference solver AFiD [8] with highperformance Message Passing Interface (MPI) and point-particle model. The size of the computational domain in dimensional form is 0.37 m (spanwise length) ×0.37 m (height) ×1.10 m (streamwise length) and is tested to be large enough to capture the cough vapour and spreading droplets. The grid points chosen are 512 × 512 × 1536 to ensure that enough resolution has been employed. The mouth is modelled as a circular inlet centred at mid-height of the domain.

Initial droplet distribution
The topic of distribution of the initial droplet sizes is in itself a subject of considerable importance and active debate [9,10,11,12]. Here, we seeded the respiratory event with droplets with initial diameters ranging between 10 µm and roughly 1000 µm, based on an experimental measurement [10]. Estimates of the droplet volume fractions give O(10 −5 ) and so droplets are assumed not to collide or coalesce. These droplets also do not set the relative humidity of the puff, since the relative humidity of the puff at the inflow is assumed to be 100%. Note that even smaller droplets could be added, but they would further extend the required CPU time and are not necessary to convey the main message of this work.

Cough properties
The cough temporal profile we apply is the average of the spirometer measurements shown in figure 2 in the manuscript. We did not impose any noise on the inflow velocity, and the transition to turbulence was not forced, but occurred naturally in the simulation. The inlet velocity, temperature and vapour mass fraction spatial profiles are prescribed as a Gaussian profile centred at the mouth location. Note that we set saturated vapour mass fraction at the inlet, which is calculated based on temperature of the human respiratory tract [13]. For the initial droplet size, we employ a similar distribution as [10] with 2500 droplets. The droplets are randomly positioned at the inlet and evenly injected in time with the velocity matching the local inlet velocity.

Initial and boundary conditions
The flow field is initialised as quiescent with uniform atmospheric pressure, temperature (20 • C) and a relative humidity of 50%. Neumann condition is applied to the boundary of the simulation domain, except for the mouth opening where the injection velocity is given. When a droplet reaches the boundary of the domain, it is flagged as 'dead' and then removed from the simulation. The simulation of a droplet also terminates when its diameter falls below 4.6µm, as simulating smaller droplets require significantly more CPU time without altering the main results.