2.1. Background

We consider the settling of small ($d_p/\eta \ll 1$, where $d_p$ is the particle diameter and $\eta$ is the Kolmogorov length scale), heavy ($\rho _p/\rho _f \gg 1$, where $\rho _p$ is particle density and $\rho _f$ is fluid density), spherical, mono-disperse inertial particles in a statistically stationary, homogeneous (and initially isotropic) turbulent flow. In the regime of a linear drag force, the particle equation of motion reduces to (see Maxey & Riley Reference Maxey and Riley1983)

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

where $\boldsymbol {x}^p(t),\boldsymbol {v}^p(t)$ are the particle position and velocity vectors, respectively, $\boldsymbol {u}(\boldsymbol {x}^p(t),t)$ is the undisturbed fluid velocity at the particle position and $\boldsymbol {g}$ is the gravitational acceleration vector. Even though idealized, this simplified equation of motion is widely used for describing particle motion in atmospheric flows (e.g. Grabowski & Wang Reference Grabowski and Wang2013), and even for this equation of motion, many unanswered questions remain concerning the settling of particles in turbulence. Moreover, our DNS will be confined to $St\leq 2$ and $Sv\leq 2$ for which we expect that the linear drag approximation will be reasonable.

In a 2WC system, a particle exerts a force on the fluid at its location equal to

(2.2)\begin{equation} \boldsymbol{f}^p(t) ={-}\frac{m_p}{St\,\tau_\eta}(\boldsymbol{u}(\boldsymbol{x}^p(t),t)-\boldsymbol{v}^p(t)), \end{equation}

which is equal and opposite to the drag force exerted by the fluid on the particles that is described by (2.1), and $m_p=(4/3){\rm \pi} (d_p/2)^3\rho _p$ is the mass of the spherical particle. The addition of the feedback force, $\boldsymbol {f}^p(t)$, breaks the isotropic symmetry of an otherwise isotropic turbulent flow, but the fluid is still statistically homogeneous since the gravitational body force is independent of position.

Statistical stationarity of the system implies that, in the vertical direction defined by $\boldsymbol {e}_z$, the ensemble average of (2.1) is (with $\boldsymbol {g}=-g\boldsymbol {e}_z$)

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

Note that we used $\boldsymbol {g}=g\boldsymbol {e}_z$ in Tom & Bragg (Reference Tom and Bragg2019), but have decided to switch to this more standard choice of coordinate system in this paper. Equation (2.3) shows that the average particle velocity may differ from the Stokes settling velocity $-St\,\tau _\eta {g}$ only if $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle \neq {0}$ (we consider flows with Eulerian average $\langle u_z(\boldsymbol {x},t)\rangle =0$, such that $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle \neq {0}$ can only arise due to non-trivial effects in the flow). Several numerical studies for both 1WC (Wang & Maxey Reference Wang and Maxey1993; Bec et al. Reference Bec, Homann and Ray2014; Ireland et al. Reference Ireland, Bragg and Collins2016b; Rosa et al. Reference Rosa, Parishani, Ayala and Wang2016) and 2WC systems (Bosse et al. Reference Bosse, Kleiser and Meiburg2006; Dejoan Reference Dejoan2011; Monchaux & Dejoan Reference Monchaux and Dejoan2017; Rosa et al. Reference Rosa, Kopeć, Ababaei and Pozorski2021) have indeed shown that $\langle {v}_z^p(t)\rangle \neq -St\,\tau _\eta {g}$, implying $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle \neq {0}$. Experimental studies in the laboratory (Aliseda et al. Reference Aliseda, Cartellier, Hainaux and Lasheras2002; Yang & Shy Reference Yang and Shy2005; Good et al. Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014; Petersen et al. Reference Petersen, Baker and Coletti2019) and more recently, settling velocity measurements in the atmosphere (Nemes et al. Reference Nemes, Dasari, Hong, Guala and Coletti2017; Li et al. Reference Li, Lim, Berk, Abraham, Heisel, Guala, Coletti and Hong2021a,Reference Li, Abraham, Guala and Hongb), both of which inherently have 2WC effects, have also observed an enhancement of settling velocity compared with the Stokes settling velocity. Hence, there is a need to explore the physical mechanism(s) responsible for generating $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <{0}$ and whether there are different mechanisms at play for the 1WC and 2WC systems.

2.2. Brief summary of the theoretical framework

