Role of pulsatility on particle dispersion in expiratory flows

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¹ while pyroclastic density currents, generated by large density differences between gas-particle mixtures, feeds explosive volcanic eruptions.² 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¹⁹ (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.²⁰ Gupta et al.²³ 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²⁴ during speech, coughing, or sneezing—bronchiolar (droplet size $< 1 − 5 μ$m); laryngeal (⁠ $5 − 20 μ$m); and oral (⁠ $> 50 μ$m).²⁵ If the site of severe infection is deeper in the lungs (bronchiolar), the expelled aerosols/droplets are generated by a “fluid-film burst” mechanism²⁶ 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.³⁵ 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 %$⁠.³⁷ 

Turbulent transport in statistically stationary, axisymmetric, free jets has been well characterized experimentally^39–41 and numerically.^42–44 Chein and Chung¹⁹ 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⁴⁵ 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.⁴⁶ 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³. The fluid is considered to be air with a density of $ρ = 1.172$ kg/m³ and kinematic viscosity of $ν = 1.62 × 10 − 5$ m²/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 × 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.⁴⁶ 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 $∼ 7 %$ of the peak inflow velocity U₀ and was observed to have a 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^13,47) corresponding to a bulk Reynolds number $Re b = U 0 D / ν = 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₀ 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.²³ we propose a self-similar profile for the bulk velocity resulting from multiple expulsions.

To accurately characterize the particle size distribution generated by coughing, we employ a lognormal distribution fit to the experimental measurements of Duguid.⁴⁸ The particle diameter ranges between $1 ≤ d p ≤ 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.⁴⁹ 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 η = τ p / τ η$⁠, may be utilized to gauge the role of particle inertia, where $τ p = ρ p d p 2 / ( 18 ρ ν )$ is the particle response time, $τ η = ( ν / ε ) 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 ∈ [ 1 , 30 ] , [ 30 , 60 ]$⁠, and $[ 60 , 100 ] μ$m, which yields the following Stokes number ranges St $∈ [ 0.0054 , 4.87 ] , [ 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 $∈ [ 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.

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.⁴⁷ 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).⁵² 

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.

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 = ∫ 0 ∞ 〈 u 〉 2 r d r$⁠, normalized by the maximum value $M 0 = ∫ 0 ∞ U 0 2 r d r = 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 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 $α < 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.

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 = ∑ 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 = ∑ 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_c and $z rms$ are seen to oscillate in the two- and three-pulse cases. Additionally, the final displacement of x_c and $z rms$ at 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⁴⁷ that the unsteadiness of expiratory events diminishes far from the mouth.

In studying the entrainment characteristics of particle-laden channel flows, Marchioli and Soldati⁵⁶ 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 ≈ 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 ∈ [ 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 ∈ [ 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 ∈ [ 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_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.

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 ∼ t 1 / 4$⁠. For the first jet-like regime, however, the penetration profiles deviate from $x c ∼ t 1 / 2$ and instead exhibit oscillations for the pulsatile cases, which is not correctly captured by the model.

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.

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.¹³ While the penetration for all three cases is well predicted by the puff-like 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 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.

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 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 × D × D$⁠, discretized using $326 × 256 × 256$ grid points (see Fig. 14). The grid spacing is chosen such that $Δ y + = Δ z + = 1.25$ and $Δ x + = 9.8$ to satisfy the resolution criteria of DNS for pipe flows,^60,61 where $( · ) + = ( · ) u τ / ν$ denotes frictional wall units with $u τ$ 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 D / U 0$⁠. A comparison of the velocity statistics against DNS experimental data from the literature⁶² 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.