2.1. Background turbulent flow

Following the previous studies (Torres et al. Reference Torres, Parishani, Ayala, Rossi and Wang2013; Ayala et al. Reference Ayala, Parishani, Chen, Rosa and Wang2014), the turbulent flow is modelled using the Eulerian approach. The three-dimensional incompressible homogeneous isotropic turbulence is simulated in a cube of size $2{\rm \pi}$ using a pseudo-spectral method. The fast Fourier transform applied to the fluid velocity field implies periodicity at the domain boundaries in all three spatial directions. The Navier–Stokes equations are discretised on a three-dimensional uniform Cartesian mesh with $N$ equally spaced grid points. The governing equations are formulated in the rotational form

(2.1)\begin{gather} \frac{\partial \boldsymbol U}{\partial t}=\boldsymbol U \times \boldsymbol \omega - \boldsymbol \nabla\left(\frac{P}{\rho}+\frac{\boldsymbol U^2}{2}\right)+\nu \nabla^2\boldsymbol U+\boldsymbol f, \end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla \cdot U}=0, \end{gather}

where $\boldsymbol U(\boldsymbol x,t)$ and $\boldsymbol \omega (\boldsymbol x,t)$ are vectors of velocity and vorticity of the flow, respectively. Here $P(\boldsymbol x,t)$ is the pressure field; $\rho$ and $\nu$ denote the density and kinematic viscosity of air; $\,\boldsymbol f(\boldsymbol x,t)$ is the external body force acting on the fluid to maintain the turbulent flow. In the present study, we used a deterministic forcing scheme akin to the one developed by Sullivan, Mahalingam & Kerr (Reference Sullivan, Mahalingam and Kerr1994). In this method, the energy of the first two shells, $0.5<|\boldsymbol {k}|<1.5$ and $1.5<|\boldsymbol {k}|<2.5$, is preset to be constant and consistent with the Kolmogorov energy spectrum, i.e. $|\boldsymbol {k}|^{-5/3}$. For a desired Reynolds number, these values cannot be set arbitrarily since for a statistically stationary turbulence the average rate of energy input is equal to the average dissipation rate determined by the fluid viscosity. In our simulations the energy in these two shells is fixed to $E(1) = 0.55544$ and $E(2)=0.159843$, respectively. The added energy is distributed between 80 low wavenumber modes specified by $|\boldsymbol {k}| <2.5$. The effect of turbulence modulation by particles is not considered because in most simulations only diluted systems with low mass loadings of inertial droplets are modelled. In this study, the computational grids are limited to $64^3$ nodes since the main focus is on the effect of aerodynamic interaction, rather than the Reynolds number.

The basic parameters and time-averaged statistics of the flow are listed in table 1. These statistics were calculated based on the data collected after the flow reached the statistically stationary state. In table 1, $N$ is the number of grid points – resolution – in each spatial direction, ${\rm \Delta} t$ is the spectral time step chosen to evolve the flow, $\nu$ is the kinematic viscosity of the turbulent air given in spectral units, $\tau _{k}$ and $\eta$ are the Kolmogorov time and length scales, $T_e$ is the large-eddy turnover time, $L_f$ is the integral length scale, $\epsilon$ is the dissipation rate of turbulent kinetic energy, $u'$ is the root mean square of the fluctuations in velocity, $R_\lambda$ is the Taylor-microscale flow Reynolds number, ${CFL}$ is the Courant–Friedrichs–Lewy number and $k_{max}\eta$ is the commonly defined resolution parameter which must be larger than unity. Further, $\mathcal {S}$ and $\mathcal {F}$ are the skewness and flatness of longitudinal velocity gradient.

Table 1.

Input parameter settings and average flow statistics in spectral units.

2.2. Lagrangian particle tracking

Modelling of particle motion begins when the turbulent flow reaches a statistically stationary state. All particles are introduced into the domain at the same time instant. The initial locations are generated randomly enforcing uniform spatial distribution. Hereafter, the velocity and location of each particle are updated by integrating the equations of motion (Maxey & Riley Reference Maxey and Riley1983) as follows:

(2.3)\begin{gather} \frac{\mathrm d \boldsymbol{V}^{(k)}}{\mathrm d t}={-}\frac{\boldsymbol{V}^{(k)}-(\boldsymbol U^{(k)}+\boldsymbol{u}^{(k)})}{\tau_{p}^{(k)}}+\boldsymbol{g}, \end{gather}

(2.4)\begin{gather}\frac{\mathrm d \boldsymbol{Y}^{(k)}}{\mathrm d t}=\boldsymbol{V}^{(k)} , \end{gather}

in which $\boldsymbol {V}^{(k)}$ and $\boldsymbol {Y}^{(k)}$ are the velocity and location of $k$th particle, where $k=1,\ldots , N_{p}$. To simplify the notation, $\boldsymbol {U}^{(k)}\equiv \boldsymbol {U}(\boldsymbol {Y}^{(k)},t)$ is the background flow velocity $\boldsymbol {U}(\boldsymbol {x},t)$ interpolated at $\boldsymbol {Y}^{(k)}$. Next, $\boldsymbol {u}^{(k)}$ is the disturbance velocity resulting from the presence of all particles, except particle $k$, sensed at the location of $k$th particle. Symbol $\boldsymbol {g}$ is the gravitational acceleration and $\tau _{p}^{(k)}=2\rho _{p}(a^{(k)})^2/9\mu$ is the Stokes inertial relaxation time of particle $k$. The quantity is a measure indicating how quickly the particle responds to the changes in the flow field. The velocity of particles in (2.3) is updated using the implicit fourth-order Adams–Moulton (AM3) scheme for $\boldsymbol {V}^{(k)}$ in conjunction with the explicit fourth-order Adams–Bashforth (AB4) for the composite flow field: $\boldsymbol {U}^{(k)}+\boldsymbol {u}^{(k)}$. Similarly, AM3 is used for integrating equation (2.4) to advance the location of particles in time. The multistep schemes require keeping the record of velocities from previous time steps in memory arrays. The arrays are initialised with zeros. The initial condition for the velocity of each particle is the flow velocity at its location, to which its terminal velocity is added if the effect of gravity is considered. After each step of integration, periodicity is applied to the particles whose new locations fall outside the domain box; this simply means adding or subtracting the spectral length of the domain, $2{\rm \pi}$, to or from their updated locations.

