On the effect of lubrication forces on the collision statistics of cloud droplets in homogeneous isotropic turbulence

Abstract

We investigate the dynamics of inertial particles in homogeneous isotropic turbulence, under one-way momentum coupling, using a new computational approach that incorporates the effect of long-range many-body aerodynamic interactions along with the short-range lubrication forces. The implementation couples hybrid direct numerical simulations (HDNS) with the analytical solutions of two rigid spheres moving in an unbounded fluid. Concerning the velocity field seen by the particles, the algorithm switches from the flow solution in terms of HDNS to analytical formulae when the separation distance between particles becomes comparable to their average radius. Standard HDNS is unable to correctly represent the short-range interactions since this method is based on the superposition of the Stokes solutions for single spheres. Our results show that for the turbulent kinetic energy dissipation rates typical of atmospheric clouds, the radial relative velocities (RRVs) of the droplets increase, and the radial distribution function (RDF) decreases in the near-contact region if the lubrication forces are taken into account. These changes are more pronounced when the effect of gravity is considered. Away from the contact region, there is not much change in RRVs and RDFs. For turbulent clouds with lower dissipation rates lubrication forces significantly enhance the average RRV in the limit of low Stokes number. This enhancement, however, is statistically insignificant because the number of particle pairs at close proximity is very small. The effect of mass loading on the collision statistics is also investigated, demonstrating an increase in RRV and a reduction in RDF with the droplet concentration.

1. Introduction

Transport of inertial particles by turbulent flows is a common phenomenon observed in nature. The examples include dust storms, dunes formation, transport of organic matter in water reservoirs, or airborne particles, such as pollen, soot and ashes. Turbulent transport is also an important process in a variety of industrial applications, e.g. filtering of solid particles, propulsion systems and oil processing. In this study, we focus on quantitative analysis of cloud microphysical processes in turbulent air. Accurate description of these complex atmospheric phenomena is central for reliable weather and climate predictions on Earth. Interaction of cloud droplets or ice crystals with the turbulent flow as well as interactions among different hydrometeors have a direct impact on the precipitation formation, i.e. the rate and amount of precipitation. The characteristic length scales of these microphysical processes are significantly smaller than those defining large-scale atmospheric flows; therefore, they cannot be resolved in numerical weather prediction (NWP) models. Nowadays, NWP models provide regular forecasts at horizontal resolutions varying from coarse resolutions, of several kilometres, up to one kilometre (Yano et al. 2018). In the standard NWP approach, the effect of cloud processes at unresolved scales is usually accounted for by parameterisation, representing only some statistical features of the modelled systems. To develop more realistic parameterisations, a detailed knowledge of the physics underlying these processes is necessary.

The importance of the cloud microphysical process resulted in a rich scientific literature with a particular focus on quantifying the growth of cloud droplets in the sizes ranging from 10 to 50 $\mathrm {\mu }$m in radius. This range is commonly called the size gap, for which neither the condensation nor the gravitational collision-coalescence mechanism is effective (Pruppacher & Klett 1997; Chen et al. 2018). A number of studies have been carried out to explain the growth of such droplets as well as the fast broadening of their size spectra. In these studies several mechanisms have been investigated such as growth by ultragiant particles (Yin et al. 2000; Van Den Heever & Cotton 2007), entrainment-induced spectral broadening (Blyth 1993), effects of pre-existing clouds, and enhanced collision rate by turbulence (Devenish et al. 2012; Grabowski & Wang 2013; Rosa et al. 2013). In recent years, the most intensely analysed and discussed has been the mechanism of turbulent collision-coalescence (Wang et al. 2008; Rosa et al. 2011c, 2013; Grabowski & Wang 2013; Rosa et al. 2016).

