Investigating the parametric dependence of the impact of two-way coupling on inertial particle settling in turbulence

Abstract Tom et al. (J. Fluid Mech., vol. 947, 2022, p. A7) investigated the impact of two-way coupling (2WC) on particle settling velocities in turbulence. For the limited parameter choices explored, it was found that (i) 2WC substantially enhances particle settling compared with the one-way coupled case, even at low mass loading $\varPhi _m$ and (ii) preferential sweeping remains the mechanism responsible for the particles settling faster than the Stokes settling velocity in 2WC flows. However, significant alterations to the flow structure that can occur at higher mass loadings mean that the conclusions from Tom et al. (J. Fluid Mech., vol. 947, 2022, p. A7) may not generalise. Indeed, even under very low mass loadings, the influence of 2WC on particle settling might persist, challenging the conventional assumption. We therefore explore a much broader portion of the parameter space, with simulations covering cases where the impact of 2WC on the global fluid statistics ranges from negligible to strong. We find that, even for $\varPhi _m=7.5\times 10^{-3}$, 2WC can noticeably increase the settling for some choices of the Stokes and Froude numbers. When $\varPhi _m$ is large enough for the global fluid statistics to be strongly affected, we show that preferential sweeping continues to be the mechanism that enhances particle settling rates. Finally, we compare our results with previous numerical and experimental studies. While in some cases there is reasonable agreement, discrepancies exist even between different numerical studies and between different experiments. Future studies must seek to understand this before the discrepancies between numerical and experimental results can be adequately addressed.


Introduction
Turbulent, dispersed, multiphase flows are not yet fully understood due to their complex, multiscale dynamics.The dispersed phase may consist of solid inertial particles, bubbles or gases and their interaction with the carrier fluid depends on a number of factors including their density (relative to the fluid density), their total mass (relative to the total fluid mass) and the nature of the turbulence.Despite their complexity, experimental, computational and theoretical investigations have made significant progress in understand various aspects of such flows, including the dispersion and clustering of particles, droplet breakdown and coalescence, and turbulent modulation (see Balachandar & Eaton 2010;Brandt & Coletti 2022;Gustavsson & Mehlig 2016, for reviews).This has led to advances in our understanding of diverse problems including cloud microphysics (Pruppacher & Klett 1997;Shaw 2003a), aerosol deposition in human lungs (Ou et al. 2020), and planetesimal formation (Cuzzi et al. 2000;Birnstiel et al. 2016).
We focus on the problem of the settling velocity of inertial particles in isotropic turbulence for particles that are small compared with the Kolmogorov length scale, and are much denser than the fluid.Many experimental and numerical studies have observed that such inertial particles settle faster in a turbulent flow than they would in a quiescent fluid (Wang & Maxey 1993;Aliseda et al. 2002;Yang & Shy 2005;Good et al. 2014;Bec et al. 2014;Rosa et al. 2016;Monchaux & Dejoan 2017;Momenifar et al. 2019;Momenifar & Bragg 2020;Dhariwal & Bragg 2019;Petersen et al. 2019;Berk & Coletti 2021).Prior to these observations, it had also been shown by Maxey & Corrsin (1986) that in random flow fields inertial particles settle faster than they would in a quiescent flow.Maxey (1987) proposed the preferential sweeping mechanism to explain how turbulence (or fluctuations in a random flow) could cause inertial particles to settle faster than they would in a quiescent flow.In particular, Maxey argued that for weakly-inertial settling particles, the combined effects of the structure of the strain-rate and vorticity fields of the flow, particle inertia and gravity cause the particles to be preferentially swept around the downward moving side of vortices in the flow.Hence, the average fluid velocity along their trajectory points down, leading to an enhancement of their settling velocity.Recently, Tom & Bragg (2019) extended the analysis of Maxey (1987) to explain how the preferential sweeping mechanism works when the particle inertia is finite, and the role that different turbulent flow scales play in the mechanism.This was called the multi-scale preferential sweeping mechanism and it successfully explained a number of previous results that could not be accounted for by the original analysis of Maxey (1987).
Many of the studies mentioned above focused on the one-way coupled (1WC) regime where the effect of the particles on the flow is ignored, in contrast to the two-way coupled (2WC) regime.Other studies have however considered the effect of 2WC on particle settling, including Bosse et al. (2006); Dejoan (2011); Monchaux & Dejoan (2017); Rosa et al. (2021).These studies all found that 2WC causes the particles to settle faster than the 1WC case.The study of Monchaux & Dejoan (2017) also emphasized that 2WC can substantially increase the particle settling velocities compared with the 1WC case even in regimes where the particle mass loading is small enough for the effect of the particles on the global fluid statistics to be negligible.They also argued that in the presence of 2WC, the preferential sweeping mechanism is not responsible for the enhanced settling speeds due to turbulence.Instead, they argued that the dominant contribution comes from a fluid dragging effect, where the particles drag the fluid surrounding them with them as the fall, lowering the drag force on the particles and thereby enhancing their settling velocities.Tom et al. (2022) investigated this in greater detail, applying the multi-scale preferential sweeping mechanism developed in Tom & Bragg (2019) to the 2WC case.In agreement with Monchaux & Dejoan (2017), Tom et al. (2022) found that 2WC can substantially increase the particle settling velocities compared with the 1WC case even in regimes where the particle mass loading is small enough for the effect of the particles on the global fluid statistics to be negligible.However, in contrast to Monchaux & Dejoan (2017), they demonstrated that preferential sweeping continues to be the mechanism responsible for the enhanced settling speeds due to turbulence.In particular they showed that the difference between the 1WC and 2WC cases is that in the 2WC case, the particles are not merely swept around the downward moving side of vortices in the flow but they also drag the fluid down with them in these regions as they fall.Hence, the fluid dragging effect does not replace the preferential sweeping mechanism, but acts with it to further enhance the particle settling velocities compared to the 1WC case.
The study of Tom et al. (2022) was, however, for a limited portion of the parameter space of the problem.In particular, while it considered a range of Stokes numbers ≡ / ∈ [0.3, 2] (where is the particle response time and is the Kolmogorov timescale), it restricted attention to a single Froude number ≡ /( ) = 1 (where is the Kolmogorov velocity scale and is the gravitational acceleration) and volume fraction Φ = 1.5 × 10 −5 .Due to this, crucial open questions remain.In particular, 1) how small must Φ be for the effects of 2WC on the particle settling to be negligible?, 2) does the preferential sweeping mechanism continue to be relevant in 2WC flows at higher volume fractions?This study aims to answer these crucial questions.
The rest of the paper is organised as follows: In Sec. 2 we outline the problem of interest, in Sec. 3 we summarise the numerical methods, in Sec.4we present and discuss the observations of our simulations.Finally, in Sec. 5 we summarise our findings.