In our numerical model, two simplifying assumptions were made concerning rotation and deformation of droplets. The rotation of particles in turbulent flow can be considered from two different perspectives. First, as the motion of particles driven by randomly oriented vortical tubes relative to a fixed point of reference. We consider this type of motion in our model. Second, as the rotation enforced by torques generated by the turbulent flow or by interactions with other droplets in close proximity. Consequences of these effects can be safely neglected for cloud droplets due to the fact that the droplet radii are much smaller than the Kolmogorov length scale (by a factor of ${\sim }50$) and their moments of inertia are relatively large (owing to the large particle-to-fluid density ratio). Thus, there is no efficient mechanism that forces the small droplets to rotate. The effect of rotation on the collision efficiency of two droplets in a quiescent fluid was studied in an earlier study by Rosa et al. (Reference Rosa, Parishani, Ayala, Wang and Grabowski2011a). They considered two scenarios, pure translational motion (rotation excluded) and translational–rotational motion. According to their results, the effect of rotation, in the absence of turbulence, is negligible for droplets of radii $30\ \mathrm {\mu }$m and larger.

The second simplifying assumption is that the droplets are treated as rigid bodies. For the considered ranges of Reynolds number and droplet sizes this approach accurately reflects the physics of the modelled processes and our simulation results do not lose their generality. A common measure used in the literature to quantify droplet deformation is the Weber number $(We)$. The number depends on the flow conditions and represents the relative importance of inertia and surface tension forces. For an isolated sphere the deformation is negligible if $We \ll 1$. In our simulations the maximal $We$ evaluated for $60\ \mathrm {\mu }$m settling droplets (in absence of turbulence) is approximately 0.035. The presence of other particles has no significant effect here because their radial relative velocity is very low.

2.3. Aerodynamic interactions among particles

In the original HDNS formulation (Ayala et al. Reference Ayala, Grabowski and Wang2007), the disturbance velocity felt by $k$th particle, $\boldsymbol {u}^{(k)}$, is computed by solving the linear system of equations

(2.5)\begin{equation} \boldsymbol{u}^{(k)}= \underbrace{\sum_{m=1}^{N_{p}}}_{m\neq{k}}{\boldsymbol{u}_{St}} (\boldsymbol{r}^{(k,m)};a^{(m)},{\boldsymbol{V}^{(m)}}- ({\boldsymbol{U}^{(m)}}+{\boldsymbol{u}^{(m)}})), \quad k=1,2,\ldots,N_{p}. \end{equation}

Here ${\boldsymbol {u}_{St}}$ is the velocity disturbance induced by an arbitrary droplet $m$ at the location of droplet $k$. In the limit of low Reynolds number, ${\boldsymbol {u}_{St}}$ can be approximated by the analytical formula representing the solution of the Stokes equation for the perturbation flow induced by a single sphere of radius $a^{(m)}$ moving at ${\boldsymbol {v}}^{(m)}$,

(2.6)\begin{align} \boldsymbol{u}_{{St}}(\boldsymbol{r}^{(k,m)};a^{(m)},{\boldsymbol{v}^{(m)}}) &=\left[\frac{3}{4}\frac{a^{(m)}}{{r}^{(k,m)}}-\frac{3}{4} \left(\frac{a^{(m)}}{{r}^{(k,m)}}\right)^{3}\right] \frac{\boldsymbol{r}^{(k,m)}}{({r}^{(k,m)})^2}(\boldsymbol{r}^{(k,m)} \boldsymbol{\cdot}\boldsymbol{v}^{(m)}) \nonumber\\ & \quad + \left[\frac{3}{4}\frac{a^{(m)}}{{r}^{(k,m)}}+\frac{1}{4} \left(\frac{a^{(m)}}{{r}^{(k,m)}}\right)^{3}\right] \boldsymbol{v}^{(m)},\quad k=1,2,\ldots,N_{p}, \end{align}

where $\boldsymbol {r}^{(k,m)} \equiv \boldsymbol {Y}^{(k)}-\boldsymbol {Y}^{(m)}$ is the vector connecting centres of particles $m$ to $k$, and its magnitude is $r^{(k,m)}=\|\boldsymbol {r}^{(k,m)}\|$. As it transpires from (2.5) and (2.6), the flow disturbance acting on each particle is a linear function of the composite velocity field, ${{\boldsymbol U}^{(m)}}+{\boldsymbol u}^{(m)}$, which implicitly depends on the disturbance velocities of all other particles. Consequently, by moving all unknown disturbance velocities, ${\boldsymbol u}^{(m)}$ u^(m), to the left-hand side of (2.5), the resulting system of equations in block matrix notation is

(2.7)\begin{equation} \underbrace{ {\begin{bmatrix} {{\boldsymbol I_3}} & \cdots & {\boldsymbol \alpha^{(1,m)}} & \cdots &{\boldsymbol \alpha^{(1,N_{\text p})}} \\ \vdots &\ddots &\vdots &\ddots &\vdots \\ {\boldsymbol \alpha^{(k,1)}} & \cdots & {\boldsymbol \alpha^{(k,m)}} & \cdots &{\boldsymbol \alpha^{(k,N_{\text p})}}\\\vdots &\ddots &\vdots &\ddots &\vdots \\ {\boldsymbol \alpha^{(N_{\text p},1)}} & \cdots & {\boldsymbol \alpha^{(N_{\text p},m)}} & \cdots &{{\boldsymbol I_3}} \end{bmatrix}}}_{{\boldsymbol A}} {\begin{pmatrix} \boldsymbol{u}^{(1)}\\\vdots\\\boldsymbol{u}^{(m)}\\\vdots\\\boldsymbol{u}^{(N_{\text p})}\end{pmatrix}} = {(\boldsymbol A - \boldsymbol I_{N_{\text p}})} {\begin{pmatrix}\boldsymbol{V}^{(1)}-\boldsymbol{U}^{(1)}\\\vdots\\\boldsymbol{V}^{(m)}-\boldsymbol{U}^{(m)}\\\vdots\\\boldsymbol{V}^{(N_{\text p})}-\boldsymbol{U}^{(N_{\text p})} \end{pmatrix}}, \end{equation}

in which ${\boldsymbol \alpha ^{(k,m)}}$ are $3\times 3$ symmetric matrices, and $\boldsymbol I_3$ is the identity matrix while $\boldsymbol I_{N_{\text p}}$ is a diagonal matrix with $\boldsymbol I_3$s on the main diagonal (see Torres et al. (Reference Torres, Parishani, Ayala, Rossi and Wang2013) for details). Accordingly, this system has $3N_{p}$ inseparable equations in scalar form. At each time step it is solved iteratively by means of a parallel solver developed based on the generalised minimal residual method, GMRES (Torres et al. Reference Torres, Parishani, Ayala, Rossi and Wang2013). In practice the matrix of coefficients is sparse since it is not necessary to consider interactions among all particles in the computational domain. This assumption is justified because the disturbance velocity induced by the particle decays as $1/r$ when $r\to \infty$, thus the effect of far-field aerodynamic interactions on the local interaction of two nearby particles may be safely neglected. Therefore, only perturbations due to droplets inside a spherical vicinity of each droplet are considered. The radius of this sphere is $H_{{tr}}\times a^{(k)}$, where $H_{{tr}}$ is a dimensionless truncation factor; hence, the larger the particle, the larger the neighbourhood scanned for perturbing droplets. Conducting a sensitivity analysis on the influence of truncation radius, Ayala et al. (Reference Ayala, Grabowski and Wang2007) demonstrated that considering perturbations of particles outside a truncation sphere with $H_{{tr}}=30$ does not make significant changes in collision statistics. Still, to more inclusively consider the neighbourhood around each droplet, $H_{{tr}}=50$ is used in this study.

