key: cord-0541058-xi1ltlq2
authors: Monroe, K.; Yao, Y.; Lattanzi, A.; Raghav, V.; Engineering, J. Capecelatro Department of Aerospace; Michigan, University of; Arbor, Ann; MI,; Engineering, Department of Mechanical; Engineering, Department of Aerospace; University, Auburn; Auburn,; AL,
title: Role of pulsatility on particle dispersion in expiratory flows
date: 2021-02-28
journal: nan
DOI: nan
sha: 7aff4f3148ad266d3065729005ffe2da63a6af6c
doc_id: 541058
cord_uid: xi1ltlq2

Expiratory events, such as coughs, are often pulsatile in nature and result in vortical flow structures that transport respiratory particles. In this work, direct numerical simulation (DNS) of turbulent pulsatile jets, coupled with Lagrangian particle tracking of micron-sized droplets, is performed to investigate the role of secondary and tertiary expulsions on particle dispersion and penetration. Fully-developed turbulence obtained from DNS of a turbulent pipe flow is provided at the jet orifice. The volumetric flow rate at the orifice is modulated in time according to a damped sine wave; thereby allowing for control of the number of pulses, duration, and peak amplitude. The resulting vortex structures are analyzed for single-, two-, and three-pulse jets. The evolution of the particle cloud is then compared to existing single-pulse models. Particle dispersion and penetration of the entire cloud is found to be hindered by increased pulsatility. However, the penetration of particles emanating from a secondary or tertiary expulsion are enhanced due to acceleration downstream by vortex structures.

Particle transport by turbulent free-shear jets plays a crucial role in many engineering and environmental applications. For example, atomization of liquid fuels leads to complex droplet size distributions and dispersion patterns that strongly influence internal combustion engine efficiency 1 while pyroclastic density currents, generated by large density differences between gas-particle mixtures, feeds explosive volcanic eruptions 2 . Of particular importance during the COVID-19 pandemic is the transmission of liquid droplets and aerosols (referred to interchangeably as particles herein) due to coughing, sneezing, or continuous speech 3 . Accurately describing particle dispersion from expiratory events is a critical aspect to defining physics-informed guidelines for social distancing best practices. While remarkable insight has been gained from analytical 4,5 , experimental 6 , and computational [7] [8] [9] works, the vast majority of studies are restricted to single expulsion events. However, realistic coughing is often characterized by multiple expulsions that lead to vortex-vortex interactions, which can have significant consequences on particle dynamics 10 (see Fig. 1 ). Experimental measurements have demonstrated that realistic coughs are pulsatile, involving a sequence of coughing events, sometimes referred to as "cough epochs." 12, 13 The flow rate associ-ated with a typical human cough is shown in Fig. 1 . Multiple pulses are observed over a duration of approximately 1 s, with the peak amplitude occurring at 0.1 s 11 . Gupta et al. 14 experimentally characterized the flow dynamics of coughs from human subjects, showing that the flow rate variation of a cough with time can be defined as a combination of gamma probability functions.

While single-pulse expiratory events have been well-studied, the influence of pulsatility on both the particle and fluid physics has received significantly less attention.

To understand the influence of pulsatility, it is important to first consider the modes of particle generation and sites of origination 15 during speech, coughing, or sneezing -Bronchiolar (droplet size < 1 − 5 µm); Laryngeal (5 − 20 µm); and Oral (> 50 µm) 16 . If the site of severe infection is deeper in the lungs (bronchiolar), the expelled aerosols/droplets are generated by a "fluid-film burst" mechanism 17 during collapse and reopening of the small airways resulting in small size particles. Since several respiratory infections, including H5N1 and SARS-CoV, replicate primarily in the bronchioles and alevoli 15, 18 , the aerosols/droplets generated in the lower airways are likely to contain higher doses of virus particles. In this case, secondary and tertiary pulses of a multi-pulse cough will expel the volume of air originating predominantly from deeper within the lungs and is expected to contain a higher concentration of virus particles. As such, based on the sites of severe infection in the respiratory tract, we hypothesize that the volume of air expelled by secondary and tertiary pulses could contain a higher viral load (illustrated in Fig. 1 ) and that the resulting vortexvortex interactions could significantly influence the dispersion of these more infectious particles.