Problem Formulation
We consider a dilute suspension (i.e. the volume-fraction Φ is small enough to ignore particle collisions) of small ( / ≪ 1, where is the particle diameter and is the Kolmogorov length scale), dense ( / ≫ 1, where is the particle density and is the fluid density), spherical, inertial particles settling under the force of gravity.Such assumptions are often considered to be suitable for modeling droplet transport in atmospheric clouds (Shaw 2003b;Grabowski & Wang 2013), and dust transport in the atmosphere (Richter & Chamecki 2018).
Assuming that the particle Reynolds number is small, the evolution equation for particles in this regime is given by (Maxey & Riley 1983) where ( ), ( ) denote the particle position and velocity, respectively, ( ( ), ) is the undisturbed fluid velocity at the particle position, is the Kolmogorov time scale, and denotes the gravitational acceleration.
The particles settle in a statistically stationary, homogeneous turbulent flow, which is governed by the incompressible Navier-Stokes equation (written here in rotational form) Here, ( , ) and ≡ ∇ × are the fluid velocity and vorticity fields respectively, ( , ) is the pressure, is the dynamic viscosity, ( , ) is the total momentum feedback of the particles on the carrier fluid, and ( , ) is the large-scale forcing required to maintain steady-state turbulence.In our DNS we will consider the case where isotropically forces the flow, such that when = 0 the flow is statistically isotropic, but it is anisotropic for ≠ 0 due to the momentum coupling term .
We consider the particle motion in both the one-way coupled (1WC) scenario where we set = 0, and the two-way coupled (2WC) scenario where ≠ 0. The particle feedback is the sum of the hydrodynamic forces generated by all particles at location ( ) = where is the mass of each particle, and are the position and velocity of the -th particle, and is the reaction force generated by the -th particle.As discussed in Tom et al. (2022), when ≠ 0 settling particles will generate a mean flow in the direction of .In order to generate a zero mean-flow velocity in the vertical direction the pressure is accordingly modified, a method that was also used in Maxey & Patel (2001) The ensemble-averaged mean particle settling speed can be obtained from Eq. 2.1, assuming statistical homogeneity and stationarity Here, is the Stokes settling speed (the settling speed the particle would have in a quiescent fluid) and ( ( ), ) is the ensemble-average vertical fluid velocity at the particle location.It is this latter contribution that represents the contribution of turbulence to the settling speed.Even when the Eulerian average of the vertical fluid velocity is zero ( , ) = 0, the average along the particle trajectory ( ( ), ) need not be zero because inertial particles will not in general sample the flow field ergodically.As a result of this, through the contribution ( ( ), ) turbulence can lead to either enhanced settling ( ( ) < ) or loitering ( ( ) > ).Results from laboratory experiments, field measurements, and numerical simulations show that for the case of small, dense, inertial particles it is enhanced settling that has been predominantly observed (Wang & Maxey 1993;Yang & Shy 2005;Good et al. 2014;Petersen et al. 2019;Nemes et al. 2017;Li et al. 2021a,b;Bec et al. 2014;Rosa et al. 2016;Ireland et al. 2016;Tom et al. 2022).

