2. Problem formulation

We consider a dilute suspension (i.e. the volume fraction $\varPhi _v$ is small enough to ignore particle collisions) of small ($d_p/\eta \ll 1$, where $d_p$ is the particle diameter and $\eta$ is the Kolmogorov length scale), dense ($\rho _p/\rho _f \gg 1$, where $\rho _p$ is the particle density and $\rho _f$ is the fluid density), spherical, inertial particles settling under the force of gravity. Such assumptions are often considered to be suitable for modelling droplet transport in atmospheric clouds (Shaw Reference Shaw2003; Grabowski & Wang Reference Grabowski and Wang2013), and dust transport in the atmosphere (Richter & Chamecki Reference Richter and Chamecki2018).

Assuming that the particle Reynolds number is small, the evolution equation for particles in this regime is given by (Maxey & Riley Reference Maxey and Riley1983)

(2.1)\begin{equation} \ddot{\boldsymbol{x}}^p(t)\equiv\dot{\boldsymbol{v}}^p(t)=\frac{1}{\tau_p} [\boldsymbol{u}(\boldsymbol{x}^p(t),t)-\boldsymbol{v}^p(t)]+\boldsymbol{g},\end{equation}

where $\boldsymbol {x}^p(t)$ and $\boldsymbol {v}^p(t)$ denote the particle position and velocity, respectively, $\boldsymbol {u}(\boldsymbol {x}^p(t),t)$ is the undisturbed fluid velocity at the particle position, $\tau _\eta$ is the Kolmogorov time scale and $\boldsymbol {g}$ denotes the gravitational acceleration. The particle response time $\tau _p$ is given by

(2.2)\begin{equation} \tau_p = \frac{1}{18}\frac{\rho_p/\rho_f}{\nu} d_p^2, \end{equation}

where $\nu$ is the dynamic viscosity.

The particles settle in a statistically stationary, homogeneous turbulent flow, which is governed by the incompressible Navier–Stokes equation (written here in rotational form)

(2.3a,b)\begin{equation} \partial_t \boldsymbol{u} =\boldsymbol{u}\times\boldsymbol{\omega} -\boldsymbol{\nabla}\left(\frac{P}{\rho_f} + \frac{u^2}{2}\right) + \nu\nabla^2\boldsymbol{u} + \boldsymbol{F} + \boldsymbol{C},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0. \end{equation}

Here, $\boldsymbol {u}(\boldsymbol {x},t)$ and $\boldsymbol {\omega }\equiv \boldsymbol {\nabla }\times \boldsymbol {u}$ are the fluid velocity and vorticity fields, respectively, $P(\boldsymbol {x},t)$ is the pressure, $\boldsymbol {C}(\boldsymbol {x},t)$ is the total momentum feedback of the particles on the carrier fluid and $\boldsymbol {F}(\boldsymbol {x},t)$ is the large-scale forcing required to maintain steady-state turbulence. In our direct numerical simulations (DNSs) we will consider the case where $\boldsymbol {F}$ isotropically forces the flow, such that, when $\boldsymbol {g}=\boldsymbol {0}$, the flow is statistically isotropic, but it is anisotropic for $\boldsymbol {g}\neq \boldsymbol {0}$ due to the momentum coupling term $\boldsymbol {C}$.

We consider the particle motion in both the 1WC scenario where we set $\boldsymbol {C}=\boldsymbol {0}$, and the 2WC scenario where $\boldsymbol {C}\neq \boldsymbol {0}$. The particle feedback $\boldsymbol {C}$ is the sum of the hydrodynamic forces generated by all particles at location $\boldsymbol {x}^p(t)=\boldsymbol {x}$

(2.4a,b)\begin{equation} \boldsymbol{C}(\boldsymbol{x},t)=\frac{1}{\rho_f}\sum_i \boldsymbol{f}_i^p(t) \delta(\boldsymbol{x}-\boldsymbol{x}^p_i(t)),\quad \boldsymbol{f}_i^p(t) =-\frac{m_p}{St\tau_\eta} [\boldsymbol{u}(\boldsymbol{x}^p_i(t),t)-\boldsymbol{v}^p_i(t)], \end{equation}