For the 1WC system, Maxey (Reference Maxey1987) developed a theoretical framework to explain how $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <{0}$ can occur even in flows where the Eulerian flow average is $\langle {u}_z(\boldsymbol {x},t)\rangle ={0}$. Essentially, the explanation is that particles with inertia do not uniformly sample the underlying fluid velocity field, and that gravity leads to a bias for inertial particles to accumulate in regions of the flow where $\boldsymbol {e}_z\boldsymbol {\cdot } \boldsymbol {u}<0$. However, the analysis of Maxey (Reference Maxey1987) was restricted to $St\ll 1$. In a recent paper (Tom & Bragg Reference Tom and Bragg2019), we developed a theoretical framework for considering the behaviour of $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle$ for arbitrary $St$, and for revealing the scales of motion that contribute to $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle$ being finite and negative. Here, we briefly summarize the steps of the theoretical analysis and comment on how the analytical result can be used to understand how 2WC might modify the physical mechanisms governing the particle settling velocity.

The key ideas used for developing the theoretical framework to analyse $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle$ for arbitrary $St$ are averaging decompositions and probability density function (PDF) methods. The only assumptions made in our derivation are that the system is statistically stationary and homogeneous, and that the particle dynamics is governed by (2.1). The basic steps in the derivation are as follows (see Tom & Bragg Reference Tom and Bragg2019 for more details). For the system under consideration, an ensemble average $\langle\,\cdot\,\rangle$ over all possible states of the system corresponds to averaging over all realizations of $\boldsymbol {u}$, as well as over all initial particle positions $\boldsymbol {x}^p(0)=\boldsymbol {x}_0$ and velocities $\boldsymbol {v}^p(0)=\boldsymbol {v}_0$. In view of this, we introduce the averaging decomposition $\langle\,\cdot\,\rangle =\langle \langle\,\cdot\,\rangle ^{\boldsymbol {x}_0,\boldsymbol {v}_0}_{\boldsymbol {u}}\rangle ^{\boldsymbol {u}}$ (Bragg, Swailes & Skartlien Reference Bragg, Swailes and Skartlien2012), where $\langle\,\cdot\,\rangle ^{\boldsymbol {x}_0,\boldsymbol {v}_0}_{\boldsymbol {u}}$ denotes an average over all initial particle positions $\boldsymbol {x}^p(0)=\boldsymbol {x}_0$ and velocities $\boldsymbol {v}^p(0)=\boldsymbol {v}_0$ for a given realization of the fluid velocity field $\boldsymbol {u}$, and $\langle\,\cdot\,\rangle ^{\boldsymbol {u}}$ denotes an average over all realizations of $\boldsymbol {u}$. Using this operator, we construct a generalized particle velocity field $\boldsymbol {\mathcal {V}}(\boldsymbol {x},t)\equiv \langle \boldsymbol {v}^p(t)\rangle ^{\boldsymbol {x}_0,\boldsymbol {v}_0}_{\boldsymbol {u}}$ that is well defined for all $St$, including the regime where caustics occur (Wilkinson & Mehlig Reference Wilkinson and Mehlig2005) because it is formed by averaging over trajectories satisfying $\boldsymbol {x}^p(t)=\boldsymbol {x}$ (i.e. it does not assume uniqueness of these trajectories, a scenario that only applies in the limit $St\to 0$ that was considered in Maxey Reference Maxey1987). With this velocity field we then construct the result

(2.4)\begin{equation} \langle{u}_z(\boldsymbol{x}^p(t),t)\rangle=\left\langle{u}_z(\boldsymbol{x},t)\exp\left( -\int_0^t \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathcal{V}}(\boldsymbol{\mathcal{X}}(s\mid \boldsymbol{x},t),s)\,{\rm d}s \right)\right\rangle^{\boldsymbol{u}},\end{equation}

where the notation $s\mid \boldsymbol {x},t$ denotes that the variable is measured at time $s$ along a trajectory satisfying $\boldsymbol {\mathcal {X}}(t)=\boldsymbol {x}$, and where $\partial _s\boldsymbol {\mathcal {X}}\equiv \boldsymbol {\mathcal {V}}(\boldsymbol {\mathcal {X}}(s),s)$.

The result was further developed to gain insight into the multiscale nature of the problem by introducing coarse-graining decompositions $u_z(\boldsymbol {x},t)=\tilde {u}_z(\boldsymbol {x},t)+u'_z(\boldsymbol {x},t)$ and $\boldsymbol {\mathcal {V}}=\tilde {\boldsymbol {\mathcal {V}}}+\boldsymbol {\mathcal {V}}'$, where $\tilde {u}_z$ and $\tilde {\boldsymbol {\mathcal {V}}}$ denote the fields $u_z$ and ${\boldsymbol {\mathcal {V}}}$ coarse grained on the length scale $\ell _c$, while $u'_z(\boldsymbol {x},t)\equiv u_z(\boldsymbol {x},t)-\tilde {u}_z(\boldsymbol {x},t)$ and $\boldsymbol {\mathcal {V}}'\equiv \boldsymbol {\mathcal {V}}-\tilde {\boldsymbol {\mathcal {V}}}$ are the ‘sub-grid’ fields. Next, we choose $\ell _c$ to be a function of $\textit {St}$, i.e. $\ell _c(\textit {St})$, which is defined such that the scale-dependent Stokes number $St_{\ell }\equiv \tau _p/\tau _{\ell }$ (where $\tau _{\ell }$ is the eddy-turnover time scale at scale $\ell$) is $\lll 1$ for $\ell >\ell _c(St)$. Physically, this means that $\ell _c(St)$ is defined to be the scale above which the particle behaves as if it were a fluid particle. Using this choice for the filtering length, we finally obtain

