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. Thermodynamic effects, such as evaporation and buoyancy, are neglected in order to isolate the role of pulsatility on particle dispersion. 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 are found to be hindered by increased pulsatility. However, the penetration of particles emanating from a secondary or tertiary expulsion is enhanced due to acceleration downstream by vortex structures.

## I. INTRODUCTION

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,4}

Recent studies considering the airborne transmission of COVID-19 have largely determined
that the prescribed social distancing rules may not be sufficient to protect against
host-to-host transmission.^{5–8} Additionally, when incorporating environmental factors, such as temperature, humidity, and
wind speed, the traveled distance and dispersion of aerosols are seen to travel far beyond
the typical 6 ft social distance guideline.^{9,10} 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,^{11,12} experimental,^{13,14} and computational^{14–18} 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^{19} (see Fig. 1).

Experimental measurements have demonstrated that realistic coughs are pulsatile, involving
a sequence of coughing events, sometimes referred to as “cough epochs.”^{21,22} The flow rate associated 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.^{20} Gupta *et al.*^{23} 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 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^{24} during speech, coughing, or sneezing—bronchiolar (droplet size $ < 1 \u2212 5 \u2009 \mu $m); laryngeal ( $ 5 \u2212 20 \u2009 \mu $m); and oral ( $ > 50 \u2009 \mu $m).^{25} If the
site of severe infection is deeper in the lungs (bronchiolar), the expelled
aerosols/droplets are generated by a “fluid-film burst” mechanism^{26} during the 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,^{24,27} 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 multipulse cough will expel the volume of air originating predominantly from
deeper within the lungs and are 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 vortex–vortex 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.^{28–32} When the Stokes number, defined as the
ratio of particle-to-fluid 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.^{33,34} Gualtieri *et
al.*^{35} 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.^{36–38} The gravitational settling of aerosol particles
can be enhanced by this mechanism by as much as $ 50 % $.^{37}

Turbulent transport in statistically stationary, axisymmetric, free jets has been well
characterized experimentally^{39–41} and numerically.^{42–44} Chein and
Chung^{19} 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^{45} 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.*^{46} 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.

## II. SIMULATION DETAILS

### A. Flow configuration

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 $ \rho p = 998 $ kg/m^{3}. The fluid is considered to be air with a
density of $ \rho = 1.172 $ kg/m^{3} and kinematic viscosity of $ \nu = 1.62 \xd7 10 \u2212 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 = 40 D $, $ L y = 20 D $, and $ L z = 20 D $, 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

*y*- and

*z*-directions. The spanwise grid spacing varies between $ 4.98 \xd7 10 \u2212 4 $ m $ \u2264 \Delta y , \Delta z \u2264 2.1 \xd7 10 \u2212 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.

^{46}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 coflow is introduced along the positive

*x*-direction with a velocity magnitude 0.32 m/s. The coflow is $ \u223c 7 % $ of the peak inflow velocity

*U*

_{0}and was observed to have a negligible effect on the particle dynamics.