Numerical method
We adopt an Eulerian-Lagrangian approach to numerically integrate Eq. 2.1 in conjunction with Eq. 2.2.We perform DNS of statistically stationary, homogeneous turbulence in a three-dimensional, periodic domain using a pseudo-spectral scheme on 3 collocation points (Ireland et al. 2013).A second-order Runge-Kutta time-stepping scheme with an exponential integration is employed for the viscous stress, and a combination of spherical truncation and phase-shifting for the dealiasing schemes.A seventh-order, B-spline polynomial interpolation scheme (that uses eight points) is used to calculate both the fluid velocities at the particle position and to project the particle momentum feedback onto the fluid at the grid points.Equation 2.1 is solved using an exponential integrator method (Hochbruck & Ostermann 2010), which ensures stability and accuracy for low Stokes numbers particles while allowing the same time step to be used for both the fluid and particle equations of motion.We refer the reader to ≡ ′ / = 2 / 5/3 is the Taylor microscale Reynolds Number, is the Taylor microscale, L is the domain length, is the number of grid points in each direction, is the fluid kinematic viscosity, is the mean turbulent kinetic energy dissipation rate, is the integral length scale, ≡ ( 3 / ) 1/4 is the Kolmogorov length scale, ′ ≡ 2K/3 is the r.m.s of fluctuating fluid velocity, K is the turbulent kinetic energy, is the Kolmogorov velocity scale, ≡ / ′ is the large-eddy turnover time, is the Kolmogorov time scale, = 2 /3 is the maximum resolved wavenumber, is the time step.The small-scale resolution, and the total flow kinetic energy measured by ′ are approximately constant between the different simulations.These flow statistics are constructed by averaging over the spatial domain and averaging over a time period of 10 .

Simulation parameters
Table 1 summarises the parameters of our DNS, where L is the domain-length, is the number of collocation points (the grid-size is thus 3 ), and is the viscosity.These same values are used for all of the 1WC and 2WC DNS cases.The table also shows turbulence statistics from the DNS of the unladen flow: the mean dissipation rate , the root mean square (r.m.s.) fluctuating velocity ′ ≡ 2K/3 (where K is the turbulent kinetic energy), the Kolmogorov velocity and length scales, and , respectively, the integral length scale , and the Taylor microscale Reynolds number ≡ ′ / , where is the Taylor lenghtscale.For the 1WC and 2WC simulations, the Stokes number ≡ / and the Froude number ≡ /( ) are defined with respect to the unladen fluid statistics, and we consider Stokes numbers = 0.3, 0.5, 0.7, 1, 2 (referred to as St 1 , St 2 , St 3 , St 4 , and St 5 , respectively), each of which is considered for three different Froude numbers = 0.3, 1, and 3.The particle mass-loading Φ is related to the volume-fraction Φ through Φ = ( / )Φ = ( / ) 3 /6, and we consider particles with constant density / = 5000, the same as that used in Bosse et al. (2006);Monchaux & Dejoan (2017);Tom et al. (2022).Since / is fixed, varying Φ corresponds to varying Φ, and therefore we will often refer to the effect of varying Φ even though strictly speaking the momentum coupling depends on Φ and not simply Φ.For each choice of and , we consider Φ = 1.5 × 10 −6 , 7 × 10 −6 , 1.5 × 10 −5 , 7 × 10 −5 , and 1.5 × 10 −4 which we refer to by Φ 1 , Φ 2 , Φ 3 , Φ 4 , Φ 5 , respectively.However, due to computational resources, we were unable to complete the run with = 2, = 0.3, Φ = 1.5 × 10 −4 .This makes for a total of 59 different DNS, and due to this we are restricted to cases where the unladen flow has ≃ 88.The present study significantly expands on the study of Tom et al. (2022) where only = 1 and Φ = 1.5 × 10 −5 were considered, for the same and values.
3.3.Simulation approach DNS of the unladen flow were run for approximately 20 until a statistically stationary state was reached.For the 1WC runs, the DNS were initialized using the fluid data obtained at the final time step of the unladen DNS, and using random initial particle positions and initial particle velocities equal to the local fluid velocity plus the Stokes settling velocity.For the 2WC simulations, the DNS were initialized using the fluid and particle data obtained at the final time step of the corresponding 1WC DNS (i.e.having the same , , Φ).For (a) both the 1WC and 2WC simulations, we collected statistics over a period of 60 to ensure reasonable statistical convergence (this is a smaller window than that used in Tom et al. (2022) but tests showed that it was sufficient).A table presenting some of the key fluid statistics and parameters from all of the 2WC runs is provided in the Appendix for reference.