(2.5)\begin{equation} \langle{u}_z(\boldsymbol{x}^p(t),t)\rangle \approx \left\langle{u}'_z(\boldsymbol{x},t)\exp\left( -\int_0^t\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathcal{V}}'(\boldsymbol{\mathcal{X}}(s\mid\boldsymbol{x},t),s)\,{\rm d}s \right)\right\rangle^{\boldsymbol{u}}. \end{equation}

Since the right-hand side of this result only contains the sub-grid fields, it shows that the turbulent flow scales that contribute to the enhanced particle settling speeds are those with size $<\ell _c$, while scales $>\ell _c$ make a vanishingly small contribution. The physical mechanism embedded in (2.5) is a multiscale version of the original preferential sweeping mechanism described by Maxey (Reference Maxey1987) and Wang & Maxey (Reference Wang and Maxey1993). In particular, according to (2.5), $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <0$ occurs because the inertial particles are preferentially swept by eddies of size $\ell <\ell _c(St)$.

2.3. Applicability of theoretical analysis for the 2WC case

Although the result in (2.5) was derived for a 1WC system, its assumptions are such that it also applies in the 2WC regime. The only difference is one of interpretation, since in the 2WC case, all of the dynamical variables contained within (2.5) are implicitly affected by the force from the particle on the fluid. In view of this, we now consider how 2WC might affect the different quantities in (2.5) in order to understand the impact of 2WC on particle settling in turbulence.

The result in (2.5) shows that the settling modification due to turbulence associated with $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <0$ is determined by two things, namely the strength of the fluctuations of the sub-grid field ${u}'_z(\boldsymbol {x},t)$, and the intensity of the particle clustering (associated with the exponential term involving the compressibility $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\mathcal {V}}'$). It is to be noted, however, that as emphasized in Tom & Bragg (Reference Tom and Bragg2019), it is not merely the intensity of the clustering that matters but also whether it is correlated with ${u}'_z(\boldsymbol {x},t)$. For example, in the 1WC case, unless the clustering is correlated with ${u}'_z(\boldsymbol {x},t)$ then the right-hand side of (2.5) is zero because $\langle {u}'_z(\boldsymbol {x},t)\rangle ^{\boldsymbol {u}}=0$. For this reason, it is better to think of the exponential integral term in (2.5) as describing the process of preferential concentration rather than clustering, where preferential concentration describes how particles cluster preferentially in certain regions of the flow (see Bragg, Ireland & Collins (Reference Bragg, Ireland and Collins2015) for a discussion on the subtle differences between these phenomena).

In the 1WC regime, the particles do not modify ${u}'_z(\boldsymbol {x},t)$. In the 2WC regime, the particles will modify ${u}'_z(\boldsymbol {x},t)$, and previous results for non-settling particles have revealed a ‘pivoting effect’ wherein the inertial particles were observed to increase the energy content of turbulent flow scales of size less than some threshold, and decrease the energy content of flow scales larger than this (Elghobashi & Truesdell Reference Elghobashi and Truesdell1993; Sundaram & Collins Reference Sundaram and Collins1999; Bosse et al. Reference Bosse, Kleiser and Meiburg2006). The result in (2.5) shows that flow scales of size $\ell <\ell _c(St)$ contribute to the particle settling. The flows scales in the range $\ell <\ell _c(St)$ that dominate the settling will determine whether the pivoting mechanism tends to decrease or increase ${u}'_z(\boldsymbol {x},t)$. In the context of gravitational settling, Monchaux & Dejoan (Reference Monchaux and Dejoan2017) argued that the settling particles can drag the fluid in their vicinity with them, an effect referred to as the ‘fluid-drag mechanism.’ This fluid drag mechanism could lead to an increase in the negative values of ${u}'_z(\boldsymbol {x},t)$ responsible for generating $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <0$. In such a case, 2WC could enhance the particle settling speeds.

2WC can also, however, modify the preferential concentration of the particles associated with the contribution $\exp (-\int _0^t\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\mathcal {V}}'(\boldsymbol {\mathcal {X}}(s\mid \boldsymbol {x},t),s)\,{\rm d}s)$ in (2.5). It is more difficult to predict theoretically how 2WC might affect this contribution to (2.5), however, this is something that we can explore using DNS. It is important to note, however, that in a 2WC system, enhanced settling could occur even in the absence of preferential concentration. For example, even if $\exp ( -\int _0^t\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\mathcal {V}}'(\boldsymbol {\mathcal {X}}(s\mid \boldsymbol {x},t),s)\,{\rm d}s)= 1$ for all $\boldsymbol {x}^p(t)=\boldsymbol {x}$, then we would have