Air turbulence can increase the rate of droplet collisions in several different ways. It can raise the probability of collisions by increasing the relative velocity between droplet pairs as a result of shear effects and differential acceleration. In addition, droplet trajectories are biased towards the regions of lower vorticity and higher strain rate. Consequently, turbulence can increase, in average, the clustering in the distribution of droplets, thus bringing droplets closer together and increasing the probability of collision. Moreover, if the effect of gravity is taken into consideration, the interaction of particles with the turbulent flow will selectively alter their settling velocity, consequently enhancing the collision rate. Furthermore, turbulence affects the local aerodynamic interactions among particles by changing their relative motion, local acceleration and shear effects.

In most previous numerical studies of cloud processes, the continuous phase was modelled using the Eulerian approach, employing direct numerical simulations (DNS) or large eddy simulations (LES), e.g. (Rosa & Pozorski 2017), while the dispersed phase was treated using the Lagrangian approach together with the point-particle assumption and one-way momentum coupling between the continuous and dispersed phases (Wang et al. 2009; Rosa et al. 2013). This simplified approach will be sufficient only if the particles do not significantly affect the motion of the continuous phase. For larger concentrations, i.e. mass loadings, of the dispersed phase, the effect of turbulence modulation induced by particles becomes important. There are several studies where the effect of two-way momentum coupling was considered (Bosse, Kleiser & Meiburg 2006; Monchaux & Dejoan 2017; Rosa, Pozorski & Wang 2020a,b). It should be emphasised, however, that the effect of aerodynamic interaction in simulations with two-way coupling has been considered only to a limited extent. Full representation of lubrication forces, to the best of our knowledge, has so far never been included in the particle equations of motion. As for the fully resolved simulations of turbulence with finite-size particles, the two-way momentum coupling and particle aerodynamic interactions are automatically accounted for. Such studies have become possible only in recent years (Devenish et al. 2012; Gao, Li & Wang 2013; Kuerten 2016; Wang et al. 2016a,b), yet, they are limited to a relatively small number of particles, typically $O(10^5)$.

An innovative method for treating the aerodynamic interactions among cloud droplets was introduced by Wang, Ayala & Grabowski (2005a). Their improved superposition method (ISM) was tested in a set of idealised experiments aimed at computing the collision efficiency of small water droplets settling under gravity in stagnant air. Compared with the original superposition method (Pruppacher & Klett 1997), ISM is more accurate. For example, in simulations where the effect of lubrication is not dominant the relative error on the drag force can be reduced by an order of magnitude. The original formulation of the superposition method fails to satisfy the no-slip boundary condition on the surfaces of the interacting spheres, therefore, errors in the representations of drag forces can be large. Application of ISM to a many-body system of mutually interacting droplets in a turbulent flow is rather straightforward, leading to the hybrid direct numerical simulation (HDNS) approach (Ayala, Grabowski & Wang 2007). The basic idea of HDNS is to combine the DNS of the background turbulence with an analytical representation of the disturbance flows introduced by particles. This approach takes advantage of the fact that the disturbance flows due to droplets are localised in space and there is a sufficient length scale separation between the droplet size and the Kolmogorov scale of the background flow. In HDNS, the disturbance flow experienced by each droplet due to the presence of other droplets is derived from a linear system of $3N_{p}$ unknowns, namely three components of velocity in each spatial direction, where $N_{p}$ is the number of droplets in the simulation. Hybrid DNS is a significant step forward in modelling cloud processes. Nevertheless, the iterative solution of the large linear system is computationally expensive.

The first HDNS implementation was developed based on the OpenMP parallel library (Ayala et al. 2007) and therefore it could not take the full advantage of modern machines with distributed memory. This limited the scalability of the code, and hence, the simulations could be run using low resolution meshes, equivalent to low Reynolds numbers, and small numbers of particles. To increase the scalability of HDNS several attempts have been made (Rosa & Wang 2009; Rosa et al. 2011a; Ayala et al. 2014), the most recent of which (Torres et al. 2013) is a massively parallel implementation based on two-dimensional domain decomposition that use the MPI library for data communication. Under the HDNS framework, Onishi, Takahashi & Vassilicos (2013) proposed a more efficient approach named binary interaction-based superposition method, suitable for dilute systems. Compared with HDNS, calculating the hydrodynamic interactions in systems typical of atmospheric clouds requires an order of magnitude less central processing unit (CPU) time. Therefore, the method is capable of performing simulations with larger number of particles and using larger domains. As for the accuracy of the method, the error is comparable to that of HDNS if the volume fraction of the disperse phase is low.