It is now well established that interactions between turbulence and particles can give rise to preferential concentration, which describes the accumulation of particles away from highly vortical regions of the turbulent flow [19] [20] [21] [22] [23] . When the Stokes number, defined as the ratio of particle-tofluid time scales, is near unity, particles are directed by coherent vortical structures to create nonhomogeneities in concentration and the onset of clusters. Large-scale velocity gradients present in free-shear flows affect the transport of small (Kolmogorov-scale) heavy particles and the clustering process at small scales 24, 25 . Gualtieri et al. 26 showed that free shear flows generate anisotropic velocity fluctuations which, in turn, arrange particles in directionally biased clusters. In the presence of gravity, preferential concentration by turbulence has been observed to cause particles to further accumulate near the downward moving side of vortices, referred to as preferential sweeping [27] [28] [29] .

The gravitational settling of aerosol particles can be enhanced by this mechanism by as much as 50% 28 .

Turbulent transport in statistically stationary, axisymmetric, free jets has been well character-ized experimentally [30] [31] [32] and numerically [33] [34] [35] . Chein and Chung 10 demonstrated that particles with relatively small Stokes numbers disperse laterally at approximately the same rate as fluid particles, while particles with larger Stokes numbers exhibit significantly less dispersion. In particular, particles with intermediate Stokes numbers are transported laterally farther than fluid particles due to enhanced entrainment by vortex structures. Shortly after, Longmire and Eaton 36 showed that particles become clustered in the saddle regions downstream of vortex rings and are propelled away from the jet axis by the outwardly moving flow. More recently, direct numerical simulations (DNS) of particle-laden round jets by Li et al. 37 showed that all particles, regardless of their size, tend to preferentially accumulate in regions with larger-than-mean fluid streamwise velocity.

Particle dispersion was found to be directly associated with three-dimensional vortex structures.

While incredibly valuable, the aforementioned studies are restricted to jets with inflow characteristics that remain constant in time. By contrast, the transient characteristics of turbulent pulsatile jets are far less understood.

In this work, a realistic human cough is investigated computationally through DNS of pulsatile, turbulent, particle-laden jets. Fully-developed turbulence is provided at the orifice exit (mouth) using data obtained from an auxiliary simulation of turbulent pipe flow. The flow rate of the incoming turbulence is modulated in time according to a prescribed profile that controls the number of pulses, its duration, and peak amplitude. Particles are seeded in the flow with diameters sampled from a lognormal distribution informed by experimental measurements from the literature.

Two-phase statistics, in particular fluid entrainment and particle evolution, are then reported for each case, with emphasis on the effect of pulsatility on the resulting vortex structures and particle dispersion.

The present work considers a three-dimensional pulsatile jet laden with liquid droplets expelled into an ambient surrounding. Particles are considered to be well characterized as water droplets, and thus their density is held constant ρ p = 998 kg/m 3 . The fluid is considered to be air with a density of ρ = 1.172 kg/m 3 and kinematic viscosity of ν = 1.62 × 10 −5 m 2 /s. The diameter of the orifice exit (mouth) is taken to be D = 0.02 m. A Cartesian domain with length in the x (streamwise), y (spanwise, gravity-aligned) and z (spanwise) directions are L x = 40D, L y = 20D, and L z = 20D, respectively (see Fig. 2 ). The domain is discretized using N x = 1024 and N y = N z = 420 grid points, with exponential grid stretching in the yand z-directions. The spanwise grid spacing varies between 4.98 × 10 −4 m ≤ ∆y, ∆z ≤ 2.1 × 10 −3 m such that the minimum grid spacing at the jet centerline is D/40. Previous work has shown this level of resolution is sufficient for free-shear jets at similar Reynolds numbers 37 . A Dirichlet boundary condition is enforced at the jet inlet, a convective outflow is enforced at the downstream boundary, and all other boundaries are treated as slip walls. To prevent fluid recirculation within the computational domain, a co-flow is introduced along the positive x-direction with velocity magnitude 0.32 m/s. The co-flow is ∼ 7% of the peak inflow velocity U 0 and was observed to have negligible effect on the particle dynamics. 