(2.6)\begin{equation} \langle{u}_z(\boldsymbol{x}^p(t),t)\rangle \approx \langle{u}'_z(\boldsymbol{y},t)\rangle^{\boldsymbol{u}}, \end{equation}

where $\boldsymbol {y}$ corresponds to the homogeneous distribution of points in space where the particles are located. Although in a 1WC system $\langle {u}'_z(\boldsymbol {y},t)\rangle ^{\boldsymbol {u}}=0$, in a 2WC system we could have $\langle {u}'_z(\boldsymbol {y},t)\rangle ^{\boldsymbol {u}}<0$ and hence $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <0$ due to the settling particles dragging the fluid down with them. In such a scenario the preferential sweeping mechanism would be playing no role. Measuring preferential concentration in the flow can therefore be used to indirectly determine whether $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle <0$ in a 2WC system is due to preferential sweeping or not.

It is also important to point out that, according to its definition, $\ell _c(St)$ will be the same for the 1WC and 2WC cases. To see this, we note that for point particles, the flow modification through 2WC is only possible due to particle inertia, and at any scale $\ell$, the role of particle inertia can be quantified using $St_\ell$. Since $\ell _c(St)$ is defined as the scale above which $St_\ell$ is so small that particle inertia can be ignored entirely, then this must also be the scale above which the effects of 2WC vanish, since at scales $\ell >\ell _c(St)$, the inertial particles behave as if they were fluid particles, which by definition do not modify the flow.

In our previous work (Tom & Bragg Reference Tom and Bragg2019), we studied the contributions of different flow scales to the settling velocity enhancement in a 1WC system. In this work, we aim to extend this analysis to the 2WC regime. First, we will look at the effect of 2WC on the overall settling velocity enhancement, i.e. the contribution of all the flow scales to $\langle {u}_z(\boldsymbol {x}^p(t),t)\rangle$. Second, we will study the effect of 2WC at different scales and also look at the relative contribution of different scales to the overall settling velocity enhancement. Third, we will analyse the effect of the two contributions embedded in (2.5) that affects settling velocity at different scales and investigate how these are modified by 2WC at different scales of the flow. Thus, we aim to elucidate how 2WC modifies the multiscale preferential sweeping mechanism.

3.1. Fluid solver

We solve the three-dimensional, incompressible Navier–Stokes equation in rotational form

(3.1)\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}

using the HiPPSTR code (Ireland et al. Reference Ireland, Vaithianathan, Sukheswalla, Ray and Collins2013), modified to account for the particle feedback on the flow $\boldsymbol {C}$. Here, $\boldsymbol {\omega }\equiv \boldsymbol {\nabla }\times \boldsymbol {u}$ is the vorticity, $P$ is the pressure, $\rho _f$ is the constant density of the fluid, $\nu$ is the kinematic viscosity, $\boldsymbol {f}$ is a large-scale forcing and $\boldsymbol {C}$ is the force exerted by the particles on the fluid. The spatial derivatives in (3.1) are discretized through a Fourier pseudo-spectral method. The nonlinearity introduces an aliasing error that is removed by using a combination of spherical truncation and phase shifting. The pressure gradient compensates the divergence of both the nonlinear convective term and the divergence of $\boldsymbol {C}$, thus yielding a solenoidal fluid velocity field. The Fourier modes are evolved in time by means of a second-order Runge–Kutta scheme with exponential integration for the viscous stress. A large-scale deterministic forcing scheme is applied at wavenumbers with magnitude $\kappa \leq \sqrt {2}$ that maintains a constant kinetic energy of the flow.

3.2. Particle solver

The particle position and velocity are evolved using (2.1), with the fluid velocity at the particle position computed through B-spline interpolation (van Hinsberg et al. Reference van Hinsberg, Thije Boonkkamp, Toschi and Clercx2012; Carbone & Iovieno Reference Carbone and Iovieno2018; Mirigaldi, Carbone & Perrone Reference Mirigaldi, Carbone and Perrone2021). This operation is understood as a backward non-uniform fast Fourier transform (NUFFT) with B-spline basis (Beylkin Reference Beylkin1995). The fluid velocity field is projected onto the B-spline basis in Fourier space, then brought back to physical space by means of an inverse fast Fourier transform (FFT). Finally, this projected velocity field is interpolated at the particle position. We employed B-spline polynomials of order seven, which uses eight points for interpolation along each direction and has $C^6$ continuity (Ireland, Bragg & Collins Reference Ireland, Bragg and Collins2016a). The particles position and velocity are then evolved in time by means of an exponential integrator (Hochbruck & Ostermann Reference Hochbruck and Ostermann2010), that guarantees stability and accuracy even for very low Stokes numbers. For a more complete description of the DNS solver see Ireland et al. (Reference Ireland, Vaithianathan, Sukheswalla, Ray and Collins2013). Note, however, that in the current version of the HiPPSTR code used to simulate the 2WC system, the same exponential integration scheme is used for time advancement of both the fluid and particle solvers (unlike the original HiPPSTR code which used exponential integrators only for the particles), to ensure consistency between the handling of the two phases.