All the above-mentioned methods are based on the Stokes solution for a single sphere and therefore cannot correctly handle short-range or lubrication forces. Other studies in which the simple superposition method was used (Shafrir & Neiburger 1963; Shafrir & Gal-Chen 1971; Lin & Lee 1975; Pinsky, Khain & Shapiro 2001) resulted in a rough prediction of the collision efficiency. Therefore, rigorously incorporating the lubrication effect is necessary in order to obtain reliable data to quantify the microphysical processes.

There are several studies aimed at investigating the lubrication forces between spherical particles. The first rigorous analytical solution for an idealised system of two coaxial spheres moving at equal small constant velocities along their line of centres was developed by Stimson & Jeffery (1926). Another interesting example, here, is the solution elaborated by Jeffrey & Onishi (1984). They considered an unbounded low-Reynolds-number viscous flow enforced by translational and rotational motion of two unequal rigid spheres. The forces and torques acting on the particles were derived in terms of an infinite series making use of the twin multipole expansions method. While relatively easy for numerical implementation, their solutions are valid for all separation distances between spheres and accurately represent the asymptotic behaviour of the forces when the gap distance between them goes to zero, i.e. the lubrication effects. Very recently, Goddard, Mills-Williams & Sun (2020) used a bipolar coordinate system to develop an alternative representation of interacting forces for a squeezing flow of two unequal rigid spheres and a shearing flow of two equal-sized spheres. They claimed that for different size particles their solution is more precise compared with the asymptotic solutions of other existing methods. There are also some studies in which lubrication interactions among three spherical particles are considered, e.g. Cichocki, Ekiel-Jeżewska & Wajnryb (1999). Also, the dynamics as well as stability of three spheres sedimenting under gravity was examined by Ekiel-Jeżewska & Wajnryb (2006). Unfortunately, the computational cost of such simulations is so large that the method cannot be directly used for modelling systems with many particles.

In 2011, Rosa et al. (2011a,b) developed a computationally efficient approach for treating the aerodynamic interactions between two spheres moving in still air. Compared with the exact solution of Jeffrey & Onishi (1984), the forces and torques were accurately represented in the entire range of separation distances. The method combines the force–torque–stresslet (FTS) multipole expansion of Durlofsky, Brady & Bossis (1987) for the long-range interaction, with the lubrication expansion of Jeffrey & Onishi (1984) at short separation distances. For the intermediate range where no simple method is available, a third-order polynomial fitting was used to bridge the treatments for long-range and short-range interactions. This composite approach resulted in a large reduction in the computational cost when compared with the full treatment of Jeffrey & Onishi (1984), while keeping the relative error in calculated collision efficiency typically within only $O(1\,\%)$. This approach, while numerically efficient, allows the handling of the interaction between two particles, not among many particles. Hence, further advancement is needed in order to adopt the method for computing collision statistics of many droplets in turbulent flows.

Accordingly, the objective of this study is to investigate dynamics of cloud droplets under one-way momentum coupling in homogeneous isotropic turbulence (HIT), considering full representation of interaction forces, i.e. including the effects of lubrication. To achieve this goal, the exact analytical solution developed by Jeffrey & Onishi (1984) has been incorporated into the numerical code for HDNS. The simulations of turbulent flows are handled using the standard pseudo-spectral method, as detailed in § 2. The hybrid approach is rigorous and convenient for studying collision statistics, preferential concentration and collision efficiency of cloud droplets.