Fully-developed turbulence is fed into the jet orifice using an auxiliary simulation of a turbulent pipe flow. The auxiliary simulation was performed using 256 grid points across the diameter with a bulk velocity of U 0 = 4 m/s (a typical peak velocity associated with expiratory events 6, 38 ), corresponding to a bulk Reynolds number Re b = U 0 D/ν = 4938. Further details on the pipe flow simulation are provided in Appendix A. Here, we note that the turbulent pipe simulation is statistically stationary and evolves to a constant bulk velocity U 0 defined by an imposed pressure gradient. To obtain a pulsatile turbulent inflow in the main simulation, the fluid velocity at the jet inlet u u u(x = 0, y, z;t) is adjusted dynamically to control the volumetric flow rate in time. Building upon the experimental observations of Gupta et al. 14 , we propose a self-similar profile for the bulk velocity resulting from multiple expulsions.

The proposed functional form for the pulsatile volumetric flow rate Q(t) is given by a damped sine wave according to

where Q 0 = U 0 A 0 and A 0 = πD 2 /4 is the area of the orifice exit. The fluid velocity is read in from the auxiliary pipe flow simulation and rescaled such that u u u(x = 0, y, z,t) · n n n dydz = Q(t), where n n n = [1, 0, 0] T is the outward surface normal. In the present study, we consider three profiles corresponding to one, two, and three pulses (see Fig. 3 ). The relaxation time is chosen to be τ = 

To accurately characterize the particle size distribution generated by coughing, we employ a lognormal distribution fit to the experimental measurements of Duguid 39 The turbulence Stokes number, St η = τ p /τ η , may be utilized to gauge the role of particle inertia, where τ p = ρ p d 2 p /(18ρν) is the particle response time, τ η = (ν/ε) 1/2 is the Kolmogorov time scale of the fluid phase, and ε is the turbulence dissipation rate. When analyzing the results in Sec. III, particles are demarcated into three size ranges: d p ∈ [1, 30] , [30, 60] 16] . We note that the Stokes numbers are defined using values taken at the orifice exit and therefore will decay as particles evolve in time and downstream. Nevertheless, these Stokes number ranges provide insight into the relative inertia of the fluid and particle phases.

Specifically, particles will behave ballistically in the large Stokes limit but act as fluid tracers in the small Stokes limit.

In real expiratory events, exhaled particles will exhibit a distribution of velocities that may deviate from the local air flow due to complex interactions in the upper respiratory tract. Motivated by the fact that the majority of particles lie within the first size bin, where Stokes numbers are small St ∈ [0.0054, 4.87], we treat the particles as fluid tracers at the orifice exit and specify their initial velocity to be the fluid velocity interpolated at the particle position. Further details are provided in Appendix A. We emphasize that this assumption of zero interphase slip velocity at the orifice exit is a significant assumption made within the present work. Future experimental studies will be required to find the extent at which this assumption can be considered valid.

The simulations are solved in Eulerian-Lagrangian framework, where individual particles are treated in a Lagrangian manner, and the gas phase is solved on a background Eulerian mesh. Due to the low concentrations considered in this study, volume fraction effects and two-way coupling between the phases are neglected. The governing equations for the incompressible carrier phase are given by

and

where u u u = [u, v, w] T is the fluid velocity, p is the hydrodynamic pressure, and g g g = [0, −g, 0] T is the gravitational acceleration with g = 9.8 m/s 2 . The equations are implemented in the framework of the NGA code 41 . The Navier-Stokes equations are solved on a staggered grid with second-order spatial accuracy for both convective and viscous terms, and the semi-implicit Crank-Nicolson scheme is used for time advancement maintaining overall second-order accuracy.

Particles are treated in a Lagrangian manner where the translational motion of an individual particle 'i' is given by

where

p ] T are the instantaneous particle position and velocity, respectively, and m p = ρ p πd 3 p /6 is the particle mass.

is the resolved fluid stress at the particle location and f f f (i) drag accounts for unresolved stress due to drag. In this work, the classic Schiller and Naumann drag correlation 42 is used to account for finite Reynolds number effects, given by

where

p ] is the fluid velocity at the location of particle 'i' and

p d p /ν is the particle Reynolds number. The particle equations are advanced in time using a second-order Runge-Kutta scheme.