3.3. Flow parameters

The particles simulated are monodisperse, having the same radius and density. In order to be able to compare our results with those of Monchaux & Dejoan (Reference Monchaux and Dejoan2017) we use the same density ratio that they employed, i.e. $\rho _p/\rho _f=5000$, the same Froude number, $\textit {Fr}=1$ and the same volume fraction, $\varPhi _v=1.5\times 10^{-5}$, corresponding to a mass fraction $\varPhi _m= 0.075$. We consider Stokes numbers $\textit {St}=0.3$, 0.5, 0.7, 1, 2 in order to explore the behaviour in the regimes of weak to moderate inertia. With these values of $\textit {Fr}$ and $\textit {St}$, the artificial periodicity effects discussed in Ireland et al. (Reference Ireland, Bragg and Collins2016b) do not arise, thus avoiding the need to use very large flow domains in the direction of gravity.

It is also desirable to consider the effects of the Taylor Reynolds number $Re_\lambda$ on the flow, since as shown in Tom & Bragg (Reference Tom and Bragg2019), the range of scales in the flow is of crucial importance in determining the particle settling behaviour via the multiscale preferential sweeping mechanism. We found that it was necessary to run the DNS for very long times, of the order of 100 eddy-turnover times in order for the average vertical velocities of the particles to statistically converge, and this is similar to the convergence times reported in Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) and Rosa et al. (Reference Rosa, Parishani, Ayala and Wang2016). This then places severe limitations on the $Re_\lambda$ that can be simulated (due to limited computational resources), and as a result, we considered the range $Re_\lambda \in [31,142]$. Although this range is too small to observe how the settling process might approach an asymptotic state with increasing $Re_\lambda$ for sufficiently large $Re_\lambda$ (as suggested in Tom & Bragg Reference Tom and Bragg2019), it will nevertheless be shown to be sufficient to reveal significant effects of $Re_\lambda$ on the settling process.

3.4. Two-way coupling

In the framework of the point-particle model, the force from the particles on the fluid is a superposition of Dirac delta functions centred at the particle position

(3.2)\begin{equation} \boldsymbol{C}(\boldsymbol{x},t) = \frac{L^3}{N_p}\varPhi_m \sum_p \frac{\boldsymbol{v}^p(t) - \boldsymbol{u}(\boldsymbol{x}^p,t)}{\tau_p} \delta(\boldsymbol{x} - \boldsymbol{x}^p(t)). \end{equation}

Note that, in this term, only the drag force on the particle is accounted for and not the gravitational body force experienced by the particle, as was also done in previous works (e.g. Bosse et al. Reference Bosse, Kleiser and Meiburg2006; Monchaux & Dejoan Reference Monchaux and Dejoan2017). The reason for this is that this term is supposed to be equal and opposite to the force that the fluid applies to the particles located at $\boldsymbol {x}$, and the latter is associated with the hydrodynamic stresses on the particle surfaces (represented by the drag force in the point-particle model being used) and not with the gravitational body force.

The Fourier transform of (3.2) needed by the pseudo-spectral solver for the fluid phase is computed by means of the forward NUFFT with B-spline basis (Beylkin Reference Beylkin1995; Carbone & Iovieno Reference Carbone and Iovieno2018; Carbone, Bragg & Iovieno Reference Carbone, Bragg and Iovieno2019). This is motivated by the fact that a FFT is not directly applicable to (3.2), while a direct discrete Fourier transform would be computationally too expensive. A detailed description of the NUFFT algorithm for particles in turbulence can be found in Carbone & Iovieno (Reference Carbone and Iovieno2018, Reference Carbone and Iovieno2020). Briefly, the algorithm can be summarized through the equation

(3.3)\begin{equation} \mathcal{F}[\boldsymbol{C}] = \frac{\mathcal{F}[\boldsymbol{C}*B_n]}{\mathcal{F}[B_n]}, \end{equation}

where $\mathcal {F}$ indicates the Fourier transform, $*$ denotes a convolution, and $B_n$ is the B-spline polynomial of order $n$.