The content of this paper is structured as follows. In § 2, we describe the numerical tools that are used to simulate homogeneous isotropic turbulence and particle tracking. Results from HDNS are discussed in § 3. These include sensitivity of the kinematic and dynamic collision statistics to the time step size and location of matching point. This matching point indicates the separation distance at which standard HDNS switches to the short-range analytical solution. Extended analysis of the radial distribution function (RDF) and radial relative velocity (RRV) at different droplet mass loadings is presented in § 4. Section 5 contains a summary and main conclusions.

2. Methodology

2.1. Background turbulent flow

Following the previous studies (Torres et al. 2013; Ayala et al. 2014), 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 (1994). 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 1983) 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. (2011a). 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. 2007), 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. (2013) 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. 2013). 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. (2007) 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. 2005b; Ayala et al. 2007; Wang, Ayala & Grabowski 2007; Wang et al. 2008, 2009). 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. 2005a) 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 1984) 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 1996). Comparison of these analytically predicted forces with those obtained from the ISM (Wang et al. 2005a) 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 (1984). This strategy is somewhat similar to the one proposed by Durlofsky et al. (1987, (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 (1984).

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 1984)

(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 (1984). 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. 1992). 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 1996). 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 (1970). 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. Collision statistics – sensitivity of numerical results to ${\rm \Delta} t$ and $\delta$

The new implementation was used to compute collision statistics of monodisperse particle systems relevant to atmospheric clouds. The droplet radii considered in these simulations varied between $20\text {--}60\ \mathrm {\mu }$m. The main focus is on two-point particle statistics such as the radial distribution function and the average radial relative velocity. These quantities were computed ‘on the fly’ during simulations and then averaged over time at postprocessing. The RRV between two particles is defined as $w_{r}=\boldsymbol {w\cdot r}/r$, where $\boldsymbol {w}=\boldsymbol {V}^{(2)}-\boldsymbol {V}^{(1)}$ is their relative velocity vector. This quantity is a function of separation distance $r$ and can be computed for every pair of droplets in the computational domain. In further analysis, we only consider the absolute value of the average RRV, i.e. $\langle |w_{r} (r/R)|\rangle$ for all particle pairs computed for a sufficiently long simulation time; in most cases 100$T_e$. The spherical neighbourhood that is scanned around each droplet is limited to $1\leq r/R\leq 10$. This range is divided into 180 spherical shells of thickness $\delta _{sh}=0.05R=0.1a$.

A similar method was used to compute the RDF. The RDF characterises preferential concentration, or clustering (Bragg, Ireland & Collins 2015), of the particles at different length scales. For a system with a monodisperse set of droplets, RDF is defined as

(3.1)\begin{equation} g_{r} (r/R; t)=\frac{n_{pair}/V_{sh}}{\frac{1}{2}N_{p}(N_{p}-1)/V_{box}}, \end{equation}

where $V_{box}$ is the volume of the computational domain and $V_{sh}(r/R)$ is the volume of each spherical shell within which fall $n_{pair}(r/R; t)$ pairs of droplets.

The product of the RRV and RDF of aerodynamically non-interacting particles at contact, i.e. $r=R$, is proportional to the kinematic collision kernel, $\varGamma ^{K}$, namely

(3.2)\begin{equation} \varGamma^{K}=2{\rm \pi} R^2 \langle|w_{r} (r=R)|\rangle g_{r} (r=R).\end{equation}

In NWP, the collision kernel is central for mathematical models that represent grid scale cloud processes and precipitation formation. An alternative formulation for $\varGamma$, the so-called dynamic collision kernel, can be evaluated by counting collisions of droplets, $N_c(t)$, in the entire computational box for a given time period. Having the volumetric rate of collisions, $\dot {N}_c$, the dynamic collision kernel, $\varGamma ^{D}$, can be evaluated using the following formula (Ayala et al. 2008b):

(3.3)\begin{equation} \dot{N}_c=\varGamma^{D} \tfrac{1}{2}N_{p}(N_{p}-1)/V_{box}.\end{equation}