where $m_p$ is the mass of each particle, $\boldsymbol {x}^p_i$ and $\boldsymbol {v}^p_i$ are the position and velocity of the $i$th particle and $\boldsymbol {f}^p_i$ is the reaction force generated by the $i$th particle. As discussed in Tom et al. (Reference Tom, Carbone and Bragg2022), when $\boldsymbol {C}\neq \boldsymbol {0}$, settling particles will generate a mean flow in the direction of $\boldsymbol {g}$. 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 (Reference Maxey and Patel2001), Bosse et al. (Reference Bosse, Kleiser and Meiburg2006), Monchaux & Dejoan (Reference Monchaux and Dejoan2017) and Rosa et al. (Reference Rosa, Kopeć, Ababaei and Pozorski2021).

The ensemble-averaged mean particle settling speed can be obtained from (2.1), assuming statistical homogeneity and stationarity

(2.5)\begin{equation} \langle{v}_z^p(t)\rangle=\langle{u}_z(\boldsymbol{x}^p(t),t)\rangle-St\tau_\eta{g}.\end{equation}

Here, $St \tau _\eta g$ is the Stokes settling speed (the settling speed the particle would have in a quiescent fluid) and $\langle \boldsymbol {u}(\boldsymbol {x}^p(t),t) \rangle$ 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 $\langle {u}_z(\boldsymbol {x},t)\rangle ={0}$, the average along the particle trajectory $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle$ need not be zero because inertial particles will not in general sample the flow field ergodically. As a result of this, through the contribution $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle$ turbulence can lead to either enhanced settling ($\langle {v}_z^p(t)\rangle < -St\tau _\eta {g}$) or hindering ($\langle {v}_z^p(t)\rangle > -St\tau _\eta {g}$). 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 Reference Wang and Maxey1993; Yang & Shy Reference Yang and Shy2005; Bec et al. Reference Bec, Homann and Ray2014; Good et al. Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014; Ireland, Bragg & Collins Reference Ireland, Bragg and Collins2016; Rosa et al. Reference Rosa, Parishani, Ayala and Wang2016; Nemes et al. Reference Nemes, Dasari, Hong, Guala and Coletti2017; Petersen et al. Reference Petersen, Baker and Coletti2019; Li et al. Reference Li, Lim, Berk, Abraham, Heisel, Guala, Coletti and Hong2021a,Reference Li, Abraham, Guala and Hongb; Tom et al. Reference Tom, Carbone and Bragg2022).

3. Direct numerical simulations

3.2. Simulation parameters

Table 1 summarises the parameters of our unladen DNS, where $\mathcal {L}$ is the domain length, $N$ is the number of collocation points (the grid size is thus $N^3$) and $\nu$ 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 $\langle \epsilon \rangle$, the root mean square (r.m.s.) fluctuating velocity $u'\equiv \sqrt {2\mathcal {K}/3}$ (where $\mathcal {K}$ is the turbulent kinetic energy), the Kolmogorov velocity and length scales, $u_{\eta }$ and $\eta$, respectively, the integral length scale $L$ and the Taylor microscale Reynolds number $Re_{\lambda } \equiv u'\lambda /\nu$, where $\lambda$ is the Taylor length scale.

Table 1.