Among many numerical methods, HDNS is the one that allows us to quite accurately represent the many-body interactions among widely separated particles. This desirable feature is important, especially for modelling systems with a large number of droplets. It is worth recalling that HDNS was used in several previous studies for modelling systems that are relevant to typical atmospheric clouds (Wang et al. Reference Wang, Ayala, Kasprzak and Grabowski2005b; Ayala et al. Reference Ayala, Grabowski and Wang2007; Wang, Ayala & Grabowski Reference Wang, Ayala and Grabowski2007; Wang et al. Reference Wang, Ayala, Rosa and Grabowski2008, Reference Wang, Rosa, Gao, He and Jin2009). Nevertheless, it should be pointed out that HDNS does not accurately represent the short-range lubrication forces. The method fails when the gap between particles is significantly smaller than the radii of the particles. This shortcoming is an indirect effect of using the ISM (Wang et al. Reference Wang, Ayala and Grabowski2005a) to construct the numerical algorithm for computing particle drag. The discrepancy can be quantitatively evaluated in a simplified case of two spheres in translational motion. In figure 1, the normalised drag forces, $F/6{\rm \pi} \mu a V$, acting on two droplets approaching or receding from each other along their centreline, at velocity $V$, are shown at various dimensionless gap distances: $s-2$, where $s\equiv 2r/(a_1+a_2)$ is the dimensionless separation distance, i.e. the distance between the centres of two arbitrary interacting particles. Defining the collision radius $R\equiv a_1+a_2=2a$, then $s=2r/R$.

Figure 1.

Normalised magnitude of the drag force acting on two same-size spherical particles, whether approaching or receding, along their line of centres as a function of non-dimensional gap between their surfaces.

In figure 1, the exact analytical prediction of drag forces (Jeffrey & Onishi Reference Jeffrey and Onishi1984) indicates that the forces asymptotically grow as the gap distance decreases, i.e. $s-2\to 0$. This, however, is valid only until the separation distance between the droplets is of the order of the mean free path of air molecules beyond which the molecular nature of the fluid becomes important (Sundararajakumar & Koch Reference Sundararajakumar and Koch1996). Comparison of these analytically predicted forces with those obtained from the ISM (Wang et al. Reference Wang, Ayala and Grabowski2005a) shows that the ISM does not yield an accurate representation of lubrication forces between two droplets when their gap distance is comparable to their average radii, i.e. $s-2<1$. To overcome this problem, a new numerical approach is proposed here in which the standard HDNS is only used to represent the interactions between widely separated droplets while the drag of nearly touching droplets is computed using the exact analytical solution derived by Jeffrey & Onishi (Reference Jeffrey and Onishi1984). This strategy is somewhat similar to the one proposed by Durlofsky et al. (Reference Durlofsky, Brady and Bossis1987, (2.18) therein).

In the new approach the disturbance velocity acting on $k$th droplet, $\boldsymbol {u}^{(k)}$, is a superposition of two components: $\boldsymbol {u}^{(k)}_{{HDNS}} + \boldsymbol {u}^{(k)}_{{lub}}$. The first component, $\boldsymbol {u}^{(k)}_{{HDNS}}$, is evaluated using the HDNS algorithm with a modification that excludes the pairs separated at some preset distance $\delta$ (we call it the matching point), namely $s\le \delta$. The second component is only computed for all nearby pairs, $s\le \delta$, using the exact analytical formula derived by Jeffrey & Onishi (Reference Jeffrey and Onishi1984).

In most of simulations we used $\delta =3$. This means that the numerical method switches from HDNS to the exact solution when the gap size between droplet surfaces is smaller than the average size of the radii of the droplets. Sensitivity study of the modelled collisional statistics for different $\delta$ was also carried out and the results are presented and discussed in § 3. Contrary to other numerical models, the new approach allows us to simultaneously include two effects, namely many-body interaction and lubrication forces. The only assumption made is that the binary interactions between nearby particles are not affected by the presence of other particles. This simplification has a rather negligible impact on the modelled systems, since the magnitude of short-range lubrication forces significantly exceeds the magnitude of interaction forces between particles separated by more than $s = 3$. While the lubrication forces can be influenced by many-body interactions, such a scenario is unlikely in dilute systems such as those considered in this study.

The forces acting on two nearby particles in the Stokes approximation are expressed in terms of the resistance matrix (Jeffrey & Onishi Reference Jeffrey and Onishi1984)

(2.8)\begin{equation} \frac{1}{3{\rm \pi} \mu(a_1+a_2)} {\begin{pmatrix} {\boldsymbol{F}^{(1)}_{{lub}}} \\ {\boldsymbol{F}^{(2)}_{{lub}}} \end{pmatrix}} = {\begin{pmatrix} {\boldsymbol{u}^{(1)}_{{lub}}} \\ {\boldsymbol{u}^{(2)}_{{lub}}} \end{pmatrix}} = {\begin{bmatrix} \boldsymbol{\mathsf{A}}_{11} & \boldsymbol{\mathsf{A}}_{12} \\ \boldsymbol{\mathsf{A}}_{21} & \boldsymbol{\mathsf{A}}_{22} \end{bmatrix}} {\begin{pmatrix} \boldsymbol{V}^{(1)}-\boldsymbol{U}^{(1)} \\ \boldsymbol{V}^{(2)}-\boldsymbol{U}^{(2)} \end{pmatrix}}, \end{equation}

where $\mu$ is the dynamic viscosity of the fluid. Here, each component of the resistance matrix, $\boldsymbol{\mathsf{A}}_{\alpha \beta }$, is a non-dimensional second-rank tensor consisting of resistance functions $X^{\boldsymbol{\mathsf{A}}}_{\alpha \beta }(s,\lambda )$ and $Y^{\boldsymbol{\mathsf{A}}}_{\alpha \beta }(s,\lambda )$ that depend on the dimensionless separation distance and particles radii ratio $\lambda ={a_2}/{a_1}$. In the present study only monodisperse systems are modelled, hence $s=r/a$ and $\lambda =1$. Moreover, the effect of particle rotation is not considered, therefore torques and other rotational terms are already excluded from the resistance matrix in (2.8). Defining particle velocity relative to the background flow $\boldsymbol {v}^{(\beta )}\equiv \boldsymbol {V}^{(\beta )}-\boldsymbol {U}^{(\beta )}$, each tensor multiplication in (2.8) takes the following form:

(2.9)\begin{equation} \boldsymbol{\mathsf{A}}_{\alpha\beta} \boldsymbol{v}^{(\beta)} = X^{\boldsymbol{\mathsf{A}}}_{\alpha\beta}\frac{\boldsymbol{r}}{r^2} (\boldsymbol{r \cdot v}^{(\beta)}) + Y^{\boldsymbol{\mathsf{A}}}_{\alpha\beta} \left(\boldsymbol{v}^{(\beta)}-\frac{\boldsymbol{r}}{r^2} (\boldsymbol{r \cdot v}^{(\beta)})\right),\end{equation}

where $\boldsymbol {r}\equiv \boldsymbol {Y}^{(2)}-\boldsymbol {Y}^{(1)}$. The first term is related to the drag force due to the squeezing motion of spheres along the line connecting their centres. The drag enforced by the shearing motion of spheres perpendicular to the line connecting their centres is handled by the second term. The resistance coefficients $X^{\boldsymbol{\mathsf{A}}}_{\alpha \beta }(s,1)$ and $Y^{\boldsymbol{\mathsf{A}}}_{\alpha \beta }(s,1)$ can be computed analytically using the procedure of twin multipole expansion developed by Jeffrey & Onishi (Reference Jeffrey and Onishi1984). The method is rather complex because the coefficients are given in terms of infinite series summations. Moreover, each term of the series can be determined only by solving complex recurrence relations. Further, the convergence rate of the series varies and depends on the specific configuration of the system such as the relative locations of the spheres. To obtain a satisfactory precision, more terms are needed when the particles are in close proximity. Due to high complexity, the method cannot be used efficiently for modelling systems with many pairs of particles. Therefore, in our simulations the coefficients were pretabulated and then interpolated ‘on the fly’ using the Bulirsch–Stoer rational interpolation algorithm: ratint.f77 (Press et al. Reference Press, Teukolsky, Flannery and Vetterling1992). This interpolation scheme has been adopted for a specific range of $s$ since, as is shown in figure 1, the drag force resulting from the Stokes solution tends to infinity when $s-2\to 0$. In real applications this drag force is lower since the continuum assumption of the fluid breaks down when the gap is of the order of the mean free path of air molecules (Sundararajakumar & Koch Reference Sundararajakumar and Koch1996). Besides, there would be potential numerical instabilities as a result of extremely large drag forces if the gap between droplets were extremely small. To address these problems, we made use of ‘the finite-gap model’ proposed by Hocking & Jonas (Reference Hocking and Jonas1970). Basically, the algorithm for collision detection has been modified in a way that collision between two particles is assumed even if the minimum gap between their surfaces is greater than zero but smaller than $\epsilon _R$. In other words, the collision radius is slightly enlarged to $R=(2+\epsilon _R)a$, where $\epsilon _R=10^{-3}$. To assess sensitivity of our approach we performed a set of test simulations with droplets of radii $40\ \mathrm {\mu }$m using different values of $\epsilon _R$ in the range $(0.25\text {--}2.5) \times 10^{-3}$. Based on the obtained results, not presented here, we conclude that this approach guarantees the numerical stability of the code while the differences in collision statistics are within statistical uncertainty. It should be added that after collision, one of the particles is relocated and the problem with the large drag force acting on the particles is in all situations eliminated.

3.1. Sensitivity of droplet collision statistics to time step size

Including the lubrication forces into the model has a direct impact on the numerical properties, e.g. stability of integration, of the ordinary differential equations for tracking particles, i.e. (2.3) and (2.4). Due to strong nonlinearity of these short-range aerodynamical forces, the resulting set of governing equations is stiff. As such, to guarantee the numerical stability, the integration scheme has to be properly adjusted. The time step size is of special importance in this regard and hence must be carefully optimised. If ${\rm \Delta} t$ is too large, the algorithm may be unstable. Moreover, some collision statistics, for example the collision kernels, may be overestimated. This is because for a large ${\rm \Delta} t$ the lubrication forces are sampled with an excessively coarse spatial resolution. On the other hand, if ${\rm \Delta} t$ is very small the numerical cost of computations will be needlessly large. Therefore, to examine this problem, a number of simulations have been performed to find an optimal ${\rm \Delta} t$. The largest time step, ${\rm \Delta} t_{max}=4.24\times 10^{-4}$ s, is the same as the one used to evolve the turbulent flow without particles, i.e. the spectral ${\rm \Delta} t$ in table 1. Then, in subsequent simulations, the time step was consistently reduced by half up to ${\rm \Delta} t_{max}/32$. It should be added that the total simulation time, roughly $100T_e$, was the same in all modelled cases.

The numerical experiments were divided into two groups. The first set of simulations was performed for droplets of relatively low inertia. In these simulations the Stokes number, defined as $St\equiv \tau _{p}/ \tau _{k}$, did not exceed $St$ = 0.36. For such systems the effects of gravitational settling can be safely neglected. On the other hand, by referring to some former studies, e.g. Hocking & Jonas (Reference Hocking and Jonas1970) and Rosa et al. (Reference Rosa, Wang, Maxey and Grabowski2011c), it is expected that the effect of lubrication forces should be important. The droplet inertia was set by changing the value of energy dissipation rate. The droplet radii were fixed to $a_{p} = 40\ \mathrm {\mu }$m (according to the notation above, $a_{p}=a$). In the second group of simulations we considered several different sizes of particles at fixed energy dissipation rate equal to $\epsilon = 400\ \mathrm {cm}^2\ \mathrm {s}^{-3}$. This setting is representative for moderate to strong convection in a cloud. For such systems gravity is important, hence the simulations were performed for both settling droplets and without gravity.

The effect of the time step size is assessed by comparing two-point collision statistics of almost touching water droplets. Figure 2(a) shows the RRV normalised by the Kolmogorov velocity scale, computed in simulations with different ${\rm \Delta} t$, which are normalised by the particle inertial response time. Corresponding values of the RDF are presented in figure 2(b). In all these simulations the lubrication forces were considered. The main conclusion resulting from these simulations is that the statistics do not change much for ${\rm \Delta} t < 0.04 \tau _{p}$. This seems to be a trade-off between accuracy and numerical complexity.

Figure 2.

Two-point collision statistics at contact computed using different time steps: (a) normalised RRV and (b) RDF.