The assessment of uncertainty in the dynamic collision kernel is explained in (Rosa et al. 2013). The dynamic formulation of the collision kernel is based on the actual collision events and therefore is more accurate. Besides, the effect of aerodynamic interactions is inherently accounted for, while the standard definition of kinematic collision kernel (3.2) needs to be revised for that purpose. This is largely due to the non-overlapping condition applied to the pairs of droplets that collide. After each collision one of the droplets must be relocated, resulting in a particle-free zone behind the collision sphere of the collided pair. To address this issue, Wang et al. (2005b) proposed a correction based on subtracting the influx through the surface area of each shell that is in the shadow of the collision sphere from the total influx. Such a treatment is justified since no particles can be found behind the collision sphere due to the imposed relocation condition. Consequently, the resulting correction factors for $g_{r} (r/R)$ and $\langle |w_{r} (r/R)|\rangle$ in each spherical shell, with inner and outer radii of $r_1/R$ and $r_2/R$, are

(3.4)\begin{gather} C_g\left(\frac{r_1}{R},\frac{r_2}{R}\right)=\frac{1}{2}+ \frac{1}{2}\frac{((r_2/R)^2-1 )^{3/2}-((r_1/R)^2-1 )^{3/2}}{(r_2/R)^3 -(r_1/R)^3}, \end{gather}

(3.5)\begin{gather}C_w\left(\frac{r_1}{R},\frac{r_2}{R}\right)=\left.\left( 1-\frac{3}{2}\frac{r_2/R-r_1/R}{(r_2/R)^3-(r_1/R)^3} \right) \right/ C_g\left(\frac{r_1}{R},\frac{r_2}{R}\right), \end{gather}

respectively, such that the corrected RDF and RRV will be obtained by

(3.6)\begin{gather} g_{r}^{c} (r/R) =g_{r}/C_g, \end{gather}

(3.7)\begin{gather}\langle|w_{r}^{c} (r/R)|\rangle =\langle|w_{r}|\rangle/C_w \end{gather}

and the kinematic collision kernel

(3.8)\begin{equation} \varGamma^{K}=2{\rm \pi} R^2 \langle|w_{r}^{c} (r=R)|\rangle g_{r}^{c} (r=R).\end{equation}

From now on, the superscript $c$ for corrected RRV and RDF will be dropped.

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 (1970) and Rosa et al. (2011c), 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 (1984) 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 (1984) 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. Effects of aerodynamic interactions on kinematic and dynamic collision statistics

In this section, we discuss more broadly the effect of lubrication interaction on the droplet collision statistics. Particularly, we analyse the role of lubrication forces in systems with droplets of different inertia. The range of droplet radii considered in these simulations varies from 20 to 60 $\mathrm {\mu }$m. Moreover, the simulations have been performed for a wide range of energy dissipation rates: $10\text {--}400$ cm$^2$ s$^{-3}$. The results from the new implementation of HDNS with lubrication forces are compared with the simulations from the standard HDNS, those without aerodynamic interactions, and the results from other studies as well as theoretical estimations. First, we address the effect of the energy dissipation rate on the kinematic and dynamic collision statistics in the absence of gravity. Figure 10 shows the at-contact RDFs and RRVs of monodisperse systems computed in six series of simulations each at a different $\epsilon$. In all these simulations the number of droplets was fixed and equal to $N_{p} = 50\ 000$. Additionally, the results are plotted against the Stokes number, which is a function of both the size of the droplet and turbulent energy dissipation rate. The Stokes number is a quantitative measure of the droplet inertia. With such scaling all the lines collapse to one U-shaped line with a minimum roughly at $St = 0.2$. The three additional curves shown there correspond to: (i) the current implementation for the same set of droplets at $\epsilon = 100\ \mathrm {cm}^2\ \mathrm {s}^{-3}$ but without aerodynamic interactions; (ii) the data from the study of Rosa et al. (2013, figure 13a), and values from the theoretical model by Saffman & Turner (1956) valid in the limiting case of zero-inertia particles. The results in figure 10(a) show that the RRVs are highly sensitive to $\epsilon$ in the limit of both very small and large droplets, but in a different manner. As for the large droplets, e.g. 60 $\mathrm {\mu }$m, this increase can be explained as the effect of larger inertia. The motion of such droplets is weakly correlated with the background flow, and consequently the RRV is larger.