Flow parameters in DNS for the unladen flow. Here, $Re_{\lambda } \equiv u'\lambda /\nu = 2k/\sqrt {5/3\nu \langle \epsilon \rangle }$ is the Taylor microscale Reynolds number, $\lambda$ is the Taylor microscale, $\mathcal {L}$ is the domain length, $N$ is the number of grid points in each direction, $\nu$ is the fluid kinematic viscosity, $\langle \epsilon \rangle$ is the mean turbulent kinetic energy dissipation rate, $L$ is the integral length scale, $\eta \equiv (\nu ^{3}/\langle \epsilon \rangle )^{1/4}$ is the Kolmogorov length scale, $u' \equiv \sqrt {2\mathcal {K}/3}$ is the r.m.s. of the fluctuating fluid velocity, $\mathcal {K}$ is the turbulent kinetic energy, $u_{\eta }$ is the Kolmogorov velocity scale, $\tau _L \equiv L/u'$ is the large-eddy turnover time, $\tau _{\eta }$ is the Kolmogorov time scale, $k_{max} = \sqrt {2N/3}$ is the maximum resolved wavenumber and $dt$ is the time step. The small-scale resolution, $k_{max}\eta$ and the total flow kinetic energy measured by $u'$ 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\tau _L$.

For the 1WC and 2WC simulations, the Stokes number $St\equiv \tau _p/\tau _\eta$ and the Froude number $Fr\equiv u_\eta /(\tau _\eta g)$ are defined with respect to the unladen fluid statistics, and we consider Stokes numbers $St=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 $Fr=0.3$, $1$ and $3$ (referred to as $Fr_1$, $Fr_2$ and $Fr_3$, respectively). The particle mass loading $\varPhi _m$ is related to the volume fraction $\varPhi$ through

(3.1) \begin{equation} \varPhi_m=\frac{\rho_p}{\rho_f}\varPhi=\frac{\rho_p}{\rho_f} \frac{{\rm \pi} d_p^3}{6} \frac{N_p}{L^3}, \end{equation}

and we consider particles with constant density $\rho _p/\rho _f= 5000$, the same as that used in Bosse et al. (Reference Bosse, Kleiser and Meiburg2006), Monchaux & Dejoan (Reference Monchaux and Dejoan2017) and Tom et al. (Reference Tom, Carbone and Bragg2022). Since $\rho _p/\rho _f$ is fixed, varying $\varPhi _m$ corresponds to varying $\varPhi$, and therefore we will often refer to the effect of varying $\varPhi$ even though strictly speaking the momentum coupling depends on $\varPhi _m$ and not simply $\varPhi$. For each choice of $St$ and $Fr$, we consider $\varPhi = 1.5\times 10^{-6}$, $7\times 10^{-6}$, $1.5\times 10^{-5}$, $7\times 10^{-5}$ and $1.5\times 10^{-4}$, which we refer to as $\varPhi _1,\varPhi _2,\varPhi _3,\varPhi _4,\varPhi _5$, respectively. However, due to computational resources, we were unable to complete the run with $St=2, Fr=0.3, \varPhi =1.5\times 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 $Re_{\lambda } \simeq 88$. The present study significantly expands on the study of Tom et al. (Reference Tom, Carbone and Bragg2022), where only $Fr=1$ and $\varPhi = 1.5\times 10^{-5}$ were considered, for the same $St$ and $Re_\lambda$ values.

5. Comparison with previous work

So far, we have considered the impact of 2WC on particle settling and certain flow statistics by comparison with the corresponding results from our DNS for 1WC simulations. We also discussed the ways in which our results compare qualitatively with those of Monchaux & Dejoan (Reference Monchaux and Dejoan2017). In this section we turn to considering a quantitative comparison of our DNS results with selected prior experiments and numerical studies that have parameters closely aligned with the current simulations. This comparison includes the studies by Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002), Bosse et al. (Reference Bosse, Kleiser and Meiburg2006), Monchaux & Dejoan (Reference Monchaux and Dejoan2017), Berk & Coletti (Reference Berk and Coletti2021) and Hassaini & Coletti (Reference Hassaini and Coletti2022). Table 2 presents a summary of the different simulations and experiments in these studies. It is important to note that, despite their significance, recent works by Petersen et al. (Reference Petersen, Baker and Coletti2019), Rosa et al. (Reference Rosa, Kopeć, Ababaei and Pozorski2022) and Hassaini et al. (Reference Hassaini, Petersen and Coletti2023) were not included in this comparison due to the lack of sufficient parameter overlap. For example, in the study by Petersen et al. (Reference Petersen, Baker and Coletti2019), $Re_\lambda \geq 300$, in Hassaini et al. (Reference Hassaini, Petersen and Coletti2023), the flow Reynolds number is higher $Re_\lambda \gtrsim 150$ and we have poor overlap in $Fr\sim 0.1, 0.7$ and, in Rosa et al. (Reference Rosa, Kopeć, Ababaei and Pozorski2022), we have poor overlap in the Stokes number $St$ and Rouse number $R=u_\eta Sv/u^\prime$.