We briefly note that the present work does not consider thermodynamic effects, such as evaporation and buoyancy, in order to isolate the role of pulsatility on particle dispersion and minimize the parameter space under study. Recent experiments showed that thermal effects are small until the jet speeds are reduced to ambient speeds 38 . Thus, the present work focuses on the near-mouth region, where the unsteadiness of the expiratory events is expected to have a more pressing role on particle dynamics.

A. Pulsatile free-shear jet Figure 5 shows visualizations of the single-, two-, and three-pulse jets at t = 0.75 s, immediately after the final pulse is complete (cf., Fig. 3 ). Inspection of the vorticity magnitude ω ω ω , with ω ω ω = [ω x , ω y , ω z ] T , reveals distinct differences between the three cases. Vortical structures are visualized using the Q-criterion 43 , defined as the second invariant of the velocity gradient tensor, given by

where Ω Ω Ω = (∇u u u − ∇u u u T )/2 and S S S = (∇u u u + ∇u u u T )/2 are the antisymmetric and symmetric components of the velocity gradient, respectively. Physically speaking, the Q-criterion represents a local balance between shear strain and vorticity, with vortices being defined by regions where rigid body rotation is greater than the rate-of-strain.

To demonstrate the effect of pulsatility on the fluid phase, we first consider the location of vortical structures at the end of the third pulse (see Fig. 5 ). For the single-pulse case, a primary vortex ring structure is generated at the downstream edge of the jet while vorticity is minimal near the orifice exit. By contrast, the two-pulse and three-pulse cases exhibit multiple vortex ring structures, corresponding to the number of pulses, with comparatively higher regions of vorticity upstream near the orifice. Vortical structures in the near-orifice region, which are absent from the single-pulse case, will impact the transport of low inertia particles. Specifically, the higher vorticity levels observed in two-and three-pulse cases are expected to accelerate and entrain latestage injected particles to a larger degree when compared to the single-pulse case. However, the strength of the leading vortex structure for two-and three-pulse cases is significantly attenuated from the single-pulse case. The role of pulsatility on particle dynamics is reserved for § III B.

Entrainment of the surrounding air into the jet plays a key role on its transport properties. In the seminal paper by Morton et al. 44 , it was suggested that entrainment, defined as the mean radial velocity at the edge of an axisymmetric boundary layer (in this case edge of the jet), is

proportional to the axial velocity, i.e., u r = α u , where u r = (vy + wz)/r is the radial velocity with r = y 2 + z 2 the radial position. Due to the lack of statistical stationarity in the present configuration, angled brackets used herein denote an average in the circumferential direction but not in time. Pham et al. 45 suggested that the coefficient of entrainment, α, can be estimated as

where δ is the momentum balance length scale defined as

Taub et al. 35 showed that α is approximately constant when the jet has achieved self similarity. For the pulsatile transient jets considered here, the centerline velocity u | r=0 is zero near the orifice when the pulses complete, resulting in an ill-defined α. To this end, we propose an alternative definition based on the jet bulk velocity U 0 according to

s for each case is summarized in Fig. 6 . The corresponding specific momentum M = ∞ 0 u 2 r dr, normalized by the maximum value M 0 = ∞ 0 U 2 0 r dr = 8 × 10 −4 m 4 /s 2 , is also reported to indicate the instantaneous location of each pulse. It can be seen that the single-pulse case generates significantly more momentum downstream at the jet front, while the two-pulse and three-pulse cases exhibit a wider distribution in momentum in the streamwise direction. In addition, the momentum profile is bi-modal for the two-pulse case and tri-modal for the three-pulse case at t = 1 s, where each peak coincides with the peak amplitude of each pulse. As a result, larger values of momentum are observed near the orifice for the two-pulse and three-pulse cases, which proceed to decay as the flow propagates downstream.

The entrainment coefficient is seen to be positive for all three cases at x/D < 20 when t = 1 s, which indicates a net entrainment of ambient air into the jet. Beyond this point (x/D > 20) α becomes negative with the largest magnitude in concert with the first pulse. Note that unlike in the more traditional statistically stationary jet, where α is positive for all x, here the primary vortex ring generated by the first pulse induces a net momentum flux from the jet into the ambient air, resulting in α < 0. By comparing the three cases at t = 1 s, it is observed that α is larger in the upstream regions (x/D < 16) for the multi-pulse cases compared to the single-pulse case due to the vortical structures generated by subsequent pulses. In addition, the three-pulse case exhibits significantly larger entrainment at the jet front where the primary vortex ring resides.