Results
4.1.Discussion on fluid statistics Before we explore the parametric dependence of the impact of 2WC on particle settling in turbulence, it is informative to first consider how 2WC affects the global flow statistics.We first consider the probability density function (PDF) P (Q) ≡ ( ( , ) − Q , where denotes a fixed point in the flow, is the second-invariant of the velocity gradient tensor and Q is the sample-space coordinate.In Fig. 1 we show the results for P (Q) for the 2WC simulations for varying , Φ at = 0.3, and the corresponding 1WC result (which is the same as for the unladen flow) is also shown for reference.Here, Q has been normalised using the standard-deviation of (referred to as Q ) so that the PDFs are in standard form, enabling us to compare the shape of the PDFs for the different cases.P (Q) for the 1WC simulations, which is the same as for the unladen flow, are shown in black solid lines.The mean value of Q for 1WC and for all the 2WC simulations is zero, and for the 1WC flows, P (Q) is negatively skewed due to the vorticity field being more intermittent than the strain-rate field.Fig. 1 (a) shows that P (Q) is only weakly affected by 2WC when Φ Φ 3 , which was also observed in Tom et al. (2022) for Φ 3 .For larger Φ, P (Q) becomes more symmetric, with the probability of the tail for Q < 0 reducing, while it increases for Q > 0. This means that 2WC is suppressing the probability of regions of strong enstrophy and enhancing the probability of regions of intense straining motions.Figure 1 (b) shows the dependence of P (Q) for volume fraction Φ 4 .We see that changing the particle inertia has in general little effect in altering the skewness of P (Q), but increases the probability of high-strain and high-rotation regions alike.
Figure 2 shows the turbulent kinetic energy (TKE) dissipation rate for the different 2WC simulations, normalised by the value of from the 1WC runs.We observe that does not change appreciably for small Φ, and shows negligible dependence on the Stokes number .(a) However, for volume fractions Φ > Φ 3 , generally increases with increasing Φ (especially in simulations where the Froude number is = 0.3) and reduces slightly with increasing .It is important to note, however, that these results depend on the type of forcing used in the DNS.In the present simulations, the flow is forced by keeping the flow TKE constant by re-injecting the TKE that is dissipated back into the flow at the lowest wavenumbers, see Ireland et al. (2013).In Carbone et al. (2019) where it was the energy injection rate that was controlled, was found to decrease with increasing Φ in 2WC DNS with = ∞.We now consider the energy spectrum in order to consider the impact of 2WC over the entire range of scales in the flow.Figure 3 shows (a) ( ) the turbulent kinetic energy spectrum and (b) ( )/ 33 ( ) the ratio of the total to the vertical component of the turbulent kinetic energy spectrum versus the wave-number , for 2WC simulations with = 1, = 0.3, and for all volume fractions Φ.The bottom panel in Fig. 3 shown corresponding plots of ( ) (c) and ( )/ 33 ( ) for 2WC simulations with = 1, = 3.0.It is evident from Fig. 3 (a) that at the highest volume fractions (Φ 4 , Φ 5 ), the energy content at high wavenumbers increases significantly.This increase in the high-wavenumber energy content is reflected in the corresponding increase in the turbulent kinetic energy dissipation rate observed in Fig. 2 (a).There is also a corresponding reduction in the energy at intermediate wavenumbers since in our DNS the TKE is held constant across all the cases due to the forcing scheme used (this behavior at intermediate wavenumbers need not occur for alternative forcing schemes).Figure 3 (b) shows that the fraction of the kinetic energy contained in the vertical motions of the flow increases at all wavenumbers, and ( )/ 33 ( ) approaches one as Φ increases.This occurs due to the transfer of the potential energy of the particles to the vertical component of the kinetic energy, with the amount of potential energy increasing with increasing Φ.The energy build in the high-wavenumber modes and increase in the share of vertical fluctuations is less pronounced in simulations with = 3.0; correspondingly, the increase in energy dissipation rate with increasing Φ is also less pronounced (Fig. 2 (c)).

Settling enhancement
We now consider the settling velocity enhancement, ( ( ), ) , and its dependence on , .In the limit → ∞ and for a finite Reynolds number flow (so that the range of fluid velocity scales is finite), to leading order, the particles settle in vertical paths and ( ( ), ) = 0.In the opposite limit → 0, ( ( ), ) = 0 also occurs because the symmetry breaking effect responsible for generating ( ( ), ) ≠ 0 (namely gravitation settling) has vanished.As a result, ( ( ), ) is expected to obtain a maximum value for some intermediate value 0 < < ∞, with the value of at which the maximum occurs depending on the other flow parameters.Tom & Bragg (2019) discussed this in detail from the perspective of the multi-scale preferential sweeping mechanism that they developed.In Fig. 4 we show the normalized settling velocity enhancement, − ( ( ), ) / , in 1WC (dashed) and the 2WC (solid) simulation regimes, as a function of , for volumefractions (a) Φ 1 and (b) Φ 4 , respectively (plots for other volume fractions are shown in the Appendix).Settling enhancements for simulations with = 0.3, 1, and 3 are shown in red, blue, and green respectively.For a given Φ we see that the settling enhancement due to turbulence generally increases with (i) increasing , for a given , and (ii) decreasing , for a fixed , both of which correspond to increasing .However, for = 0.3, a decrease in the settling-enhancement is observed when going from = 1 to = 2 for both the 1WC and 2WC simulations at Φ 1 and for the 1WC simulations at Φ 3 .For the 1WC case (the explanation for its occurrence in the 2WC case is given below), this reduction can be understood by noting that for fixed , increasing corresponds to increasing , and as already discussed, ( ( ), ) is expected to exhibit a non-monotonic dependence on .The values of and at which − ( ( ), ) / is observed to decrease must therefore correspond to values of that exceed the value of at which − ( ( ), ) / is maximum.
Figure 5  As Φ is increased, we observe increasingly large differences between the 1WC and 2WC cases, as expected (note that physically, the 1WC results should be independent of Φ. Slight variations of the 1WC results when varying Φ must therefore be due to statistical convergence effects due to the varying numbers of particles in the flow in each case).In Monchaux & Dejoan (2017) it was argued that the enhancement of particle settling speeds in the presence of 2WC compared with the 1WC case is because for a 2WC system, the settling inertial particles drag the fluid down with them which reduces the drag force on the particles and enables them to settle faster than in the 1WC case.As Φ is increased, this dragging effect on the fluid by the settling particles becomes stronger because of the increased mass-loading of the particles, and hence − ( ( ), ) / becomes larger.It was argued in Monchaux & Dejoan (2017) that this fluid dragging effect takes over from the preferential sweeping mechanism as being the dominant effect leading to − ( ( ), ) / > 0. However, in Tom et al. (2022) it was demonstrated that this is not the case, and that at least for the parameters they considered, the fluid drag effect actually works together with the preferential sweeping mechanism to further enhance − ( ( ), ) / compared to its value for the 1WC system.This will be investigated further later in this paper.
The effect of 2WC on − ( ( ), ) / is seen to become stronger for increasing (at fixed ) and/or decreasing (at fixed ), i.e. for increasing .This occurs because an increase in leads to greater potential energy of the particles, so that there is more energy available to be converted into fluid turbulent kinetic energy through the fluid dragging effect, which can lead to an increase for − ( ( ), ) / .An implication of this is that whereas − ( ( ), ) / is expected to exhibit a non-monotonic dependence on in a 1WC system, in the presence of 2WC it may monotonically increase with increasing when Φ is sufficiently large for the effects of 2WC to be strong.The data in Fig. 5 is consistent with this, except for Φ 1 and = 0.3 where a decrease of − ( ( ), ) / is observed in going from = 1 to = 2.However, this is likely because for Φ 1 , the effects of 2WC are relatively weak and so the non-monotonic dependence of − ( ( ), ) / on that is expected for a 1WC system (and which is also observed in Fig. 5(a)) would also be expected for the 2WC case.This also explains the decrease in − ( ( ), ) / observed in Fig. 4 (a) for Φ 1 and = 0.3 for the 2WC case when going from = 1 to = 2, which is not observed for the 2WC at the higher volume fraction Φ 3 in Fig. 4 (b).

Preferential sampling: , , Φ dependence
Having considered the effect of the parameters , , Φ on the enhancement of the particle settling velocity, we now turn to consider the mechanism underlying this enhancement.In Monchaux & Dejoan (2017) it had been argued that in the presence of 2WC, the enhanced settling velocities due to turbulence is no longer due to the preferential sweeping mechanism but due to the fluid dragging effect, i.e. the particles drag the fluid down with them which reduces the drag force acting on them leading to enhanced settling speeds.However, in a later study by Tom et al. (2022), some aspects of this argument was brought into question and the authors shed light on a nuanced perspective.It was shown that the preferential sweeping mechanism still operates in the presence of 2WC, but the particles drag the fluid with them as they are swept down, so that preferential sweeping and the fluid dragging effect act together   to enhance the particle settling velocities.We want to explore whether this argument remains true as the parameters are varied, especially when Φ is larger.
The preferential sweeping mechanism is associated with the preference for inertial particles to move in downward moving strain-dominated regions of the flow.To explore the role this mechanism is playing in governing the settling enhancement, we therefore first consider statistical measures that can quantify the preferential sampling of the flow field.We do this by analyzing the PDF (Q) ≡ ( ( ) − Q , where ( ) ≡ ( ( ), ) is the velocity gradient invariant measured along the inertial particle trajectory.Figure 6 (a) shows (Q) normalized using the standard deviation of (referred to as Q ), for 1WC (dashed) and 2WC (solid) simulations with = 1 particles and with = 0.3, for all volume fractions Φ.For small Φ, the 2WC results are quite similar to the 1WC results, as was also observed in Tom et al. (2022) for the case Φ 3 with = 1.As Φ is increased, the effect of 2WC on (Q) becomes apparent, with the primary impact being on the left tail, corresponding to enstrophy dominated regions of the flow.It is important to note that these differences between the 1WC and 2WC results are not only due to differences in the way that the particles are sampling the flow, but also because the statistics of the flow field itself changes due to 2WC, as was previously shown in Fig. 1 for the global statistics of ( , ).
In Fig. 6 (b), we show the normalised PDF (Q) at a constant volume fraction Φ 4 , and = 0.3 for different .Here, solid lines correspond to 2WC simulations, whereas dashed lines correspond to 1WC simulations.For the 1WC simulations, the results show that (Q) depends strongly on , indicating that the way that the particles sample the flow changes significantly as is varied.Interestingly, both the left and right tails of the PDF change significantly as is varied, whereas the results in Tom et al. (2022) for = 1 show a much weaker dependence on of the right tail of the PDF than the left tail.Since the results in Fig. 6 (b) are for = 0.3, this suggests that in some regimes of , the particles not only move away from regions of strong enstrophy, but also migrate significantly towards regions of stronger strain-rate.
To provide clearer quantitative insight into the degree to which the inertial particles are preferentially sampling the flow we can consider the the difference between the probability that the particles are in a strain-dominated or vorticity dominated regions, denoted by However, even for fluid particles (which sample the flow uniformly), P(Q > 0) − P(Q < 0) ≠ 0, and therefore simply considering P(Q > 0) − P(Q < 0) will not provide a direct measure of the degree to which the inertial particles are preferentially sampling the flow.Therefore, the measure of preferential sampling that we will use is where [P(Q > 0) − P(Q < 0)]| =0 denotes P(Q > 0) − P(Q < 0) evaluated along the trajectories of fluid particles (and since the flow is incompressible, these single-time statistics are the same as those based on measured at fixed locations in the flow).When > 1 this indicates that the inertial particles exhibit a preference to move in strain-dominated regions of the flow.
In Fig. 7 we plot for different parameter choices.Most strikingly, the results show that the preferential sampling of strain-dominated regions of the flow becomes stronger as Φ is increased.This is contrary to the argument of Monchaux & Dejoan (2017) that 2WC weakens the preferential sampling of the flow and the associated centrifuge mechanism.A potential reason for this discrepancy is that Monchaux & Dejoan (2017) based their conclusion on results for the joint PDF of : and : (see §4.1 for notation) measured along the inertial particle trajectories, and compared results for different Φ.However, as Φ is varied, this PDF is affected both by changes in the particle motion and also changes in the fluid velocity field due to 2WC.As a result, this joint PDF of : and : measured along the inertial particle trajectories does not give a direct measure of the preferential sampling of the flow by the inertial particles.To test for preferential sampling of the fluid-velocity-gradient field, one must compare the properties of the fluid velocity gradients measured along the inertial particle trajectories to those measured along fluid particle trajectories.In the Appendix, Figs. 12 and 13 show separately the quantities P(Q > 0) − P(Q < 0) and [P(Q > 0) − P(Q < 0)]| =0 , respectively.The results for P(Q > 0) − P(Q < 0) show that this quantity reduces as Φ is increased, and this is consistent with the results from (Monchaux & Dejoan 2017) which show that the joint PDF of : and : measured along the inertial particle trajectories becomes more symmetric about the line : = : as Φ is increased, i.e. that the enhanced probability to be in strain-dominated regions reduces as Φ is increased.However, the results also show that [P(Q > 0) − P(Q < 0)]| =0 decreases as Φ is increased, and this is why increases as Φ increases, even though P(Q > 0) − P(Q < 0) decreases as Φ increases.

Preferential sweeping: , , Φ dependence
The preferential sweeping mechanism not only argues that inertial particles preferentially sample strain dominated regions of the flow (which we have just demonstrated is indeed the case), but also that due to gravity, they will preferentially sample strain-dominated regions of the flow where the vertical fluid velocity points down.Tom et al. (2022) tested this by computing the quantities (solid) cases as a function of Φ at = 1 and for (a) = 0.3, (b) = 1, (c) = 3.For the 1WC results, is seen to become more negative as is reduced, indicating that the preferential sweeping is becoming stronger as is reduced over the range of considered.The results also show that < 0 is also negative, implying that there is also a contribution to the settling enhancement arising from particles in vorticity dominated regions of the flow.While this may seem surprising, it was argued in Tom et al. (2022) that this is not inconsistent with the preferential sweeping mechanism (the reader is referred to that paper for the explanation why).For = 1 and = 3 the values of and are almost identical for the 1WC and 2WC cases at the lowest Φ, but for = 0.3 there is an enhancement in the magnitudes of both and due to 2WC even for the lowest Φ.For each , however, both and become increasingly negative as Φ is increased, while always preserving Tom et al. (2022), this implies that as Φ is increased and 2WC becomes increasingly important, the preferential sweeping mechanism remains active and is the mechanism responsible for the turbulent enhancement of the particle settling speeds.The difference 2WC makes is that when the particles are swept into downward moving, strain-dominated regions of the flow, they are not only swept downwards but also drag the fluid down with them, and this is why is more negative for the 2WC cases than the 1WC cases at all , , Φ combinations considered.The bottom row of Fig. 8 illustrates the dependence of and at a fixed volume fraction Φ 4 = 7 × 10 −5 and for (e) = 0.3, (f) = 1, (g) = 3.The results show that the increased negativity of (as well as ) occurs at all considered.The results also indicate that the difference between the 1WC and 2WC values for and becomes larger as is increased, for a fixed Φ and a given .This is again because for a given , increasing corresponds to increased potential energy for the particles, and hence a greater ability to modify the fluid velocity field.
To examine more carefully the relative contributions of both and and their dependence on Φ, Fig. 9 shows the ratio / for 2WC cases as a function of Φ for different with (a) = 0.3, (b) = 1, and (c) = 3.Large values of / indicate cases where preferential sweeping plays a significant role, whereas / = 1 corresponds to the limit where the contribution to ( ( ), ) from strain-and rotation-dominated regions are equal, and therefore preferential sweeping is playing no role.The results show that the ratio decreases (except for a few cases where it initially increases, but then decreases) with increasing Φ, showing that the role of preferential sweeping becomes less important as Φ increases, in partial agreement with Monchaux & Dejoan (2017).However, it is still important for all of the cases; the ratio is greater than 1 even at the largest Φ, and for many of the cases it is closer to 2 at this Φ.An interesting and important question is whether preferential sweeping will continue to play an important role as Φ is increased even further.The results for = 0.3 seem to indicate that / becomes independent of Φ as Φ is increased, which could suggest that preferential sweeping remains important even at much larger Φ.The results for = 1, 3 do not show evidence of a regime where / becomes independent of Φ, but this could be simply because this asymptotic regime occurs at larger Φ than we have simulated.
Taken together, these results provide strong evidence that even in regimes where the mass loading is high enough for the particles to substantially modify the global fluid statistics, the preferential sweeping mechanism remains important for governing how turbulence enhances particle settling speeds.The difference compared with the 1WC case is that in the 2WC case, the particles are not merely swept around the downward moving side of vortices in the flow but they also (if the mass loading is sufficient) drag the fluid down with them in these regions as they fall.Our results also suggest that preferential sweeping may continue to be important at mass loadings much greater than those considered here, but the evidence is not conclusive.

Conclusions
We have conducted 1WC and 2WC DNS of settling, sub-Kolmogorov scale, heavy, inertial particles in homogeneous turbulence over a wide range of Froude numbers , Stokes numbers , and volume fractions Φ (the particle density was held constant, so that the mass loading Φ is proportional to Φ).For each combination of , , five different values of Φ were considered that span two orders of magnitude.We first considered the impact of 2WC on various global fluid statistics, including the properties of the strain-rate and vorticity fields and TKE dissipation rates, as well as the energy spectrum to consider the impact of 2WC on different scales in the flow.At the smallest Φ, the PDF of (the second invariant of the velocity gradient tensor) retains the characteristic skewness that is seen for the unladen flow which is associated with the vorticity field being more intermittent than the strain-rate field.However, the PDF becomes increasingly symmetric as Φ is increased, indicating that the momentum feedback from the particles suppresses regions of high vorticity while amplifying regions of high strain-rate.This effect of 2WC depends much more strongly on than on for the parameter regimes considered.For the energy spectrum we observe that as Φ is increased there is an accumulation of energy in the high wave-number modes, which also corresponds to an increase in the TKE dissipation rate.Moreover, as Φ is increased, the fraction of TKE associated with the vertical velocity fluctuations also increases, and this is due to the conversion of particle potential energy to the vertical component of TKE through the effect of the particles dragging the surrounding fluid down with them as they settle, and thereby doing work on the surrounding flow.
We then turned to consider how the settling velocity of the particles depends on , and Φ.We find that even for the smallest Φ considered (Φ = 1.5 × 10 −6 , corresponding to a mass loading of 7.5 × 10 −3 ), the mean settling speed of the particles can be appreciably larger for the 2WC case than the 1WC case for some choices of and .For both 1WC flows and 2WC flows at low Φ, the largest enhancement of the mean settling velocity compared with the Stokes settling velocity occurs for = 1 and = 0.3, and the enhancement exhibits a non-monotonic dependence on which is expected since the settling enhancement must go to zero in the two limits → ∞ and → 0. However, for 2WC flows with higher Φ this non-monotonic dependence does not occur since as is increased, the particles are increasingly effective in dragging the fluid down with them which reduces the drag force on the particle and enhances their settling velocity.Due to this, the difference between the particle settling speeds in the 1WC and 2WC flows increases with increasing .
The next step was to understand how the mechanisms governing the enhanced particle settling depends on Φ.In particular, it was to understand whether preferential sweeping remains relevant in 2WC flows when Φ becomes large enough for the global fluid statistics to be strongly affected by the particles.The preferential sweeping mechanism is based on the argument that the following two things should occur, 1) the particles should preferentially sample strain-dominated regions of the flow (which they would do even in the absence of gravity), and 2) the particles should preferentially sample strain-dominated regions where the vertical fluid velocity points down (which they would not do in the absence of gravity, if the flow is isotropic).To test the first part, we computed the PDF of measured along the trajectories of the settling inertial particles and compared them to results based on measured along fluid particle trajectories.We also compared the probability of particles being in > 0 regions compared with < 0 regions, normalized by the corresponding probability for fluid particles.The results showed that not only is preferential sampling still active in the 2WC flows, but that it actually becomes stronger as Φ is increased.To test the second part, we computed the settling velocity enhancement conditioned on the particles being in > 0 or < 0 regions.The results showed that at all , , Φ combinations considered, the strongest contribution to the enhanced particle settling velocities comes from particles in strain dominated regions of the flow, consistent with the preferential sweeping mechanism.The results do indicate that the imbalance in the contribution from > 0 and < 0 regions reduces as Φ is increased, for a given , .However, the imbalance does not vanish and remains significant at the highest Φ considered.The results also indicate that at least for low , the dominance of the contribution from strain dominated regions might persist to higher Φ than we have considered.
In conclusion, our results strongly support the idea that preferential sweeping remains the key mechanism governing the enhanced settling velocity of inertial particles in turbulent flows with 2WC, even though it was originally presented as a mechanism in the context of 1WC flows by Maxey (1987).This was demonstrated in Tom et al. (2022) for low values of Φ where the global fluid statistics are weakly affected by the particles, and only the flow in the vicinity of the particles is strongly affected.We have now demonstrated that this conclusion also holds in 2WC flows where Φ is large enough for the particles to significantly affect the global fluid statistics.
In our DNS, a relatively small Taylor Reynolds number = (10 2 ) was considered, a restriction imposed due to the computational expense of exploring a significant portion of the , , Φ parameter space.It is important in future work to consider 2WC flows with much larger where there is a large range of scales in the flow.The results could then be analyzed from the perspective of the multiscale preferential sweeping mechanism that was presented in Tom & Bragg (2019), which would provide a way to understand how flow scales of different sizes contribute to the enhanced settling velocity of inertial particles in turbulent flows.Such an attempt was done in Tom et al. (2022), however, that study also used = (10 2 ) and so the range of scales in the flow was quite small.New DNS of 2WC flows with high would provide insight into how 2WC modifies the sweeping contribution from different scales in the flow which is important for atmospheric contexts, for example, where typically (10 4 ).