Table 2.

Details of previous experiments (denoted by E) and numerical simulations (denoted by S). Here, ${\rm AL}\,{\rm E1}$ and ${\rm AL}\,{\rm E2}$ denote the measurements of Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002) in isotropic, decaying grid turbulence (squares and circles, respectively, in figure 14(a) of ${\rm AL}$); ${\rm YS}$ indicates the measurements of Yang & Shy (Reference Yang and Shy2005) in stationary, homogeneous, isotropic turbulence generated by a pair of counter-rotating fans (refer to figure 16(a) for data); ${\rm BM}\,{\rm S1}\unicode{x2013}{\rm S3}$ denote the results of DNS from Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) using a stochastic forcing scheme and a computational particle approach (in the Eulerian–Lagrangian framework); see table II (${\rm BM}\,{\rm S1}$) and table V, figure 14 (diamond for ${\rm BM}\,{\rm S2}$ and triangle for ${\rm BM}\,{\rm S3}$); ${\rm GW}$ indicates the studies presented in Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014): ${\rm GW}\,{\rm E1}$ denotes the experiment of particle settling in soccer-ball apparatus with 32 loudspeaker jets (labelled E1 in ${\rm GW}$ table 1), whereas ${\rm GW}\,{\rm S1}$ marks the results obtained from DNS using a deterministic forcing scheme (simulation labelled S2 in ${\rm GW}$ table 1 interpolated to the same $Fr$ as ${\rm GW}\,{\rm E1}$); see figure 3(a) in ${\rm GW}$ for data; ${\rm MD}$ refers to the numerical studies presented in Monchaux & Dejoan (Reference Monchaux and Dejoan2017) using DNS. The simulations ${\rm S1}$ and ${\rm S2}$ correspond to simulations, for Rouse number $R=\tau _pg/u^\prime =0.25$, with $St=0.36$ and $St=1$, respectively; see table II for additional details; ${\rm BC}\,{\rm E1},2$ indicate the results from the experiments in a turbulence chamber (Berk & Coletti Reference Berk and Coletti2021) with facing jet arrays at two different $Re_\lambda$ (see table 2 and figure 8 in ${\rm BC}$ for data); ${\rm HC}\,{\rm E1},2$ denote the experiments in a similar set-up as ${\rm BC}$, but with the jet arrays activated in a randomised sequence (Hassaini & Coletti Reference Hassaini and Coletti2022); refer to table 1 and figure 3 in ${\rm HC}$ for data.

Figure 10 compares the ratios of the r.m.s. fluid velocities in different directions from our 2WC DNS simulations with data from the DNS of Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) and experiments of Hassaini & Coletti (Reference Hassaini and Coletti2022), hereafter ${\rm BM}$ and ${\rm HC}$, respectively. See table 2 for the details of these studies. Panel (a) shows the ratio of the horizontal to the vertical r.m.s. velocities, $u^\prime _1/u^\prime _3$ (solid lines) and $u^\prime _2/u^\prime _3$ (dashed lines), as a function of the volume fraction $\varPhi$. The results of ${\rm HC}\,{\rm E1}$ (${\rm HC}\,{\rm E2}$) are denoted by black circles (diamonds), and those of ${\rm BM}\,{\rm S1}$ are denoted by black squares. Data from the simulations with parameters closest to ${\rm HC}\,{\rm E1}$ and ${\rm HC}\,{\rm E2}$ are indicated by blue circles and red diamonds, respectively. We see that the r.m.s. velocity in the horizontal direction decreases compared with that in the vertical direction as $\varPhi$ increases. This is because the amount of particle potential energy that is converted to the vertical component of the TKE though the fluid drag mechanism increases with volume fraction $\varPhi$. This observed behaviour is opposite to the trend observed in the experiments of ${\rm HC}$, but in agreement with the DNS results of ${\rm BM}$. Panel (b) shows the ratio of the horizontal r.m.s. velocities for 2WC simulations to that of the corresponding 1WC simulations, $u^\prime _1/u^\prime _{1,\varPhi =0}$ (solid lines) and $u^\prime _2/u^\prime _{2,\varPhi =0}$ (dashed lines), as a function of the volume fraction $\varPhi$. We observe that these ratios decrease as $\varPhi$ increases, which is again in keeping with ${\rm BM}$ but in contrast to ${\rm HC}$.