The B-spline polynomials are piecewise polynomial functions with compact support, usually defined as the inverse Fourier transform of powers of the sinc function. The zeroth-degree B-spline consists of a box function,

(3.4)\begin{equation} B_0(\boldsymbol{x}) =\begin{cases} 1/\Delta x, & -\Delta x/2 < x_i<{+}\Delta x/2 \\ 0, & \textrm{otherwise} \end{cases} \end{equation}

and higher-order B-splines are obtained by taking recursive convolutions with $B_0$, namely $B_{n+1}=B_n*B_0$. These B-spline polynomials guarantee high-accuracy interpolation (Beylkin Reference Beylkin1995) and are efficient for algorithm parallelization purposes (Carbone & Iovieno Reference Carbone and Iovieno2018, Reference Carbone and Iovieno2020).

For an ideal computational grid with infinitely many grid points, (3.3) is just a tautology following from the convolution theorem. However, for a finite grid resolution, it is of great use in aiding the numerical implementation of a FFT of the field $\boldsymbol {C}$, consisting of a superposition of Dirac delta functions. In our implementation, we first compute the numerator of the right-hand side of (3.3), computing the convolution of $\boldsymbol {C}$ with $B_n$ in physical space, thus spreading the field over the $n+1$ adjacent grid points in each direction. This convolution procedure introduces a strong (artificial) non-locality, and the resulting field depends sensitively on the details of the smoothing function. The smoothed coupling field $\boldsymbol {C}*B_n$ can be accurately represented on the structured Cartesian grid, and it is then transformed to Fourier space by means of a standard FFT. Crucially, in the NUFFT algorithm the artificial spreading and non-locality introduced by the convolution $\mathcal {F}[\boldsymbol {C}*B_n]$ is subsequently removed in Fourier space by dividing by $\mathcal {F}[B_n]$, as shown in (3.3). The way in which this operation works can be easily understood as follows. A convolution (spreading) in physical space corresponds to multiplication in Fourier space, so that dividing by the Fourier transform of the B-spline amounts to taking an inverse convolution. Therefore, the smoothing is introduced only for numerical convenience as an intermediate step, and the outcome of the NUFFT calculation depends very weakly on the details of the kernel employed for convolution (Beylkin Reference Beylkin1995; Mirigaldi Reference Mirigaldi, Carbone and Perrone2021), the dependence being weaker as the number of grid points in the support of the kernel increases.

The NUFFT is a well-established algorithm in optics (Mirigaldi Reference Mirigaldi, Carbone and Perrone2021) and has been recently used to investigate the thermal coupling between a passive temperature field and inertial particles (Carbone et al. Reference Carbone, Bragg and Iovieno2019), and the results presented in this work have been crosschecked with the results from the code used in Carbone et al. (Reference Carbone, Bragg and Iovieno2019).

Most previous works used tri-linear extrapolation to distribute the contribution of the particle on the nearest grid points (Bosse et al. Reference Bosse, Kleiser and Meiburg2006). One of the main reasons for this choice is that a low-order polynomial reduces the artificial non-localities that arise due to spreading the effect of the particle momentum onto the surrounding grid points. In the NUFFT, the initial convolution in physical space is subsequently removed in Fourier space, so that the method does not suffer from the artificial non-locality issues, which is a major advantage of the method. As a consequence, this approach allows to use higher-order B-spline polynomials to obtain a more accurate representation of the Dirac delta functions in (3.2). Moreover, many recent works (e.g. Rosa et al. Reference Rosa, Kopeć, Ababaei and Pozorski2021) used high-order schemes for the interpolation of the fluid velocity at the particle position combined with low-order schemes for the computation of the coupling term. Since this introduces errors in the energy balance of the system (Sundaram & Collins Reference Sundaram and Collins1996), we employ the same B-spline polynomial basis of degree 7 (that uses 8 points in each direction) for interpolation of the velocity field and computation of the coupling term $\boldsymbol {C}$.

3.5. Numerical considerations

Within the context of the point-particle model, our simulations ignore the effects of the Basset history force, the Faxen curvature corrections, the pressure gradient and added mass contributions that appear in the full form of the equation derived in Maxey & Riley (Reference Maxey and Riley1983). However, the neglect of these terms is well justified for the particle sizes and densities being studied here (Daitche Reference Daitche2015). We also do not consider nonlinear drag corrections since our DNS are confined to the portion of the parameter range where a linear drag coefficient is a reasonable approximation.