The increase of the RRV at larger ${\rm \Delta} t$ is a consequence of a not sufficient representation of the short-range aerodynamic interaction. It should be recalled that the lubrication forces not only prevent collisions of approaching particles, but also restrain the pairs from receding from each other. This reversibility is an intrinsic feature of the solutions to the Stokes flow. As a result, the relative motion of the interacting droplets is more strongly correlated. In other words, the lubrication forces act similarly to a hydraulic shock absorber and hence their effect must be thought of as that of a viscous damper, rather than a spring, in a simplified mass–spring–damper model. This can be readily confirmed by taking a look at (2.8) in which the drag coefficients are multiplied by the velocities. This is a general feature of aerodynamic interaction forces, but in the case of lubrication forces this barrier is harder for particles to break through or escape from.

The core results analysed in the present study were computed assuming the energy dissipation rate $\epsilon = 400\ \mathrm {cm}^2\ \mathrm {s}^{-3}$. Therefore, the sensitivity of these simulations to the time step is analysed more deeply. First, in figure 3(a) we compare kinematic collision statistics computed in simulations without gravity at different separation ranges. As expected, for larger separation distances the RRV is not sensitive to the time step size. This can be explained by the fact that the lubrication forces are effective only at short ranges (see figure 1). In turn, at near-contact separation distances, $1<r/R<1.5$, some minor changes in the RRV can be noticed. Accordingly, this region is enlarged and displayed in figure 3(b). The data show slightly lower relative velocities for smaller time steps. This effect stems from the improvement in capturing short-range forces which tend to decrease the magnitude – whether positive or negative – of the relative velocities between particles. The smaller the time step size, the more accurate the representation of the lubrication forces. Interestingly, for ${\rm \Delta} t < 5.07 \times 10^{-3} \tau _{p}$ the differences are relatively low. In the next step, we address sensitivity of the particle clustering to ${\rm \Delta} t$. The RDFs computed in the same numerical simulations are shown in figure 4. As before, figure 4(b) presents the RDF for near-contact separation distances. Similar to the RRV, in general the RDF does not change much for different time step sizes. The only difference is observed at separation distances comparable to the particles radii. The RDF increases as the time step decreases, which is a consequence of including the lubrication forces. As already discussed, the relative velocity decreases at shorter ${\rm \Delta} t$, which in turn enhances the tendency of pairs of droplets to stay together, thereby augmenting their average concentration, especially at short distances.

Figure 3.

Radial relative velocity computed using different time step sizes for the normalised separation distance in the range (a) $1 < r/R < 10$ and (b) $1 < r/R < 1.5$ with the set of droplets $a_{p} = 40\ \mathrm {\mu }$m and $N_{p} = 50\ 000$. The black rectangle in panel (a) marks the region that is enlarged and shown in panel (b).

Figure 4.

Radial distribution function computed using different time step sizes for the normalised separation distance in the range (a) $1 < r/R < 10$ and (b) $1 < r/R < 2$ with the set of droplets $a_{p} = 40\ \mathrm {\mu }$m and $N_{p} = 50\ 000$. The black rectangle in panel (a) marks the region that is enlarged and shown in panel (b).

The above analyses were performed for systems with a fixed liquid water content and droplets of the same radii, $a_{p} = 40\ \mathrm {\mu }$m. Here we extend the analysis and consider systems with particles of different radii with various or identical mass loadings. Since the lubrication forces are important mainly at short separation distances, the current analysis is restricted to the kinematic collision statistics of the droplets at contact. Figure 5(a) shows the RRV, normalised by the Kolmogorov velocity scale, of the sets of droplets with five different radii simulated with different ${\rm \Delta} t$. It is important to note that the simulations in figure 5(a,c) were performed with the same number of droplets, so that the mass loading in individual series is different. Even though the number of droplets is the same, the difference in the size of droplets changes the mass loading. Consequently, simulations with identical mass loading, $\phi = 1.22\times 10^{-2}$, were performed as well, and the results are presented in figure 5(b,d), with ${\rm \Delta} t$ being normalised by $\tau _{k}$, directly corresponding to the set of six time steps chosen and listed on the top horizontal axes. The results shown in figure 5(a,b) allow us to conclude that RRVs increase with ${\rm \Delta} t$. As expected, the larger increase is observed for particles of lower inertia, which is consistent with figure 2. It is also noteworthy to mention that the sensitivity of these statistics to ${\rm \Delta} t$ depends on both the droplet inertia and gravity. The increase in RRV is not strictly related to the aerodynamic interactions, since in both series of simulations the same mass loading was considered. The enhancement of RRV at larger ${\rm \Delta} t$ should be understood as a more effective decorrelation of droplet motion caused by a coarser sampling of the background fluid velocity. In simulations without gravity the RRV is dominated by so-called path history or non-local effects. Gravity alters this mechanism and reduces the correlation between fluid and particles. This is most likely the reason why for droplets of radii $30\text {--}40\ \mathrm {\mu }$m, a slightly larger sensitivity is observed in simulations with gravity. For heavier droplets, i.e. $50\ \mathrm {\mu }$m, the dependency of the kinematic statistics, especially the RRV, to ${\rm \Delta} t$ is larger if gravity is not considered. This is because gravity reduces the droplets’ velocity fluctuations (relative to the fluid). Thus, the RRV is less dependent on ${\rm \Delta} t$. In the limiting case of particles with infinite inertia, the turbulence effects on the RRV are negligible, and the constraints on the time step become insignificant.

Figure 5.

The effect of time step size on the RRV at contact for (a) the same number of droplets, $N_{p} = 50\ 000$ and (b) the same mass loading, $\phi = 1.22\times 10^{-2}$, as well as RDF at contact for (c) the same number of droplets, $N_{p} = 50\ 000$ and (d) the same mass loading, $\phi = 1.22\times 10^{-2}$, considering various sets of same-size droplets with different radii; simulations with gravity are labelled with a ($g$).

Corresponding to figure 5(a,b), RDFs of systems with the same number of droplets and the same mass loading are shown in figures 5(c) and 5(d), respectively. An opposite decreasing trend is observed for the at-contact RDFs. Again, the observed dependencies can be interpreted as consequences of including the lubrication forces. Large time steps do not allow adequate capturing of the short-range aerodynamic forces, so that the effect of particle hampering is smaller and consequently the RRV is overestimated. Larger RRV results in higher probability of particle collision. This in turn has an indirect effect on reducing the RDF, since after every collision one of the droplets is moved to a new location so that, in an average sense, the droplets have a more uniform distribution, leading to smaller RDFs. Since the effect of the time step size is different for the RRV and RDF, the question arises about the net effect of ${\rm \Delta} t$ on the collision kernel – the parameter directly proportional to the collision rate. The values of the dynamic collision kernel with respect to different time step sizes are presented in figure 6. Contrary to the kinematic statistics, the dynamic collision kernel does not show considerable changes with the time step size. This may be explained as a result of the mutual compensation of these two opposite effects. Here, it is worth referring to the kinematic definition of the collision kernel, namely $\varGamma ^{K}$ (see (3.8)), which is proportional to the product of the RRV and RDF. This experiment shows that even though the dynamical characteristics are similar for different choices of ${\rm \Delta} t$, the kinematic statistics may be different. In figure 6(b) we observe an increase in numerical uncertainty with the droplet radii. This larger uncertainty especially in the collision kernels of $50\ \mathrm {\mu }$m droplets is due to fewer number of detected collisions, which itself is a direct result of lower number of drops used in these simulations to maintain the same mass loading.