Figure 10. The effect of turbulent energy dissipation rate on the at-contact RRV for different (a) droplet sizes and (b) Stokes numbers, and its effect on the at-contact radial distribution function for different (c) droplet sizes and (d) Stokes numbers. In panel (b) the theoretical model of Saffman & Turner (1956) is shown. The DNS results of Rosa et al. (2013, figures 8a and 13a there) are additionally included in panels (b) and (d).

More interestingly, an opposite trend is observed for the smaller droplets. The average RRV at contact increases when the energy dissipation rate decreases. In figure 10(b), for $St < 0.2$ there is a sharp increase in the RRVs since only pairs of droplets with considerably large relative velocities are able to break through the barrier of lubrication forces. This, however, does not happen frequently, as can be deduced from figure 10(d). The values of the RDF at contact for $St < 0.2$ are very low, which means that only a small number of droplets can become nearly touching. Figure 10(c) shows that the values of at-contact RDFs computed using different $\epsilon$ follow a bell-shaped line. As the energy dissipation rate decreases, the bell stretches and the maximum shifts towards the larger droplets. In figure 10(d) the RDFs are plotted also as a function of the Stokes number, exhibiting a maximum approximately at $St = 0.8$.

The key conclusion, however, concerns differences between the standard simulations under one-way momentum coupling and those with full representation of aerodynamic forces. As for the RRV we observe a certain similarity of both models in the range $St>0.2$. A minor increase in RRV for droplets with $St<1.9$ can be seen, as well as a little reduction for droplets of larger inertia if the interaction forces are considered. For $St<0.2$ the differences between different approaches are large, but statistically insignificant, due to rarity of situations when two particles are close to each other. The lubrication forces as well as the long range many-body interactions have significant effects on the RDF. For droplets of $St=0.8$ the difference in the RDF reaches almost 50 %. There are two reasons for this reduction. First, the lubrication forces prevent the particles from approaching each other. Second, in the simulations with aerodynamic interactions after each collision one of the droplets is relocated, decreasing local droplet number density, whereas in simulations under one-way momentum coupling droplets are allowed to overlap. In addition to the two-point collision statistics, the dynamic collision kernel, $\varGamma ^{D}$, is computed for the same sets of simulations and presented in figure 11. According to the data in figure 11(a,c), the collision kernel increases monotonically both with the droplet radius and the energy dissipation rate. Since both parameters correspond to larger Stokes numbers, or equivalently larger inertia, the droplets collide more frequently due to higher relative velocities. The normalised values of the dynamic collision kernel are presented in figure 11(b,d) as a function of the Stokes number. Within a numerical uncertainty, all values collapse to a single curve showing an increase with the Stokes number until $St = 1.6$. For larger $St$, a slight reduction is noticeable. The increasing trend is largely due to the augmentation of relative velocities with the Stokes number, see figure 10(b), while the decreasing trend stems from the reduction in the clustering of droplets, see figure 10(d).

Figure 11. The dynamic collision kernel as function of droplet size in flows with different energy dissipation rates for (a) $N_{p} = 50\ 000$ and (c) $\phi = 1.22\times 10^{-2}$, and the normalised dynamic collision kernel for corresponding values of the Stokes number, also for (b) $N_{p} = 50\ 000$ and (d) $\phi = 1.22\times 10^{-2}$. The results of Rosa et al. (2013) are added to panel (b) for comparison.