With respect to the point-particle approach itself, it is known that for some parameter regimes the method has significant limitations (Eaton Reference Eaton2009; Balachandar & Eaton Reference Balachandar and Eaton2010; Brandt & Coletti Reference Brandt and Coletti2022). Only particle-resolved DNS allows for capturing all the flow details in the vicinity of the particles, and thereby overcoming the limitations of point-particle models (Vreman Reference Vreman2016). However, particle-resolved simulations are computationally extremely expensive for particles with $d_p\ll \eta$, and in view of this a point-particle approach is often used as an optimum choice (Brandt & Coletti Reference Brandt and Coletti2022). Moreover, in at least some parameter regimes, the point-particle approach should be valid provided that $d_p\ll \Delta x$ is satisfied and that the flow around the particle does in fact satisfy Stokes flow. However, we currently lack careful and extensive comparisons between point-particle DNS and experiments that would reveal exactly the parameter regimes for which the point-particle approach is valid (Petersen et al. Reference Petersen, Baker and Coletti2019; Brandt & Coletti Reference Brandt and Coletti2022). Hence, our DNS results which are based on the point-particle approach are to be placed within the context of the potential limitations of this modelling approach. But, as emphasized in § 2.1, this simple point-particle model is still widely used for describing particle motion in atmospheric flows (Grabowski & Wang Reference Grabowski and Wang2013), and even for particle motion governed by this simple model there remain many aspects that are not well understood. This again motivates the use of a point-particle model, before additional complexities such as finite particle size effects are introduced.

In the particle equation of motion (2.1), $\boldsymbol {u}(\boldsymbol {x}^p(t),t)$ is the undisturbed fluid velocity at the particle position, meaning the fluid velocity that does not include the disturbance of the particle under consideration (but can include the disturbance from all other particles in the flow). As a result, simply interpolating the fluid velocities from the surrounding grid points to evaluate $\boldsymbol {u}(\boldsymbol {x}^p(t),t)$ introduces an error because the interpolated velocities include the disturbance effect of the particle under consideration. This error can be significant for particles of size comparable to the grid spacing employed in the numerical simulation, i.e. when $d_p=O(\Delta x)$ (Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998). Moreover, this issue is particularly evident for particles settling in turbulence, and even in a still fluid (Horwitz & Mani Reference Horwitz and Mani2020). Correction schemes have been recently developed to approximately remove this error by providing a way to retrieve the undisturbed fluid velocity (Horwitz & Mani Reference Horwitz and Mani2018). However, in the regime $d_p\ll \Delta x$ the error will be small, and in our simulations $d_p/\Delta x \leq 0.05$. As a result, a correction scheme is not required, and we simply evaluate $\boldsymbol {u}(\boldsymbol {x}^p(t),t)$ in (2.1) by interpolating the fluid velocities at the surrounding grid point to the particle position, as described earlier.

The numerical resolution of DNS is usually quantified by $k_{max} \eta$, and the resolution of our DNS ($k_{max} \eta \approx 1.5$) is similar to or slightly better than the works of Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) and Monchaux & Dejoan (Reference Monchaux and Dejoan2017) against which compare our results. Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) used $k_{max} \eta \approx 1.29\unicode{x2013}1.55$ (see tables I and III) and Monchaux & Dejoan (Reference Monchaux and Dejoan2017) used $k_{max} \eta \approx 1.32$ (see table 1). While higher resolution might be required to capture intermittent events, our study focuses on first- and second-order moments and these are not very sensitive to intermittent events. Furthermore, if we were to increase the resolution of our DNS by reducing the grid spacing (for a given $St$ and $\varPhi _v$), then $d_p/\Delta x$ would enter the range where we would have to apply corrections to the fluid velocity in the particle equation of motion (Horwitz & Mani Reference Horwitz and Mani2018), and as mentioned above, this would add to the computational cost of the simulations. Hence, the current resolution is an optimal choice for tackling the questions we are trying to address in this work and for the parameter choices considered here, with reasonable computational resources.

Finally, following previous works (Bosse et al. Reference Bosse, Kleiser and Meiburg2006; Monchaux & Dejoan Reference Monchaux and Dejoan2017; Rosa et al. Reference Rosa, Kopeć, Ababaei and Pozorski2021), at each time step we set to zero the average fluid velocity in the direction of gravity. This is equivalent to applying a mean vertical pressure gradient (Bosse et al. Reference Bosse, Kleiser and Meiburg2006; Monchaux & Dejoan Reference Monchaux and Dejoan2017; Rosa et al. Reference Rosa, Kopeć, Ababaei and Pozorski2021) without which the total kinetic energy in the flow would grow without bound due to the finite settling velocity of the particles continuously injecting kinetic energy into the flow. While introduced in Maxey & Patel (Reference Maxey and Patel2001) and used for 2WC simulations by Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) to remedy the numerical issue that arises for an unbounded flow, its physical validity could be called into question. However, a possible justification can be given as follows. Suppose one is considering the settling of heavy particles through a fully developed wall-bounded turbulent flow with 2WC at distances from the wall where the flow is quasi-homogeneous. In this case, although the settling particles will continuously force the flow in the vertical direction as they fall, the impermeability of the wall prevents a mean flow from developing in this direction, and this effect is mediated to the upper, quasi-homogeneous region of the flow through the mean forces in the fluid. Hence in this situation, the settling particles would not generate a mean flow as they settle (although the mean fluid velocity evaluated at the particle positions can still be finite, and indeed will be if the flow modifies the particle settling velocity compared with the Stokes settling velocity). This potential justification is also hinted at in Maxey & Patel (Reference Maxey and Patel2001). In addition to this potential justification, we also use this aforementioned method of enforcing a zero mean fluid velocity at each time step in order to be able to meaningfully compare our results with those of the previous studies by Bosse et al. (Reference Bosse, Kleiser and Meiburg2006), Monchaux & Dejoan (Reference Monchaux and Dejoan2017) and Rosa et al. (Reference Rosa, Kopeć, Ababaei and Pozorski2021).