Figure 10.

Ratio of the (a) horizontal to vertical turbulence intensities $u^\prime _1/u^\prime _3$, $u^\prime _2/u^\prime _3$ and (b) horizontal turbulence intensities for 2WC in unladen simulations $u^\prime _1/u^\prime _{1,\varPhi =0}$, $u^\prime _2/u^\prime _{2,\varPhi =0}$ for DNS (coloured symbols) and data from previous work (in black). Experimental results from Hassaini & Coletti (Reference Hassaini and Coletti2022) figure 3(a) are shown in black circles ($St=0.3$) and black diamonds ($St=2.6$), denoted by ${\rm HC}\,{\rm E1}$ and ${\rm HC}\,{\rm E2}$, respectively. Results from Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) table II, denoted by ${\rm BM}\,{\rm S1}$, are given by black squares. Solid and dashed lines for the current and previous numerical results simply correspond to data along to two different horizontal directions. Details of the previous works are given in table 2.

A plausible explanation for the discrepancy between the DNS (both ours and ${\rm BM}$) and the experimental results of ${\rm HC}$ is that the particles introduce a certain level of anisotropy into the flow in the measurement area of the experiments which is not present in the DNS. In particular, since inertial particles carry a memory of their interaction with the flow (Bragg & Collins Reference Bragg and Collins2014a,Reference Bragg and Collinsb), they may carry into the measurement area a memory of their interaction with the large-scale flow which is anisotropic in the experiment even in the absence of the particles due to the jets used by ${\rm HC}$ to force the flow. This would then cause the particles to transfer momentum to the flow in the measurement region in an anisotropic manner, even though the measurement region is small enough such that in the absence of the particles the flow would be approximately locally isotropic. By contrast, our DNS and that of ${\rm BM}$ force the large scales of the flow in an isotropic manner, and all anisotropy in the flow emerges solely due to the particle settling.

In our DNS, TKE dissipated at the small scales is re-injected into the flow at the large scales in an isotropic manner. This means that, for the subset of Fourier modes that are forced, the flow anisotropy that would otherwise arise for these Fourier modes is effectively suppressed. However, this is not likely to be the cause of the discrepancy between the DNS and experiments noted above for the following reason: in the absence of forcing, some of the vertical kinetic energy injected to the flow due to the particle settling would be transferred to the horizontal fluid motions due to the pressure-strain redistribution effect. If anything, our isotropic forcing would artificially enhance this redistribution of energy from the vertical to horizontal flow motions, and therefore cause the ratios $u^\prime _1/u^\prime _3$ and $u^\prime _2/u^\prime _3$ to be artificially enhanced by the isotropic forcing, not reduced.