Figure 1 :
Figure 1: Plot of P (Q) for 2WC simulations with (a) fixed 4 = 1, = 0.3 and varying Φ, and (b) fixed Φ 4 = 7 × 10 −5 and = 0.3 and varying .Black solid lines show P (Q) for the 1WC simulations, which is the same as for the unladen flow.Q denotes the standard deviation of .

Figure 3 :
Figure 3: (Turbulent kinetic energy spectrum ( ) versus the wave-number for 2WC simulations with = 1 particles, for different volume fractions for fixed numbers (a) = 0.3, (c) = 3.0.Ratio of the turbulent kinetic energy spectrum ( ) to the spectrum of the vertical component of the turbulent kinetic energy 33 ( ) versus the wave-number for 2WC simulations with = 1 particles for different volume fractions for fixed numbers (b) = 0.3, (d) = 3.0.
shows the volume fraction Φ dependence of the normalized settling enhancement for simulations with (a) = 0.3, and (b) = 3. Dashed and solid lines denote the Rapids articles must not exceed this page length 11 results for the 1WC and 2WC simulations, respectively, whereas different line colors denote different .The results show that even for the smallest volume fraction considered (Φ 1 ), the momentum-coupling results in a finite, additional settling enhancement compared to the 1WC case, indicating that the particles are modifying the flow field in their vicinity.As shown in §4.1, for Φ 1 the global fluid statistics are almost identical to those of the unladen flow and this is because for low Φ, modifications to the flow in the vicinity of the particles makes a negligible contribution to the global fluid behavior.

Figure 6 :
Figure 6: Results for (Q) in 1WC and 2WC simulations, along particle trajectories (a) for St 1 = 1, = 0.3 showing Φ dependence, and (b) at a constant volume fraction Φ 4 = 7 × 10 −5 , and = 0.3 showing dependence.Dashed and solid lines mark the trends for 1WC and 2WC simulations respectively.Q denotes the standard deviation of .

Table 1 :
Flow parameters in DNS for the unladen flow.

Table 2 :
Flow parameters in DNS of 2WC simulations for different volume fraction Φ, Stokes number and Froude number .The Stokes number and Froude number are defined with respect to the unladen DNS statistics.(Continued) Continued on next page Table 2: Flow parameters in DNS of 2WC simulations for different volume fraction Φ, Stokes number and Froude number .The Stokes number and Froude number are defined with respect to the unladen DNS statistics.(Continued)