Figure 6.

The effect of time step size on the dynamic collision kernel for (a) the same number of droplets, $N_{p} = 50\ 000$ and (b) the same mass loading, $\phi = 1.22\times 10^{-2}$, considering various sets of same-size droplets with different radii.

The above considerations lead to the final conclusion that for time steps smaller than ${\rm \Delta} t=1.06\times 10^{-4}$ s, the changes in the kinematic and dynamic statistics are rather insignificant. Thus, the rest of simulations in the present study are conducted with this time step size.

3.2. Setting an optimal location for the matching point $\delta$

As already pointed out, standard HDNS is not capable of accurately representing the aerodynamic forces acting on spherical particles if the separation distance between them is smaller than or comparable to their mean radius (see figure 1). To overcome this problem, we proposed a new numerical model in which computation of the drag forces switches from HDNS to the exact formulation of Jeffrey & Onishi (Reference Jeffrey and Onishi1984) at a heuristically chosen matching point. To maximise fidelity of our model, the location of the matching point for these two approaches should be properly optimised.

It should be recalled that the exact formulation of Jeffrey & Onishi (Reference Jeffrey and Onishi1984) has been derived for two rigid spheres immersed in a viscous fluid and is limited to the Stokes regime. Therefore, this approach is suitable for modelling a binary system of interacting particles. Still, this binary treatment can be accurate for simulating multiphase flows if the volume fraction of the dispersed phase is relatively small. In such systems there is a low probability of three or more droplets simultaneously being at a separation distance comparable to their radii. On the other hand, this exact treatment is binary and thus unable to capture the effects of many-body interactions among particles, while these effects can be faithfully handled under the HDNS approach. Considering the upsides and downsides of these two approaches we have to find an optimal location for the matching point of them.

The optimisation procedure consists of comparing collision statistics, both kinematic and dynamic, from a series of simulations performed with different choices of $\delta$. First, we analyse the RRV and RDF at different separation ranges, namely, $1<r/R<10$. This is a necessary step, since the choice of $\delta$ may have some impact on the statistics. Figure 7 shows the normalised RRV and RDF of two monodisperse systems containing droplets of different radii/inertia but the same mass loading, $\phi = 1.22\times 10^{-2}$. The range of considered $\delta$ extends from 2 to 15. In figure 7(a,b) a noticeable reduction in the RRV is seen as soon as the binary formulation is applied within separation distances $s<\delta =3$. This is more obvious in the near-contact region which is marked and enlarged for each figure. A slight but systematic growth in the relative velocity is observed in simulations with larger $\delta$. Surprisingly, for a larger $\delta$, or equivalently a larger range of exact binary treatment, the RRV approaches the solution computed with the standard HDNS. This result may be counterintuitive and has to be interpreted along with the RDF. The corresponding values of the RDF increase with $\delta$, which can be explained as an effect of stronger droplet hampering due to the lubrication effect. This in turn results in a more intense clustering at droplet size scales. Stronger clustering and shorter separation distances enhance the effect of relative motion due to many-body interactions. It is important to note that $g_r(r=R)$ computed using the coupled model when $\delta = 3$ is significantly larger than the corresponding value from the standard HDNS, i.e. $\delta = 2$. The results suggest that it is justified to use the exact binary treatment for $s<3$. In figure 8, we present at-contact RRV and RDF computed for three systems with different droplet radii but the same mass loading. Again, the binary treatment is applied within different ranges of separation distances such that $\delta =2$ indicates the standard HDNS approach. The RRVs at contact show a decline when the lubrication forces are included for separation distances smaller than $\delta =3$. This confirms that the lubrication interaction prevents the particles from both approaching and receding from each other. The decrease in $\langle |w_r(r=R)| \rangle$ depends on the droplet inertia, being more noticeable for the larger droplets. This effect can be attributed to a stronger decorrelation of the fluid and particle velocities. Looking at figure 8(b), at-contact RDFs indicate stronger clustering when the lubrication forces are considered because the reduction in RRV allows the droplets to stay longer at closer separations.

Figure 7.

Variations in the RRV (a,b) as well as RDF (c,d) for the sets of droplets (a,c) $a_{p} = 40\ \mathrm {\mu }\mathrm {m}, N_{p} = 50\ 000$, and (b,d) $a_{p} = 30\ \mathrm {\mu }\mathrm {m}, N_{p} = 118\ 500$ when the lubrication forces are considered within different ranges of interaction $\delta$.

Figure 8.

Variations in the (a) at-contact RRV and (b) at-contact RDF for three sets of droplets with the same mass loading, $\phi = 1.22\times 10^{-2}$, when the lubrication forces are considered within different ranges of interaction.

For completeness of the analysis, we present the dynamic collision kernels of these three systems modelled with different $\delta$ in figure 9. As expected, the dynamic collision kernel is reduced when the lubrication forces are considered. In other words, ignoring the effect of lubrication leads to an overestimation of the collision rate. This reduction is mainly due to smaller values of the RRV between droplets. It is interesting to note that for $\delta =3$ the relative reduction in $\varGamma ^{D}$ is smaller than that in the RRV. This mitigation is mainly due to stronger clustering, demonstrated by the larger values of the RDF at contact presented in figure 8(b).

Figure 9.

Variations in the dynamic collision kernel for three sets of droplets with the same mass loading, $\phi = 1.22\times 10^{-2}$, when the lubrication forces are considered within different ranges of interaction.

4.1. Comparison of the new lubrication-included HDNS with standard HDNS and simulations without aerodynamic interactions