Figure 11 compares the TKE dissipation rate (normalised by 1WC values) of ${\rm BM}\,{\rm S1}$ (black), which corresponds to $St=1$ and $Fr=1$, with the corresponding values from the current 2WC DNS (red) with the closest parameters. In contrast to ${\rm BM}\,{\rm S1}$, which predicts an initial decrease of $\epsilon$ with $\varPhi$ before increasing again, the current simulations predict a monotonic increase in the dissipation rate. One possible reason for this difference is due to differences in the Reynolds numbers of the flow, since ours has $Re_\lambda \approx 88$ while that of ${\rm BM}\,{\rm S1}$ has a much lower value of $Re_\lambda \approx 42$. Indeed, figure 3 shows that, for some choices of $St,Fr$, $\epsilon$ is also observed to vary non-monotonically with $\varPhi$ in our DNS, and variations in $Re_\lambda$ as well as $St,Fr$ may also lead to similar regimes. Another possible explanation is due to the different forcing methods used in the DNS (discussed more at the end of this section), with our DNS using a deterministic forcing that maintains a constant flow TKE, while Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) use a stochastic method that injects a statistically constant amount of energy into the flow. An alternate cause of these disparities could be the different interpolation schemes used to compute the Stokes drag and the particle-momentum feedback (2.4a,b); Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) implemented a tri-linear interpolation method whereas our study utilised a seventh-order B-spline interpolation. Exactly how these different methods would impact the results in figure 11 is not obvious, however.

Figure 11.

Comparison of the TKE dissipation rate (normalised by 1WC values) for the current simulations (red circles) and Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) table II (black diamonds), denoted by ${\rm BM}\,{\rm S1}$, vs volume fraction $\varPhi$.

Figure 12 compares the normalised Lagrangian PDF $\mathscr {P}(\mathcal {Q})$ for the current simulations with experimental results of Berk & Coletti (Reference Berk and Coletti2021). Note that the $\mathcal {Q}$ defined here is the negative of that used in ${\rm BC}$, and that the PDFs have been normalised by their standard deviation to facilitate the comparison. The black solid and dashed lines show, for reference, the Eulerian PDFs from the 1WC (or unladen) simulations and the experiment ${\rm BC}$, respectively. Black circles and diamonds correspond to measurements along Lagrangian particle trajectories in ${\rm BC}\,{\rm E1}$ ($St=1.2$) and ${\rm BC}\,{\rm E2}$ ($St=2.1$), respectively, for $\varPhi \approx 10^{-7}$; details of these experiments are given in table 2. The results from our 1WC and 2WC simulations with parameters closest to the experiments are shown by the (dashed for 1WC and solid for 2WC) yellow ($St=1$) and purple ($St=2$) lines, respectively. It is noteworthy that the Eulerian PDF from ${\rm BC}$ is quite close to the 1WC (effectively unladen) DNS results for the left tail (corresponding to rotation-dominated regions) but predicts a much larger probability of straining events (see the heavier right tail), compared with the DNS. One possible explanation for this significant disagreement is that the difference arises due to the $Re_\lambda$ of the DNS being much lower ($Re_\lambda \approx 88$) than that in the experiments ($Re_\lambda \geq 280$). However, this does not seem likely because it is known that the vorticity field is more intermittent than the strain-rate field (Tsinober Reference Tsinober2001), and so the effect of $Re_\lambda$ on the left tail of the PDF should be stronger than that on the right tail, yet the DNS and experiments agree well for the left tail. A more plausible explanation for the difference is that it results from large-scale shear in the experiments due to the jets, which may have a significant effect on the velocity-gradient field during extreme events in the flow where large- and small-scale motions can exhibit significant coupling.

Figure 12.

Comparison of the current DNS results for normalised Lagrangian PDFs $\mathscr {P}(\mathcal {Q})$ (symbols) in 1WC and 2WC simulations and the experiments of Berk & Coletti (Reference Berk and Coletti2021), hereafter ${\rm BC}$. Black solid line corresponding to the 1WC simulations, which is the same as for the unladen flow, and the black dashed line corresponds to the red line from figure 8 of ${\rm BC}$ (corresponding to the unladen fluid statistics), are shown for reference. Black circles and diamonds correspond to measurements along Lagrangian particle trajectories in ${\rm BC}\,{\rm E1}$ and ${\rm BC}\,{\rm E2}$, respectively, for $\varPhi \approx 10^{-7}$. See table 2 for the details of these experiments. Here, $\sigma _{\mathcal {Q}}$ denotes the standard deviation of $Q$.