### B. Pulsatile inflow

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^{13,47}) corresponding to a
bulk Reynolds number $ Re b = U 0 D / \nu = 4938 $. Further details on the pipe flow simulation are provided
in Appendix. 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 ( 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.*^{23} we propose a self-similar profile for the
bulk velocity resulting from multiple expulsions.

*Q*(

*t*) is given by a damped sine wave according to

^{23}The total duration of each profile varies but the inputs are defined to yield the same volume of expelled air so a fair comparison can be drawn between cases with different numbers of pulses.

### C. Particle injection

To accurately characterize the particle size distribution generated by coughing, we
employ a lognormal distribution fit to the experimental measurements of Duguid.^{48} The particle diameter ranges between $ 1 \u2264 d p \u2264 100 $ *μ*m with a mean of 24 *μ*m and a standard deviation of
17.9 *μ*m as shown in Fig. 4. At
each simulation time step, particles are introduced at the inflow plane by assigning them
a random position within the orifice and a diameter that is sampled from the
aforementioned lognormal distribution. The number of particles per time step is adjusted
dynamically to achieve the same mass flow rate used for the fluid. Given the prescribed
expulsion volume and particle size distribution, approximately 15 000 particles are
generated at the end of a coughing spell, representative of the typical quantity observed
in experiments.^{49} For a three-pulse
case, this corresponds to roughly 8200, 4600, and 2200 particles being injected during the
first, second, and third pulses, respectively.

The turbulence Stokes number, $ St \eta = \tau p / \tau \eta $, may be utilized to gauge the role of particle inertia,
where $ \tau p = \rho p d p 2 / ( 18 \rho \nu ) $ is the particle response time, $ \tau \eta = ( \nu / \epsilon ) 1 / 2 $ is the Kolmogorov timescale 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 \u2208 [ 1 , 30 ] , \u2009 [ 30 , 60 ] $, and $ [ 60 , 100 ] \u2009 \mu $m, which yields the following Stokes number ranges St $ \u2208 [ 0.0054 , 4.87 ] , \u2009 [ 4.87 , 19.50 ] $, and $ [ 19.50 , 54.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 airflow 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 $ \u2208 [ 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. 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.

### D. Governing equations

*p*is the hydrodynamic pressure, and $ g = [ 0 , \u2212 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.

^{50}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.

*i*” is given by

^{51}is used to account for finite Reynolds number effects, given by

*i”*and $ Re p = \Vert u [ x p ( i ) ] \u2212 v p ( i ) \Vert d p / \nu $ 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. Thus, this study does not
consider particle–particle interactions, such as coalescence and the size of individual
particles are held constant throughout the duration of the simulated cough. Recent
experiments showed that thermal effects are small until the jet speeds are reduced to
ambient speeds.^{47} 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. In addition, it was recently
demonstrated that preferential concentration can increase local humidity that in turn
prevents evaporation and extends the lifetime of droplets significantly (by as much as two
orders of magnitude).^{52}

## III. RESULTS AND DISCUSSION

### A. Pulsatile free-shear jet

#### 1. Flow visualization

*t*=

*0.75 s, immediately after the final pulse is complete (cf. Fig. 3). Inspection of the vorticity magnitude $ \Vert \omega \Vert $, with $ \omega = [ \omega x , \u2009 \omega y , \omega z ] T $, reveals distinct differences between the three cases. Vortical structures are visualized using the*

*Q*-criterion,

^{53}defined as the second invariant of the velocity gradient tensor, given by

*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 late-stage 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 Sec. III B.

#### 2. Entrainment

*et al.*,

^{54}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., $ \u3008 u r \u3009 = \alpha \u3008 u \u3009 $, where $ u r = ( v y + w z ) / 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.*

^{55}suggested that the coefficient of entrainment,

*α*, can be estimated as

*δ*is the momentum balance length scale defined as

*et al.*

^{44}showed that

*α*is approximately constant when the jet has achieved self-similarity. For the pulsatile transient jets considered here, the centerline velocity $ \u3008 u \u3009 | 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

The entrainment coefficient *α* as a function of streamwise location *x*/*D* at *t *=* *1 s,
1.5 s, and 2 s for each case is summarized in Fig. 6. The corresponding specific momentum $ M = \u222b 0 \u221e \u3008 u \u3009 2 r \u2009 d r $, normalized by the maximum value $ M 0 = \u222b 0 \u221e U 0 2 r \u2009 d r = 8 \xd7 10 \u2212 4 \u2009 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 trimodal 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 $ \alpha < 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 multipulse 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. 6(a)–6(c), 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. Sections III B–III C seek to understand how these combined effects influence particle entrainment and dispersion.

### B. Role of pulsatility on particle dispersion

Visualizations of the particle cloud at three instances in time are shown in Fig. 7. Particles are colored by the corresponding pulse
they were injected into. At each time step, the total number of particles associated with
each pule, $ N p ( t ) $, is identified and used to define the geometric center of
the cloud according to $ x c = \u2211 i = 1 N p x p ( i ) / N p $. Qualitative differences in the cloud evolution can be seen
between the three cases. Careful inspection of Fig. 7 (Multimedia view) reveals that the overall penetration 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 particles
injected earlier. The observed increase in penetration by secondary and tertiary pulses
has important connotations as it is expected that later expulsions could contain higher
viral concentrations depending on the location in the respiratory tract where the
infection resides. Therefore, the aforementioned phenomena may prove to be significant
when determining distances at which an infectious person can pose a risk to others; as the
enhanced transport of high viral load droplets is expected to increase the probability of
transmission.

Particle dispersion is characterized herein by the root mean square (rms) of lateral
particle position, given by $ z rms = \u2211 i = 1 N p z p ( i ) 2 / N p $. The temporal evolution of the cloud penetration, *x _{c}*, and dispersion, $ z rms $, associated with each pulse are shown in Figs. 8 and 9.
During the early-stage injection,

*x*and $ z rms $ are seen to oscillate in the two- and three-pulse cases. Additionally, the final displacement of

_{c}*x*and $ z rms $ at

_{c}*t*=

*2 s is lower than that of the single-pulse case. These observations can be attributed to the decreasing volumetric flow rate between each pulse, and as a consequence, the average injection velocity is lower compared to the single-pulse case. After the pulsatile injection completes, however, the streamwise displacement of the cloud associated with later pulses is seen to “catch up” with the earlier pulses despite being injected at later time and with significantly lower velocity. For example, the geometric center of the second-pulse particles nearly coincides with the overall particle cloud for the three-pulse case at*

*t*=

*2 s, indicating once again that the penetration of potentially more contagious particles from later pulses is accelerated by earlier pulses. On the contrary, such an effect of pulsatility is not observed for the lateral dispersion for which particles from earlier pulses disperse further.*

The velocity of the geometry center of each pulse associated with the three-pulse case, $ u c = d x c / d t $, is shown in Fig. 10.
Due to the finite inertia of the particles, *u _{c}* lags the fluid
velocity during early-stage injection. It can be seen that the peak velocity of the third
pulse is only reduced by a factor of two compared to the first pulse, despite its
injection velocity being reduced by a factor of four. In addition, the velocity associated
with the second pulse exceeds the velocity of the first pulse when

*t*>

*0.3 s, further demonstrating the ability of particles injected at later times to catch up to particles near the front of the cloud. Further downstream of the inlet, the velocities of each geometric center converge, consistent with observations from human subjects*

^{47}that the unsteadiness of expiratory events diminishes far from the mouth.

In studying the entrainment characteristics of particle-laden channel flows, Marchioli
and Soldati^{56} 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 corresponding
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 \u2248 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 particles ( $ d p \u2208 [ 1 , 30 ] $ *μ*m) equally sample all values of $ Q crit $, exhibiting no preferential location in terms of fluid
vortical structures. Mid-sized particles ( $ d p \u2208 [ 30 , 60 ] $ *μ*m) tend to sample regions of high strain rate, indicating their ejection
from vorticity-dominated regions, i.e., classical preferential concentration. Large
particles ( $ d p \u2208 [ 60 , 100 ] $ *μ*m) are seen to sample regions of constant strain ( $ Q crit = 0 $) as a consequence of them falling out of the cloud due to
gravity. Note that the distribution of mid-sized particles are skewed to negative values
of *v _{p}*, while this is not observed 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*for a wider range of $ Q crit $ compared to the one-pulse case. The distribution in the gravity-aligned (

_{p}*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*over a wider range of negative $ Q crit $. These aforementioned differences in

_{p}*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*) and dispersion ( $ z rms $) with increasing pulsatility, as seen in Figs. 8 and 9.

_{c}### C. Theoretical modeling of respiratory emissions

*et al.*

^{13}proposed a theoretical model for the cloud penetration based on the notion of conservation of cloud momentum. It has been generally observed from their experiments that two phases of the cloud evolution exist. The first phase is dominated by jet-like dynamics, corresponding 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

*C*is a constant coefficient of the first regime. The second phase is dominated by “puff-like” dynamics, characterized by the self-similar growth of the puff cloud. The puff penetration evolution is given by

_{M}*C*is a constant coefficient of the second regime, and the total specific momentum of the cloud,

_{I}*I*, is defined as

*r*vs

_{c}*x*as shown in Figs. 12(a)–12(c). The values are consistent with the range previously reported by experiments of Bourouiba

_{c}*et al.*

^{13}( $ 0.13 \u2212 0.24 $) and Gupta

*et al.*

^{23}( $ \u223c 0.21 $). The entrainment coefficient of the second regime, $ \alpha 2 \u2032 $, is observed to be larger and more sensitive to pulsatility than the entrainment coefficient of the first regime, $ \alpha 1 \u2032 $, 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 $ \alpha 1 \u2032 $ and $ \alpha 2 \u2032 $. In contrast to Bourouiba

*et al.*,

^{13}here

*M*is time-varying due to the pulsatile nature of the expiratory flow. To this end, the average specific momentum $ M \xaf $ is used in Eq. (11), given by

*M*is assumed to follow the pulsatile profile given by Eq. (1), i.e., $ M ( t ) = M 0 | e \u2212 t / \tau \u2009 sin \u2009 ( \omega t ) | $. The coefficients

*C*and

_{I}*C*are determined by least-square fitting of the puff regime and solving $ ( 4 C M 2 M / \alpha 1 \u2032 2 ) 1 / 4 t cr 1 / 2 = ( 4 C I I / \alpha 2 \u2032 3 ) t cr 1 / 4 $, respectively, with $ t cr $ the intersection time of the two regimes.

_{M}The model predictions are displayed in Figs. 12(d)–12(f). It can be seen that the second puff-like regime for all three cases scales as $ x c \u223c t 1 / 4 $. For the first jet-like regime, however, the penetration profiles deviate from $ x c \u223c t 1 / 2 $ and instead exhibit oscillations for the pulsatile cases, which is not correctly captured by the model.

*t*

_{1},

*t*

_{2}, and

*t*

_{3}denote the time when the first, second, and third pulse completes, the penetration of the entire cloud for the two-pulse case is modeled as the weighted average of each pulse given by

Pulsatility is now incorporated in the model by explicitly superimposing of the two-stage
dynamics of each pulse. Figures 13(d)–13(f) 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^{47,57}) to investigate other effects on particle
penetration.

## IV. CONCLUSIONS

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 multipulse
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 physical distances at which an infectious person can pose a risk to others, especially since later expulsions have been found to contain higher viral concentrations. Specifically, as it was previously observed that particles from subsequent pulses were found to penetrate the cloud front, this work recommends larger-timescale studies of realistic/pulsatile coughs to verify particle dispersion against single-pulse events. Additionally, measures to dampen the strength of vortex structures in a pulsatile cough, such as home-grade cloth masks, can prove to be beneficial in reducing the penetration distances of subsequent pulses with potentially higher viral concentrations. As such, future work should focus on investigating the efficacy of common flow barriers when exposed to pulsatile particle-laden flows.

The evolution of the particle cloud penetration was then compared to an existing
single-pulse model by Bourouiba *et al.*^{13} While the penetration for all three cases is well predicted by the
puff-like regime ( $ x c \u223c t 1 / 4 $) at late time, they deviate from the jet-like regime ( $ x c \u223c 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.

## ACKNOWLEDGMENTS

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 Award Nos. CBET 2035488 and 2035489.

### APPENDIX: TURBULENT INFLOW GENERATION

To accurately model an inlet condition resembling the expiratory turbulent flow exiting
from a human mouth, a direct numerical simulation (DNS) of single-phase flow traversing
through a cylindrical pipe is performed. The pipe diameter $ D = 0.02 \u2009 m $ is representative of typical mouth opening. The fluid-phase
equations are discretized on a Cartesian mesh, and a conservative immersed boundary (IB)
method is employed to model the cylindrical pipe geometry without requiring a body-fitted
mesh. The method is based on a cut-cell formulation that requires rescaling of the
convective and viscous fluxes in these cells and provides discrete conservation of mass
and momentum.^{58,59}

We consider a domain of size $ 10 D \xd7 D \xd7 D $, discretized using $ 326 \xd7 256 \xd7 256 $ grid points (see Fig. 14). The grid spacing is chosen such that $ \Delta y + = \Delta z + = 1.25 $ and $ \Delta x + = 9.8 $ to satisfy the resolution criteria of DNS for pipe
flows,^{60,61} where $ ( \xb7 ) + = ( \xb7 ) u \tau / \nu $ denotes frictional wall units with $ u \tau $ the friction velocity. Periodic boundary conditions are
applied in the streamwise direction. A uniform source term resembling a mean pressure
gradient is added to the right-hand side of Eq. (3) and adjusted dynamically to maintain the desired flow rate. The flow is
initialized with a bulk velocity $ U 0 = 4.0 $ m/s with $ 10 % $ sinusoidal fluctuations to accelerate the transient
process. A statistical stationary state is reached after $ 240 \u2009 D / U 0 $. A comparison of the velocity statistics against DNS
experimental data from the literature^{62} is provided in Fig. 15.

The velocity field at steady state is saved and used to prescribe the boundary condition
in the main simulation. At each simulation time step, 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.

## REFERENCES

*Bubbles, Drops, and Particles*

*Proceedings of the Summer Program*