This section gives a closer insight into the differences between three numerical approaches. We compare collision statistics computed using (i) the new model with full representation of the short-range interaction forces, (ii) standard HDNS approach and (iii) simulations under one way coupling in which the effect of aerodynamic interactions is neglected. Figure 12 shows the normalised RRVs, in panel (a), and RDFs, in panel (b), at contact as a function of the droplet radius. The three series of simulations were performed using different numerical approaches and the same settings. Both the number of droplets ($N_{p} = 50\ 000$) and the energy dissipation rate ($\epsilon = 50$ cm$^2$ s$^{-3}$) were fixed in all simulations. The gravitational settling of the droplets was not considered. The obtained results are consistent with previous observations. The RRV is significantly augmented in the range of smaller droplets if the lubricative forces are included. Although the effect of lubrication is clearly visible in the RRV, in the average sense, this sharp increase does not have much influence on the collision statistics (see figure 11) since the number of detected pairs at short separation distances is very low. This can be confirmed by looking at figure 12(b) which shows RDFs at contact for the same range of sizes. For larger droplets the RRV computed with the new approach falls in between two other models. The increase over the simulations without aerodynamic interactions is a combined effect of short-range lubrication forces and many-body interactions. The most pronounced enhancement of the RRV was obtained in hybrid DNS. This result is rather expected since in HDNS the motion of the particles is not hampered by the lubricating forces.

Figure 12.

Changes in the at-contact (a) RRVs and (b) RDFs when the effect of lubrication forces is taken into consideration compared with the standard HDNS without lubrication effects as well as the case without aerodynamic interaction for $N_{p} = 50\ 000$ droplets in a flow at a low dissipation rate $\epsilon = 50$ cm$^2$ s$^{-3}$.

Figure 12(b) shows a considerable reduction in at-contact RDFs when aerodynamic interactions among droplets are considered. Whereas, in other studies, e.g. Brunk, Koch & Lion (Reference Brunk, Koch and Lion1997) and Yavuz et al. (Reference Yavuz, Kunnen, van Heijst and Clercx2018), an opposite trend was reported indicating that aerodynamic interactions enhance RDF. We believe that this discrepancy is due to some numerical effects resulting from assumptions made here. An influential part of the code used in the present study is the droplet relocation module. The module determines a new random location for one of the droplets of a colliding pair. This treatment mimics the well known physical process, in which two colliding droplets form a larger drop. This large drop does not belong to the considered size class. At the same time a new droplet formed by diffusion or collision-coalescence is introduced into the domain. By turning this module on, even if the aerodynamic interactions are not considered, the mechanism of clustering is ‘artificially’ altered, and the droplet distribution tends to become more uniform; consequently, the RDF values are smaller. It should also be emphasised that for the largest droplets $(60\ \mathrm {\mu }\mathrm {m})$ the RDF at contact is slightly larger if lubrication is added. This demonstrates that the lubrication interaction of such droplets prevents their collisions more efficiently.

Figure 13 compares the RRVs of monodisperse systems at a higher energy dissipation rate equal to $\epsilon = 400$ cm$^2$ s$^{-3}$. Again, the simulations were performed using different numerical approaches but this time we separately considered the cases with and without gravity. Since the effect of the droplet mass loading may be important, the analysis includes both the simulations at a fixed number of droplets, figure 13(a), and at a fixed mass loading, figure 13(b). Obtained results show that, compared with the standard HDNS, the RRV decreases when the additional effect of lubrication is taken into consideration. This is due to the larger drag forces that are inherently not accounted for by the HDNS. Moreover, there is a general increase in the relative velocity with the size of droplets, which arises from the higher inertia of larger droplets that enables them to move more independently from the background flow field. There is, however, an exception in the case without aerodynamic interactions and with gravity, which shows a reduction of RRV for droplets larger than $45\ \mathrm {\mu }$m. This phenomenon was reported in several former studies, see e.g. Rosa et al. (Reference Rosa, Parishani, Ayala, Grabowski and Wang2013), Rosa et al. (Reference Rosa, Parishani, Ayala and Wang2015) and Ireland, Bragg & Collins (Reference Ireland, Bragg and Collins2016). Its physics is non-trivial, since for monodisperse systems there is no differential sedimentation and the effect of gravity on the particle-pair motion is entirely implicit. In contrast to bidisperse systems (Dhariwal & Bragg Reference Dhariwal and Bragg2018) the effect of particle settling does not offset the effect of turbulent mixing. Ireland et al. (Reference Ireland, Bragg and Collins2016) presented physical arguments explaining this mechanism, namely, gravity suppresses the non-local contribution to the particle relative velocities by reducing the correlation time scale of the flow experienced by the particles. Therefore, the average relative velocity of the droplets is lower than in simulations without gravity. The effect of gravity on the RRV is negligible for droplets smaller than $30\ \mathrm {\mu }$m. Moreover, for such small droplets, at $\phi = 1.22\times 10^{-2}$ or lower, the effect of aerodynamic interaction may be practically neglected. This is likely due to stronger correlation of such particles with the turbulent flow. The motion of low inertia particles is dominated by the turbulent flow and the aerodynamic interactions among them do not have a considerable impact.

Figure 13.

Comparison of at-contact RRV in case when there is aerodynamic interaction, both including and excluding lubrication forces, and when there is no aerodynamic interaction, all with and without gravity, for (a) the same number of droplets, $N_{p} = 50\ 000$ and (b) the same mass loading, $\phi = 1.22\times 10^{-2}$ at the dissipation rate $\epsilon = 400$ cm$^2$ s$^{-3}$.

Including aerodynamic interactions into the model has a direct impact on the kinematic collision statistics making the RRV more closely dependent on the RDF and vice versa. This can be explained with an example, namely, an enhanced RRV increases collision rate and consequently enforces redistribution of particles. In simulations without aerodynamic interactions particles can overlap, so this redistribution is not needed. To quantify this phenomenon the changes in RRVs have to be examined together with changes in RDFs. Figure 14 shows the RDF complementary to the RRV (plotted in figure 13) calculated in the same series of simulation. As expected, the largest values of the RDF were obtained in simulations without aerodynamic interactions. In turn, in simulations with aerodynamic interactions the RDFs depend on the mass loading. For systems at a fixed number of droplets, see figure 14(a), the differences between simulations with lubrication forces and standard HDNS appear for droplets larger than $30\ \mathrm {\mu }$m. According to these results there is stronger clustering at short length scales when the effect of lubrication is taken into account, especially together with gravity. In simulations at a fixed mass loading, the offset in the RDF is almost constant for all considered droplet radii.

Figure 14.

Comparison of at-contact RDF when there is aerodynamic interaction, both including and excluding lubrication forces, and when there is no aerodynamic interaction, all with and without gravity, for (a) the same number of droplets, $N_{p} = 50\ 000$ and (b) the same mass loading, $\phi = 1.22\times 10^{-2}$ at the dissipation rate $\epsilon = 400$ cm$^2$ s$^{-3}$. The empirical model of Ayala, Rosa & Wang (Reference Ayala, Rosa and Wang2008a, (84) and (85) for $R_{\lambda }=76.86$) is also included for comparison.