3.6. Simulation approach and post-processing methodology

We found that very long simulation times are necessary in order to obtain converged statistics of the average vertical fluid velocity at the inertial particle position and therefore the particle settling velocity, a requirement that was also mentioned in Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) and Rosa et al. (Reference Rosa, Parishani, Ayala and Wang2016). As such, our statistics are computed over the very long time window of 100 large-scale eddy-turnover times. We ran the fluid only simulation for 20 $\tau _L$ and obtained the flow statistics by averaging over the last 10 $\tau _L$ (see table 1). Next, we introduced particles into the flow and ran the 1WC case for 120 $\tau _L$ and computed the 1WC statistics over the last 100 $\tau _L$. Finally, we started the 2WC simulation from the instantaneous velocity and particle field of the completed 1WC run to obtain faster convergence of particle and fluid statistics and also because using the 1WC particle initial conditions in the 2WC regime do not affect the long-term statistics (Bosse et al. Reference Bosse, Kleiser and Meiburg2006). Similar to the 1WC case, the 2WC case was also ran for 120 $\tau _L$ and the 2WC statistics were computed over the last 100 $\tau _L$. Some of the previous works have used small averaging windows (for e.g. Good et al. Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014 and Ireland et al. Reference Ireland, Bragg and Collins2016b uses around 10 $\tau _L$) and as pointed out by Rosa et al. (Reference Rosa, Parishani, Ayala and Wang2016), this could be potential reason for different quantitative results between various works. The fluid statistics can change once the two-way coupling is activated and thus the 2WC fluid statistics are only known a posteriori. As is common in 2WC studies, the particle Stokes number $St$ in the 1WC and 2WC simulations are computed based on fluid statistics of the unladen flow presented in table 1.

Table 1. Flow parameters in DNS of unladen isotropic turbulence where all dimensional parameters are in arbitrary units. $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 {2k/3}$ is the root mean square of fluctuating fluid velocity, $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, $\Delta t$ is the DNS time step and $N_{proc}$ is the number of processors used for the simulations. The small-scale resolution, $k_{max}\eta$ and the total flow kinetic energy measured by $u'$ are approximately constant between the different simulations. All statistics are averaged over the last 10 $\tau _L$ of the total run time of 20 $\tau _L$ for the unladen simulations.

We investigate the contribution to the particle settling from different scales by computing relevant quantities filtered at various scales. We employ a sharp filter in Fourier space that filters out all the modes with wavenumber magnitude larger than the characteristic wavenumber $k_F$, which is associated with a filtering length $\ell _F=2{\rm \pi} /k_F$ (Eyink & Aluie Reference Eyink and Aluie2009) and obtain the sub-grid field through

(3.5)\begin{gather} {u}'_z(\boldsymbol{x},t)\equiv u_z(\boldsymbol{x},t)-\widetilde{u}_z(\boldsymbol{x},t)=\sum_{\|\boldsymbol{k}\| \geq k_F}\hat{u}_z(\boldsymbol{k},t)\,{\rm e}^{{\rm i}\boldsymbol{k\cdot x}}, \end{gather}

where $\widetilde{(\,\cdot\,)}$ denotes the coarse-grained field, while $(\,\cdot\,)'$ denotes the sub-grid field. We then interpolate $u'_z(\boldsymbol {x},t)$ to the positions of the inertial particles $\boldsymbol {x}^p(t)$ using an eight-point B-spline interpolation scheme to obtain ${u}'_z(\boldsymbol {x}^p(t),t)$. The values of ${u}'_z(\boldsymbol {x}^p(t),t)$ are then averaged over all the particles and over multiple times to obtain $\langle {{u}'_z}(\boldsymbol {x}^p(t),t)\rangle$. This process is then repeated for multiple $k_F$ in order to examine the contribution of the different flow scales. This is the same approach that was adopted in our previous paper (Tom & Bragg Reference Tom and Bragg2019) for analysing the multiscale preferential sweeping mechanism.