Finally, we address the differences between simulations with and without aerodynamic interactions. Figure 11(b) shows a comparison of $\varGamma ^{D}$ computed using the new model with results obtained from simulations under one-way momentum coupling. Based on the data, it can be claimed that the collision kernel, and consequently the collision rate of monodisperse systems, are smaller if the aerodynamic interactions are considered. The reduction is generally a function of the particle inertia. Specifically, for droplets of $St=1$ the reduction is approximately 25 %. The above results together suggest that this effect should be parameterised in future NWP models.

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 (1997) and Yavuz et al. (2018), 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. (2013), Rosa et al. (2015) and Ireland, Bragg & Collins (2016). 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 2018) the effect of particle settling does not offset the effect of turbulent mixing. Ireland et al. (2016) 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 (2008a, (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. (2005b). 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. 2013), 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. 2016). 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.

5. Conclusions

The collision-coalescence of water droplets has been quantified using the combined Eulerian–Lagrangian numerical approach incorporating the exact representation of aerodynamic interactions between droplets. The numerical method combines the standard HDNS approach with the exact analytical formula for the lubrication forces. Merging of these two approaches allows us to jointly represent two important effects, usually neglected in previous studies. The HDNS accurately captures the effect of many-body interactions among widely separated droplets. In turn, the analytical solution for two interacting rigid spheres in low-Reynolds-number flows allows us to represent the lubrication effects. The impact of gravity was also considered, and the simulations have been performed for both sedimenting droplets and droplets without sedimentation. For modelling background turbulent flows, we used a computational mesh of relatively low resolution, i.e. $64^3$ (low Reynolds number). A reasonable numerical cost of such simulations enables us to obtain particle statistics for a wide range of parameters. The main focus of this study was on the effects of aerodynamic interactions on the kinematic and dynamic collision statistics rather than sensitivity to $R_{\lambda }$.

Several important conclusions can be drawn based on the results from numerical experiments. First, in the absence of gravity the aerodynamic interactions reduce the RDF, and the reduction depends on the particle inertia (see figures 4, 10d and 12b). For example, for droplets with $St=0.8$ the RDF is approximately 50 % smaller than in analogous simulations without aerodynamic interactions. Second, the RRV is slightly enhanced, particularly in those systems with droplets of low inertia (see figures 3, 10b and 12a). Consequently, the reduction of the collision kernels of interacting droplets is also a function of the Stokes number (see figure 11b). Specifically, for droplets of $St = 1$, the reduction is approximately 25 %. Third, the effect of lubrication forces is more pronounced in systems with larger energy dissipation rates, especially if gravitational settling is considered (compare figure 12 with figures 13–15). This conclusion is drawn based on the comparison of standard HDNS with simulations performed using the new approach with full representations of aerodynamic interactions. Fourth, once the mass loading grows, the kinematic collision statistics reveal an opposite trend, namely the RRV increases with the droplet number density, while the RDF monotonically decreases (compare figure 17a with 17b). This is the combined effect of many-body interaction, gravity and the numerical assumption made in which one of the droplets is relocated after collision.

The present results are a step forward in quantifying the effect of lubrication interaction on the collision statistics but are limited to the monodisperse systems and low mesh resolutions, or equivalently low Taylor-microscale flow Reynolds numbers. Therefore, a potentially important direction of further research may be an extension of this analysis to polydisperse systems, covering a wide range of droplet sizes. In doing so, the need for random relocation after collisions is eliminated yielding more physically realistic results. Also, simulations at higher mesh resolutions may shed more light on the effect of aerodynamic interactions. It should be noted, however, that the numerical cost of simulations with larger domains is disproportionately much higher. This leap in complexity is mainly due to the need to track a large number of droplets necessary to maintain a desired mass loading. Also note that including lubrication interactions imposes additional constraints on the integration time step. Another direction possibly worthy of further investigation is the effect of lubrication when the continuum assumption of the fluid is no longer applicable. At such small separation distances the lubrication forces are lower than those predicted under the continuum fluid assumption.