Figure 15 presents the dynamic collision kernels of the cases considered above for $\epsilon = 400$ cm$^2$ s$^{-3}$. In general, the collision kernel is a monotonic and increasing function of the droplet radii. This increase of $\varGamma ^{D}$ is mainly due to larger values of RRVs which is a consequence of larger particle inertia. The lubrication forces slightly reduce the collision kernel, which again stems from the reduction in the relative velocities. The RRVs are smaller when the lubrication force between pairs is taken into account. For each size of droplets, the lowest collision kernel is usually for sedimenting droplets and the highest for non-sedimenting droplets if there is no aerodynamic interaction. When the aerodynamic interaction is considered, the collision kernels lie in-between, being slightly lower for the sedimenting droplets. This is not general, however, since the mass loading is important, too. Accordingly, the results for the mass loading $\phi = 1.22\times 10^{-2}$ are shown in figure 15(c). The effect of mass loading on the collision statistics of droplets will be discussed in more detail in § 4.3.

Figure 15.

Comparison of the dynamic collision kernel when there is aerodynamic interaction, both including and excluding lubrication forces, and when there is no aerodynamic interaction, all with and without gravity, for (a) the same number of droplets, $N_{p} = 50\ 000$ in linear scale and in (b) logarithmic scale and for (c) the same mass loading, $\phi = 1.22\times 10^{-2}$ at the dissipation rate $\epsilon = 400$ cm$^2$ s$^{-3}$.

4.2. Validation of corrected kinematic formulations

The RRVs and RDFs presented in this study are subjected to the corrections (3.4)–(3.7) provided by Wang et al. (Reference Wang, Ayala, Kasprzak and Grabowski2005b). In figure 16, the collision kernels computed using the corrected kinematic statistics at contact, (3.8), are benchmarked against the dynamic collision kernels, (3.3), computed in the simulations with and without gravity. In simulations without aerodynamic interactions, the correction is not needed because particles can overlap, hence the volumetric rate of collisions is equal to the influx of particles from the surface area of the collision sphere. Consequently, $\varGamma ^{D}$ and $\varGamma ^{K}$ in figure 16 perfectly overlap.

Figure 16.

Comparisons between the collision kernel obtained from the dynamic formulation with that computed using the kinematic formulation, both before (UC) and after (C) applying corrections, when (a) gravity is not considered and (b) when there is gravity.

However, for the cases where there is aerodynamic interaction, the relocation mechanism changes the influx through the collision sphere by creating a particle-free zone behind it. In such a case, standard kinematic formulation leads to an inconsistency with the dynamic formulations. In figure 16 the kinematic collision kernels are presented both before (uncorrected (UC)) and after (corrected (C)) applying the corrections. We find that after correction the kinematic and dynamic collision kernels are in fairly good quantitative agreement. It should be noted, however, that there is a small but visible difference between $\varGamma ^{K}$ and $\varGamma ^{D}$ in simulations performed with full representation of the lubrication forces. We hypothesise that this discrepancy may have two sources. First, numerical methods have some intrinsic limitations, so that the values of RDF and RRV at contact can be evaluated with a finite accuracy. In this particular case, the RDF may be somewhat overestimated. Second, the average statistics of these simulations were calculated using fewer samples. As shown in figure 15, $\varGamma ^{D}$ is consistently smaller in simulations with lubrication forces compared with standard HDNS. This means that particles collide less frequently, which in turn might result in larger numerical uncertainty in these statistics. Overall, the investigation shows that these geometric corrections are necessary to include in simulations when the particles are displaced after collisions.

4.3. Effect of particles number density and mass loading

Finally, we examine the effect of droplet mass loading on the kinematic and dynamic collision statistics. Figure 17 shows at-contact RRVs and RDFs together with corresponding values of the dynamic collision kernels computed in simulations with different number of droplets. The considered systems were monodisperse and contained droplets of radii between $20\text {--}60\ \mathrm {\mu }$m.

Figure 17.

Changes in the at-contact (a) RRV, (b) at-contact RDF and (c) the dynamic collision kernel for the sets of droplets with different radii and numbers, i.e. different mass loadings.

The results presented in figure 17(a) confirm observations reported in former studies, e.g. (Rosa et al. Reference Rosa, Parishani, Ayala, Grabowski and Wang2013), that inertia enhances the RRV between closely separated particles. It stems from the fact that particles of larger inertia retain their kinetic energy longer than corresponding fluid elements, therefore their relative motion is more strongly decorrelated. Gravity noticeably reduces the velocity correlation between fluid and particles and thus diminishes the non-local effects (Ireland et al. Reference Ireland, Bragg and Collins2016). Consequently, the RRV of larger droplets is significantly lower. In simulations with $N_{p} = 20\,000$ of $60\ \mathrm {\mu }$m droplets, gravity reduces the RRV by approximately three times. This is the regime in which the effect of the aerodynamic interactions is relatively weak. In simulations with larger number of droplets or equivalently larger mass loadings the aerodynamic interactions gain importance and alter dynamics of the entire systems. As the mean distance between the droplets decreases the additional forces result in a stronger decorrelation of the particle motion but mainly in simulations with gravity. Therefore, we observe a monotonic increase in RRV with the particle number density. The slopes of these dashed lines decrease with the size of the droplets, which is a consequence of smaller mass loadings (the simulations were performed at fixed $N_{p}$). Thus, the effect of aerodynamic interactions is weaker. Interestingly, this effect is much smaller in simulations without gravity. This is rather expected since the RRV in simulations without gravity is already very large such that the aerodynamic interaction does not result in further enhancement of the particle relative motion.

An opposite, i.e. decreasing, trend is observed for the RDFs visualised in figure 17(b). For most simulated cases the RDF decreases with $N_{p}$. The only exceptions are systems of the middle-sized 30 and 40 $\mathrm {\mu }$m droplets. However, the little differences fall in the range of statistical uncertainty. The decreasing trend is more pronounced for larger settling droplets. This is attributed to the enhanced relative velocities, which prevents droplet clustering, at least at distances comparable to the particle scale.

The two opposite trends present in the kinematic statistics should be also reflected in the collision kernels. Figure 17(c) shows $\varGamma ^{D}$ as a function of the droplet number density. It indicates that the relative changes in the collision kernel for increasing $N_{p}$ are smaller compared with the changes in the kinematic statistics. This likely results from the fact that $\varGamma ^{D}$ is proportional to the product of the RRV and RDF. As expected, the largest sensitivity to $N_{p}$ is observed for the largest settling droplets. Interestingly, $\varGamma ^{D}$ of $60\ \mathrm {\mu }$m droplets is not monotonic. For low droplet concentrations, i.e. up to $N_{p} = 60\ 000$, there is a systematic increase in the collision kernel, which is the effect of enhancement in RRVs. Then, for $N_{p} > 60\ 000$ the dynamic collision kernel decreases, and this corresponds to the same decreasing trend in the RDF.