In figure 13 we present a comparison of our 1WC and 2WC DNS results for the particle settling velocity enhancement with previous results, both experimental and numerical. Dashed and solid lines represent the current DNS results with 1WC and 2WC, respectively, whereas dotted lines and open symbols (no lines) represent previous experimental and numerical results, respectively. Table 2 summarises the notation for these experiments and simulations. Plots at a fixed volume fraction $\varPhi \approx 10^{-6}$, vs $St$, are shown in figure 13(a) for $Fr \approx 0.3$. The current 2WC results (solid green line) agree well with the experimental results (red dotted line) of ${\rm GW}\,{\rm E1}$ for $St \leq 0.7$, while over-predicting the settling enhancement for $St\geq 1$ by 30 %–50 %. In this parameter range, we see that the 1WC DNS result is actually in better agreement with the experimental data. It is worth pointing out that the current 1WC DNS results differ in significant ways from those of Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014) in blue (${\rm GW}\,{\rm S1}$), even though they are generated using essentially the same DNS code based on Ireland et al. (Reference Ireland, Vaithianathan, Sukheswalla, Ray and Collins2013). In Tom et al. (Reference Tom, Carbone and Bragg2022) it was discussed that very long time averages (such as we use in our DNS) are needed for the particle settling velocity enhancement results to converge, and that the time averaging window used in Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014) was relatively short. This is most likely the explanation for the difference in the results, and highlights the importance of ensuring statistical convergence of the particle settling velocities.

Figure 13.

Comparison of 1WC (dashed lines) and 2WC (solid lines) DNS results for settling velocity enhancement (normalised by $u^{\prime }$), as a function of (a,b) $St$ for different $\varPhi$ and $Fr$, (c) $\varPhi$ for a constant $St$ and $Fr$, with previous studies. See table 2 for details of experiments and numerical simulations. Experimental and numerical results from previous studies are shown by dotted lines and open symbols, respectively.

In figure 13(b), the dashed and the solid green lines show the settling enhancement vs $St$, for the current 1WC and 2WC simulations (at $\varPhi _3=1.5\times 10^{-5}$ and $Fr_2=1$), respectively. The dotted green line denotes the results from experiment ${\rm AL}\,{\rm E1}$ (see table 2). The black diamond marks the corresponding prediction of ${\rm BM}\,{\rm S2}$. Our 2WC results align well with the measured settling velocity enhancement for $St \geq 0.7$; however, the agreement is quantitatively poor for smaller $St$. Figure 13(b) also compares the settling velocity enhancement predictions from our current simulations with those from previous experiments at volume fraction $\varPhi =7\times 10^{-5}$. The dotted and the solid yellow lines correspond, respectively, to the experiment ${\rm AL}\,{\rm E2}$ and the current 2WC DNS at $\varPhi _4$ and $Fr_2$. Note that both the experiments by ${\rm AL}$ correspond to $Fr \approx 1.6$, whereas the numerical simulations are for $Fr=1$, and the impact of this difference is not known. The corresponding data for the 1WC simulation and the data from ${\rm BM}\,{\rm S3}$ are given by the dashed green line (1WC trend for $\varPhi _4$ is the same as $\varPhi _3$ and hence not shown separately) and the black star, respectively. Compared with ${\rm BM}\,{\rm S3}$, the current results are in better agreement with the experimental data, however, the estimated settling enhancements are consistently smaller than the measured values for $St<2$. In particular, the predicted settling enhancement is approximately three times smaller than the observed values for the $St=0.3$. Yang & Shy (Reference Yang and Shy2005), hereafter ${\rm YS}$, and ${\rm BM}$ discuss, in detail, the shortcomings of the set-up in ${\rm AL}$, which potentially explain the larger values of the settling enhancement in ${\rm AL}\,{\rm E1},2$. However, it is unclear why the largest discrepancies are observed for small $St$; disagreement with our data is only approximately $20\,\%$ for $St \approx 1$.