Pulsatility is also observed to affect late-time single-phase dynamics of the jet. By comparing Figs. 6a-6c, the primary vortex of the two-pulse jet is seen to travel at a higher speed, followed by the single-pulse case then the three-pulse case. This results in noticeable differences in the entrainment coefficient between each case at late times. The following sections seek to understand how these combined effects influence particle entrainment and dispersion.

B. Role of pulsatility on particle dispersion of the cloud is slightly hindered with increased pulsatility, although this effect is minor. More pronounced is the penetration of particles associated with subsequent pulses. This is best seen in the three-pulse case, where particles emanating from the second pulse (colored blue) penetrate to the cloud front when t > 1.5 s, despite being injected with lower velocity. In addition, particles from the third pulse (colored red) in the three-pulse case penetrate nearly as far as particles from the second pulse in the two-pulse case. In summary, the geometric centers associated with particles of later pulses travel further downstream with increased pulsatility, and advance upon par- In studying the entrainment characteristics of particle-laden channel flows, Marchioli and Soldati 46 presented a framework highlighting the use of instantaneous joint correlations of nonvanishing components of the fluctuating velocity gradient tensor to determine locations of particle preferential concentration. In addition, the relationship between particle size and their correspond- ing fluid topology was exploited to identify particle preferential sampling. The present work extends these techniques to correlate the range of particle sizes seen in a pulsatile cough with their immediate fluid environment and vortex structures. This is accomplished by correlating the value of the Q-criterion against a particle's directional velocity (u p , v p , w p ) (see Fig. 11 ). The sign of Q crit describes the fluid environment that the particle currently resides, with Q crit > 0 corresponding to regions of high vorticity, Q crit ≈ 0 corresponding to regions of no flow or constant strain, and Q crit < 0 corresponding to regions of high strain rate. The particle directional velocity describes the dispersive characteristics of the particles, with a large spread in velocity being associating with high directional dispersion.

Particles are grouped into three size ranges as depicted in Fig. 4 . It is first observed that small for small particles, indicating that gravity has an effect on the former but not the latter.

It can also be seen that the distribution in lateral velocity, w p , associated with mid-sized particles is narrower in the cases with multiple pulses compared to the single-pulse case. Specifically, mid-sized particles are preferentially sampling near-zero values of w p for a wider range of Q crit compared to the one-pulse case. The distribution in the gravity-aligned (y) direction remains unchanged between the three cases, which is expected as gravity plays a more important role for mid-and large-sized particles in this direction. For the scatter plots in x-direction, although most particles in all three cases are preferentially sampling the right half plane (u p > 0) as the net particle flux is positive in the streamwise direction, more mid-sized particles are seen to have negative u p over a wider range of negative Q crit . These aforementioned differences in x and z are likely a result of mid-sized particle being entrained by the vortices generated by the subsequent pulses as they respond most effectively to turbulent eddies due to their intermediate Stokes numbers.

Consequently, it also explains the decreasing penetration (x c ) and dispersion (z rms ) with increasing pulsatility, as seen in Figs. 8 and 9 . ing to the high-speed release of the fluid-particle mixture. In this phase, penetration is modeled as a function of specific momentum flux M according to

where C M is a constant coefficient of the first regime. The second phase is dominated by "pufflike" dynamics, characterized by the self-similar growth of the puff cloud. The puff penetration evolution is given by

where C I is a constant coefficient of the second regime, and the total specific momentum of the cloud, I, is defined as

Here, α 1 and α 2 are particle entrainment coefficients, which satisfy

/N p is the geometric mean radius of the particle cloud. The particle entrainment coefficients of the two regimes for all three cases are determined by the slope of r c versus x c as shown in Figs. 12a-12c. The entrainment coefficient of the second regime, α 2 , is observed to be larger and more sensitive to pulsatility than the entrainment coefficient of the first regime, α 1 , indicating that the vortex interactions may be more significant in the puff-like regime. Penetration can be predicted by combining Eqs. (11) and (12) using these values of α 1 and α 2 . In contrast to Bourouiba et al. 6 , here M is time-varying due to the pulsatile nature of the expiratory flow. To this end, the average specific momentum M is used in Eq. (11), given by

where M is assumed to follow the pulsatile profile given by Eq. (1) Here, we extend the current model to account for pulsatility. Instead of applying the con-servation law to the entire cloud, the particle concentration coefficient of each pulse is extracted separately for the pulsatile cases. As shown in Figs. 13a-13c, the particle cloud from each pulse follows their own two-stage evolution. In addition, α 1 is observed to be similar between different pulses, whereas α 2 of the second and third pulse are significantly larger than the first pulse, indicating a larger dispersion for late-injected particles as facilitated by the earlier pulses despite the overall dispersion is hindered. Using these values of α 1 and α 2 , the penetration of each pulse is then modeled by following the same procedure described earlier. Let x 1 c , x 2 c , and x 3 c denote the penetration of the first, second and third pulse, and t 1 , t 2 , and t 3 denote the time when the first, second and third pulse complete, the penetration of the entire cloud for the two-pulse case is modeled as the weighted average of each pulse given by

and similarly for the three-pulse case

Pulsatility is now incorporated in the model by explicitly superimposing of the two-stage dynamics of each pulse. Figures 13d-13f show the corrected model predictions. It can been seen that the oscillatory trend of the first jet-like regime is accurately predicted by leveraging the particle entrainment coefficient obtained from of each pulse instead of from the entire particle cloud. Note that this modified model can be readily applied to different coughing profiles or extended to speech patterns (which are essentially a continuous train of expulsions 38, 47 ) to investigate other effects on particle penetration.

In this work, direct numerical simulations of particle-laden turbulent pulsatile jets were conducted to assess the role of pulsatility on particle dynamics. Realistic turbulence was provided at the jet orifice using data obtained from an auxiliary simulation of a turbulent pipe flow. The flow rate of the incoming turbulence was modulated in time according to a damped sine wave that provides control over the number of pulses, their duration, and peak amplitude. Particles were injected in the flow with diameters sampled from a lognormal distribution informed by experimental measurements from the literature.

Vortex structures were analyzed for single-, two-, and three-pulse jets. Qualitative comparison of Q-criterion revealed that the two-pulse and three-pulse cases exhibit multiple vortex ring structures with high vorticity regions persisting for longer times near the orifice. Entrainment coefficients were found to be larger for the multi-pulse cases compared to the single-pulse case due to the vortical structures generated by subsequent pulses, with their largest magnitude in concert with the pulses.

Particle dispersion and penetration were found to be hindered by increased pulsatility. However, particles emanating from later pulses traveled further downstream with increased pulsatility due to acceleration by vortex structures. The observed increase in penetration by later pulses may prove to be significant when determining distances at which an infectious person can be harmful to others, especially later expulsions have been found to contain higher viral concentrations. The evolution of the particle cloud penetration was then compared to an existing single-pulse model by Bourouiba et al. 6 . While the penetration for all three cases are well predicted by the pufflike regime (x c ∼ t 1/4 ) at late time, they deviate from the jet-like regime (x c ∼ t 1/2 ) at early time and instead exhibit oscillations for the pulsatile cases. A modified model was therefore proposed to account for pulsatility by leveraging the particle entrainment coefficient of each pulse and has been shown to accurately predict the oscillatory trend of the early-stage penetration.

The data that support the findings of this study are available from the corresponding author upon reasonable request. We consider a domain of size 10D × D × D, discretized using 326 × 256 × 256 grid points (see The velocity field at steady state is saved and used to prescribe the boundary condition in the main simulation. At each simulation timestep, the velocity in the yz-plane is interpolated from the fine-scale auxiliary simulation onto the coarser domain boundary of the main simulation. The fluid velocity is then rescaled to achieve the desired time-dependent flow rate according to Eq. (1).

The number of particles injected into the main simulation is adjusted each time step to obtain the same mass flow rate as the fluid. The velocity assigned to each particle is equal to the fluid velocity interpolated to its location. This assumes that particles expelled from the orifice have zero interphase slip velocity, and thus zero initial drag.

Proc. Natl. Acad. Sci. USA

Bubbles, drops, and particles

Proceedings of the Summer Program

The computing resources and assistance provided by the staff of the Advanced Research Computing at the University of Michigan, Ann Arbor are greatly appreciated. We would also like to acknowledge the National Science Foundation for partial support from awards CBET 2035488 and 2035489.