In figure 13(b), the dashed, solid and the dotted brown lines show the results for the settling enhancement, vs $St$, for the current 1WC and 2WC DNS simulations (at $\varPhi _4$, $Fr_1$) and the results of ${\rm YS}$ figure 16, respectively. Our 2WC results for this quantity are much larger than those of ${\rm YS}$; ${\rm BM}$ also have DNS results for similar parameters (not shown here) and their results also disagreed considerably with the corresponding results of ${\rm YS}$. Figure 13(c) compares the $\varPhi$ dependence of the settling velocity enhancement results from the current 1WC (dashed lines) and 2WC (solid lines) simulations with those from ${\rm BM}$ (black square) and ${\rm MD}$. Refer to Monchaux & Dejoan (Reference Monchaux and Dejoan2017, table II) for the details of the simulations ${\rm MD}\,{\rm S1}$ and ${\rm MD}\,{\rm S2}$. The numerical simulations ${\rm S1}$ and ${\rm S2}$ from ${\rm MD}$ correspond to $St=0.36$ (orange circle) and $St=1$ (red square), at a constant Rouse number $R=\tau _pg/u^\prime =0.25$; see solid black line with downward triangle and solid blue line with circle in ${\rm MD}$ figure 3(b) for the data. For the parameter choices of our DNS results shown, the Rouse number is $R \simeq 0.29$. Overall we see quite good agreement of the current 2WC DNS results with ${\rm BM}\,{\rm S1}$, ${\rm MD}\,{\rm S1}$ and ${\rm MD}\,{\rm S2}$.

Summing up the results presented in this section, we find good qualitative agreement and generally reasonable quantitative agreement among the DNS results that have been considered. The parameters in the different DNS studies are not exactly the same (we compared those with the closest parametric overlap), and therefore exact quantitative agreement is not expected. Other than differences in the parameters, the DNS differ in how the momentum coupling term is implemented numerically as well as how the flow is forced. The particle-momentum feedback on the computational grid in the current work follows that from Carbone et al. (Reference Carbone, Bragg and Iovieno2019) and Tom et al. (Reference Tom, Carbone and Bragg2022). On the other hand, Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) and Monchaux & Dejoan (Reference Monchaux and Dejoan2017) implement a tri-linear interpolation method and the particle-in-cell method, respectively, for the momentum feedback term. The forcing scheme we use follows that from Tom et al. (Reference Tom, Carbone and Bragg2022) and the work of Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014). Monchaux & Dejoan (Reference Monchaux and Dejoan2017) also use a forcing scheme whereby the dissipated energy is reintroduced back into the system. Instead of maintaining a constant TKE, Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) used a forcing method that maintains a constant TKE injection rate following Eswaran & Pope (Reference Eswaran and Pope1988). It should be noted that the classical argument that the energy dissipation and particle dynamics are mainly governed by the small-scale dynamics does not apply, since the turbulent statistics change according to how the large-scale forcing is implemented. To determine the influence of these different numerical methodologies on the results for particle settling, the research community would benefit greatly from a future study where results from different DNS using identical parameters are compared, similar to what was done in Marchioli et al. (Reference Marchioli, Soldati, Kuerten, Arcen, Tanière, Goldensoph, Squires, Cargnelutti and Portela2008) to compare different DNS codes for inertial particle transport in wall-bounded turbulence. Comparisons of DNS with experimental results showed that, for some quantities, the DNS and experiments show good qualitative agreement but quite poor quantitative agreement, while for other quantities they even differ quantitatively. However, for the particle settling enhancement results, different experiments gave very different results. It is not clear whether these differences occur due to physical differences in the experiments (different parameters, flow types etc.) or due to limitations in the measurement and/or data post-processing methodologies. Therefore, just as it is important in future work for the results from different DNS codes to be carefully compared in order to understand their differences, so also there needs to be a concentrated effort to understand differences between published experimental results. Only when this has been done will it be possible to make progress on understanding exactly which modelling approximations made in the DNS are principally responsible for the discrepancies between the DNS and experimental results.