Study of Colliding Particle-Pair Velocity Correlation in Homogeneous Isotropic Turbulence

: This paper deals with the numerical analysis of the particle inertia and volume fraction e ﬀ ects on colliding particle-pair velocity correlation immersed in an unsteady isotropic homogeneous turbulent ﬂow. Such correlation function is required to build reliable statistical models for inter-particle collisions, in the frame of the Euler–Lagrange approach, to be used in a broad range of two-phase ﬂow applications. Computations of the turbulent ﬂow have been carried out by means of Direct Numerical Simulation (DNS) by the Lattice Boltzmann Method (LBM). Moreover, the dependence of statistical properties of collisions on particle inertia and volumetric fraction is evaluated and quantiﬁed. It has been found that collision locations of particles of intermediate inertia, St K ∼ 1, occurs in regions where the ﬂuid strain rate and dissipation are higher than the corresponding averaged values at particle positions. Connected with this fact, the average kinetic energy of colliding particles of intermediate inertia (i.e., Stokes number around 1) is lower than the value averaged over all particles. From the study of the particle-pair velocity correlation, it has been demonstrated that the colliding particle-pair velocity correlation function cannot be approximated by the Eulerian particle-pair correlation, obtained by theoretical approaches, as particle separation tends to zero, a fact related with the larger values of the relative radial velocity between colliding particles. been for particles of and intermediate


Introduction
Inter-particle collision of small particles immersed in a turbulent flow is a common micro-process found in both industrial and natural environments. Examples of physical phenomena where the effects of inter-particle interactions are revealed as relevant are: models of planet formation from dust particle collisions and coalescence in turbulent protoplanetary nebulae [1,2], rain drops' formation in warm clouds [3,4], chemical reaction rates enhancement due to active suspended particles in fluidised beds [5][6][7], particle transport and reduction of wall erosion in pneumatic conveying [8][9][10][11] and turbomachinery [12], or sand storm evolution and dune transport in deserts [13,14].
Fundamental ingredients of inter-particle collisions are the role of particles' inertia and their interaction with the carrier turbulent flow structures, which lead to a non-linear increase of collision rates [15][16][17][18]. As a result of the efforts to understand and quantify the enhancement of collision rates induced by inertia, two main mechanisms have been identified: on the one hand, preferential concentration, i.e., particle clustering, increases the probability for two particles to be at a colliding distance, and, on the other hand, detachment from fluid trajectories may enhance the relative velocity between two approaching particles. Two limiting cases can be devised: in the case of vanishingly employed successfully for predicting the effect of inter-particle collisions in pneumatic conveying and particle separation processes [9][10][11]35]. Lain et al. [36] developed a more rigorous way to determine the Eulerian particle-pair velocity correlation function between colliding particles in homogeneous isotropic turbulence, employing the model proposed in Zaichik and Alipchenkov [37] (or its refined version, Zaichik et al. [18]), showing good agreement with previous data [19]. Precisely, the adequacy of such correlation is evaluated and discussed in the present contribution in light of the new DNS results obtained in this work.
This paper deals with the numerical analysis of the collision statistics of point particles immersed in an unsteady isotropic homogeneous turbulent flow depending on particle inertia and volume fraction. The main goal is to determine the two-particle velocity correlation function of colliding particles, information that is appropriated to be used in inter-particle collision stochastic models for use in general-purpose two-phase flow Euler-Lagrange codes. Moreover, the effect of particle volume fraction on collisions' statistics is illustrated, especially on the relative motion of colliding particles. The present article is organised as follows: Section 2 presents a context and a brief summary of the performed Direct Numerical Simulations based on the lattice Boltzmann method combined with the point-particle Lagrangian tracking. The main collision statistics as functions of particle inertia and void fraction are presented and discussed in Section 3. Section 4 is devoted to the study of the two-particle velocity correlation function of colliding particles, which is connected with the relative velocity statistics. The obtained dependence of it with volume fraction as a parameter is used to build a model for velocity correlation of colliding particles appropriate for inter-particle collision stochastic models. Finally, the main conclusions and perspectives are presented in Section 5.

Context and Summary of Numerical Simulations
The simplest and most extensively studied type of turbulence is statistically homogeneous isotropic turbulence of an incompressible fluid. In such configuration, it is well known that small-scale fields of velocity can be thought of as being more or less homogeneous and isotropic. Small-scale vortex structures are responsible for dissipation of turbulent energy and play a key role in the accumulation (clustering) of particles in a turbulent flow. Study of isotropic turbulence is thus of fundamental importance for both single-phase and two-phase flows and, therefore, it has been the configuration adopted in the present work.
In this work, the simulation of the turbulent behaviour of the flow has been performed by the Lattice Boltzmann Method with the in-house code described in Reference [38] and also employed in Reference [38] to study the effect of inter-particle collisions on several clustering measures of inertial particles in turbulent flows. For the details of the implementation of the Lattice Boltzmann Method (LBM) forcing scheme for generating the homogeneous isotropic turbulence and collision detection algorithm, the reader is referred to Reference [38]. Therefore, only an abridged summary of the main flow characteristics is mentioned here.
We focus on the dynamics of colliding particles in an isotropic homogeneous flow in a triple-periodic box of size L = 0.128 m. The fluid domain is uniformly discretised by 128 3 grid points. Forcing parameters have been adjusted to guarantee a good resolution of all flow scales, providing a grid space slightly smaller than the Kolmogorov length scale η K . Additionally, particle size is much smaller than η K so they are treated as point particles; moreover, the simulated particle volume fractions are sufficiently small to keep the two-phase flow dilute, so the back-reaction force of the particles on the flow is neglected. For the same reason, only binary elastic collisions are considered in the framework of the hard sphere approach. A brief summary of the employed collision model and detection algorithm is presented in Appendix A. Particle trajectories are computed by the simplified particle equation of motion, retaining only the contribution of drag term, in combination with the Schiller and Naumann correlation [39,40]. Domain discretisation and particle volume fraction values have been chosen to keep the computational effort limited as the direct detection of particle collisions is very costly. The main parameters employed in the fluid simulation are presented in Table 1. As an illustration, Figure 1 shows the non-dimensional three-dimensional (3D) energy spectrum ε κ obtained in the present DNS. For comparison, the experimental results of Comte-Bellot and Corrsin [41] and the Kolmogorov spectrum are also presented.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 22 direct detection of particle collisions is very costly. The main parameters employed in the fluid simulation are presented in Table 1. As an illustration, Figure 1 shows the non-dimensional three-dimensional (3D) energy spectrum obtained in the present DNS. For comparison, the experimental results of Comte-Bellot and Corrsin [41] and the Kolmogorov spectrum are also presented. The considered particles have a fixed diameter, = 200 μm, which is more than five times smaller than the Kolmogorov length scale ; therefore, particle size is well in the dissipative range.   The considered particles have a fixed diameter, d p = 200 µm, which is more than five times smaller than the Kolmogorov length scale η K ; therefore, particle size is well in the dissipative range. Three different average particle volume fractions, α p = 0.01%, 0.05% and 0.1%, have been simulated. Particle volume fraction is varied by changing the number of tracked particles: for α p = 0.01% there are 50,066 particles, for α p = 0.05% there are 250,329 particles and for α p = 0.1%, there are 500,658 particles. Particle properties are summarised in Table 2.
Inertial particles' dynamics depend on Stokes number, which characterises particle inertia and is equal to the ratio between the average response time of particles, τ p , and the Kolmogorov time The response time scale of the particle is defined as: where τ pS = ρ p d 2 p /18µ is the Stokesian particle relaxation time, µ is the carrier fluid dynamic viscosity, c D is the particle drag coefficient which depends on the particle Reynolds number, Re p = ρd p u − u p /µ, whose expression extended for Reynolds number higher than 0.5 is given in Reference [40] as: whereas the Kolmogorov time microscale is expressed as: In this study, variation of particle Stokes number is achieved by modifying its density according to the values shown in Table 2. For quantifying the importance of inter-particle collisions on disperse phase transport, a collision Stokes number can be defined as St col = τ c / τ p , where τ c is the mean inter-particle collision time, i.e., the average time between two consecutive collisions. As mentioned in Reference [42], St col allows drawing a boundary between dilute and dense flow conditions: if τ c is larger than τ p , the flow is ditute (i.e, St col > 1), and dense otherwise. In the case of dilute flows, the main particle transport mechanism is driven by its interaction with the turbulence, but in case of dense systems, collisions play a role in particle motion.
Although previous works have considered the effect of particle number density on inter-particle collision frequency (e.g., Reference [15]), it was in the article of Ernst and Sommerfeld [38] where a systematic study of the effect of particle volume fraction on collision statistics was performed. The employed DNS resolution of that work was low, only 64 3 , and therefore the grid size was a factor two larger than the Kolmogorov length scale. Nevertheless, some interesting results were found: as expected, inter-particle collision frequency augmented with increasing volume fraction and Stokes number, the Lagrangian integral time scale of particles decreased with volume fraction but increased with St K due to the decorrelation of particle motion induced by the increase in the number of collisions, also, the particle pair relative velocity at collision rose with growing St K but it decreased for higher values of α p , with this last effect being more pronounced for large inertia particles. Such findings are revised in the next section with the results obtained in the present, more refined, DNS.
Finally, our previous work (Ernst et al. [39]) studied the effect of particle volume fraction and collisions on particle segregation and clustering in homogeneous isotropic turbulence in the same configuration and numerical set up as those employed in the present study. It was found that for particle volume fractions lower than 0.1%, the effect of increasing collision frequency with α p was marginal as all evaluated clustering measures (segregation parameter, correlation dimension, radial distribution, Voronoï diagrams and Minkowski functionals) were only slightly modified; however, for higher volume fractions, of about 1%, such measures were noticeably altered due to the re-dispersing effect of inter-particle collisions.

Statistical Properties of Inter-Particle Collisions
As an illustration, Figure 2 shows typical snapshots of a two-dimensional slice of particle distribution in the turbulent flow field depending on their inertia. Following Reference [39], particles tend to form distinct patterns in the flow, which depend on its inertia. For instance, for the case of α p = 0.1%, Figure 2a presents the case of St K = 1, where well-differentiated narrow particle ropes can be observed in the background coloured by fluid strain rate (average strain rate in this configuration is 61 s −1 ). As it is clearly noticed, particles tend to preferentially sample zones with a higher than mean fluid strain rate (green to red colours), and also, although not shown, lower than mean fluid vorticitya view that agrees with the previous literature. Figure 2b shows, on the other hand, the distribution of particles with St K = 10; here, no distinct clusters can be appreciated, and particles tend to be spread much more evenly than in Figure 2a, although some voids are noticed. As commented before, the reason for such distribution patterns lies in the particle-turbulence interaction: for intermediate particle inertia, St K ∼ 1, the centrifuging effect of particles towards the rim of vortical structures as well as the sling effect impulse them to show preferential concentration; on the other side, large inertia particles are only marginally affected by its interaction with turbulent vortices and tend to describe straight trajectories (the so-called inertia filtering effect), which are only modified by collisions, resulting in a much more uniform distribution along the fluid field.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 22 = 0.1%, Figure 2a presents the case of = 1, where well-differentiated narrow particle ropes can be observed in the background coloured by fluid strain rate (average strain rate in this configuration is 61 s −1 ). As it is clearly noticed, particles tend to preferentially sample zones with a higher than mean fluid strain rate (green to red colours), and also, although not shown, lower than mean fluid vorticity-a view that agrees with the previous literature. Figure 2b shows, on the other hand, the distribution of particles with = 10; here, no distinct clusters can be appreciated, and particles tend to be spread much more evenly than in Figure 2a, although some voids are noticed. As commented before, the reason for such distribution patterns lies in the particle-turbulence interaction: for intermediate particle inertia, ~1, the centrifuging effect of particles towards the rim of vortical structures as well as the sling effect impulse them to show preferential concentration; on the other side, large inertia particles are only marginally affected by its interaction with turbulent vortices and tend to describe straight trajectories (the so-called inertia filtering effect), which are only modified by collisions, resulting in a much more uniform distribution along the fluid field. respectively. The diameter of the spheres is enlarged by a factor or 3 for the sake of better visibility. Figure 3 illustrates the clustering behaviour of inertial particles as a function of their Stokes number by means of the so-called particle segregation parameter Σ [39]. It compares the standard deviation of the actual particle number density distribution with that of a purely random distribution (which is described by a Poisson process). The value Σ = 0 corresponds to a random particle distribution, Σ < 0 corresponds to a uniform particle distribution, while Σ > 0 indicates particle clustering.  The diameter of the spheres is enlarged by a factor or 3 for the sake of better visibility. Figure 3 illustrates the clustering behaviour of inertial particles as a function of their Stokes number by means of the so-called particle segregation parameter Σ p [39]. It compares the standard deviation of the actual particle number density distribution with that of a purely random distribution (which is described by a Poisson process). The value Σ p = 0 corresponds to a random particle distribution, Σ p < 0 corresponds to a uniform particle distribution, while Σ p > 0 indicates particle clustering.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 22 = 0.1%, Figure 2a presents the case of = 1, where well-differentiated narrow particle ropes can be observed in the background coloured by fluid strain rate (average strain rate in this configuration is 61 s −1 ). As it is clearly noticed, particles tend to preferentially sample zones with a higher than mean fluid strain rate (green to red colours), and also, although not shown, lower than mean fluid vorticity-a view that agrees with the previous literature. Figure 2b shows, on the other hand, the distribution of particles with = 10; here, no distinct clusters can be appreciated, and particles tend to be spread much more evenly than in Figure 2a, although some voids are noticed. As commented before, the reason for such distribution patterns lies in the particle-turbulence interaction: for intermediate particle inertia, ~1, the centrifuging effect of particles towards the rim of vortical structures as well as the sling effect impulse them to show preferential concentration; on the other side, large inertia particles are only marginally affected by its interaction with turbulent vortices and tend to describe straight trajectories (the so-called inertia filtering effect), which are only modified by collisions, resulting in a much more uniform distribution along the fluid field. respectively. The diameter of the spheres is enlarged by a factor or 3 for the sake of better visibility. Figure 3 illustrates the clustering behaviour of inertial particles as a function of their Stokes number by means of the so-called particle segregation parameter Σ [39]. It compares the standard deviation of the actual particle number density distribution with that of a purely random distribution (which is described by a Poisson process). The value Σ = 0 corresponds to a random particle distribution, Σ < 0 corresponds to a uniform particle distribution, while Σ > 0 indicates particle clustering.  As it is clearly seen from Figure 3, Σ p is positive for all Stokes numbers, showing a peak for St K around one and values closer to zero for Stokes numbers of 0.1 and 10, where particles tend to be distributed randomly along the fluid domain.
In the following, we present the effects of increasing particle volume fraction on the statistical collision properties of the particles in dependence on the Stokes number, which implies taking a second look at the results of Reference [38] in light of the present, more refined, DNS. It should be mentioned that the averaging period for computing the presented collision statistics was larger than 25 eddy turnover times (defined as L int /u ).

Analysis of Time Scales Related with Collisions
The inter-particle collision time, τ c , is defined as the average of elapsed time between two consecutive collisions and it is computed as the inverse of inter-particle collision frequency (which indicates the number of collisions per unit time). Figure 4 presents such quantity for the three considered volume fractions, α p , and it has been made non-dimensional with the Kolmogorov time scale, τ K .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 22 As it is clearly seen from Figure 3, Σ is positive for all Stokes numbers, showing a peak for around one and values closer to zero for Stokes numbers of 0.1 and 10, where particles tend to be distributed randomly along the fluid domain.
In the following, we present the effects of increasing particle volume fraction on the statistical collision properties of the particles in dependence on the Stokes number, which implies taking a second look at the results of Reference [38] in light of the present, more refined, DNS. It should be mentioned that the averaging period for computing the presented collision statistics was larger than 25 eddy turnover times (defined as / ′).

Analysis of Time Scales Related with Collisions
The inter-particle collision time, , is defined as the average of elapsed time between two consecutive collisions and it is computed as the inverse of inter-particle collision frequency (which indicates the number of collisions per unit time). Figure 4 presents such quantity for the three considered volume fractions, , and it has been made non-dimensional with the Kolmogorov time scale, . As shown in Figure 4, inter-particle collision time is a decreasing function of the Stokes number, although for the largest Stokes numbers considered, it seems to reach an asymptotic value. This trend is similar to those presented in References [22] and [38]. Moreover, for a fixed value of , decreases almost linearly with increasing volume fraction, which is verified for all the considered Stokes number range. In the case of the smallest = 0.1, there are many repeated collisions among the same particle pair with very small radial relative velocity and nearly parallel velocities. In order to exclude possible false collision detections due to numerical errors in such small , only collision intervals larger than the Kolmogorov time scale are counted in the statistics [22].
Another point to be mentioned is that for the smallest values of = 0.1 and 0.5, not all the particles are involved in collisions. In fact, for the Stokes number of 0.1, at = 10 , only 17% of all tracked particles are experiencing at least one collision; in this case, particles tend to be nearly homogeneously distributed along the flow field and inter-particle collisions happen between particle pairs that have shared the same flow structures during long periods. As a result, the velocities of both particles are correlated, and the collision mechanism is that of fluid velocity shear. For = 0.5, where particle accumulation is already noticeable, the number of particles involved in collisions increases up to 65%, whereas for the rest of the higher , such value is over 94%. Therefore, essentially all particles are participating in collisions for the cases with intermediate and large inertia.
In the next figure, the collision Stokes number, , is plotted as a function of the particle Stokes number for the three volume fractions ( Figure 5). As shown in Figure 4, inter-particle collision time is a decreasing function of the Stokes number, although for the largest Stokes numbers considered, it seems to reach an asymptotic value. This trend is similar to those presented in References [22] and [38]. Moreover, for a fixed value of St K , τ c decreases almost linearly with increasing volume fraction, which is verified for all the considered Stokes number range. In the case of the smallest St K = 0.1, there are many repeated collisions among the same particle pair with very small radial relative velocity and nearly parallel velocities. In order to exclude possible false collision detections due to numerical errors in such small St K , only collision intervals larger than the Kolmogorov time scale are counted in the statistics [22].
Another point to be mentioned is that for the smallest values of St K = 0.1 and 0.5, not all the particles are involved in collisions. In fact, for the Stokes number of 0.1, at α p = 10 −3 , only 17% of all tracked particles are experiencing at least one collision; in this case, particles tend to be nearly homogeneously distributed along the flow field and inter-particle collisions happen between particle pairs that have shared the same flow structures during long periods. As a result, the velocities of both particles are correlated, and the collision mechanism is that of fluid velocity shear. For St K = 0.5, where particle accumulation is already noticeable, the number of particles involved in collisions increases up to 65%, whereas for the rest of the higher St K , such value is over 94%. Therefore, essentially all particles are participating in collisions for the cases with intermediate and large inertia.
In the next figure, the collision Stokes number, St col , is plotted as a function of the particle Stokes number for the three volume fractions ( Figure 5). Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 22 It can be seen that, for all particle volume fractions and for all inertias, the flow can be considered dilute, as the collision Stokes number is always above unity. A pronounced decreasing of is observed for the three considered volume fractions and the curves run pretty much parallel, suggesting a linear dependence of with for a fixed particle Stokes number. The collision Stokes number decreases continuously with , following a nearly straight line in the double logarithmic diagram for the three values of particle volume fraction considered. Only for the most inertial particles did the values of become near the dilute-dense flow conditions boundary; however, in such case, the increment of the number of collisions is mainly ascribed to the increase of the caustics contribution (or random uncorrelated motion) to inter-particle collision rate [28].
According to Figure 5, all the considered cases are in the dilute regime, where 〈 〉 < ; therefore, the main mechanism affecting the particles' motion is the aerodynamic transport, and particles have time to adapt to its local fluid environment between inter-particle collisions. Increasing has the effect of reducing inter-particle distance, so collisions are more frequent. However, such collisions are not dominant, they must be viewed instead as a modulating mechanism. This fact was demonstrated in Reference [39], where the cases with and without inter-particle collisions were compared: the differences between both cases in the considered clustering measures were marginal for < 0.1%. As a consequence, the increase in collision frequency with is expected to have only a moderate effect on the statistical quantities computed for colliding particle events, such as relative velocity, two-particle velocity correlation function or strain rate at collision locations, among others.
The importance of collisional effects on the turbulence-induced particle transport can be illustrated by plotting the behaviour of the particle Lagrangian integral time scale, , versus particle inertia.
is computed as the integral of the particle Lagrangian autocorrelation function and it measures the time scale in which particle velocities are correlated in time along their trajectories [38]. Figure 6a shows the curves of the quotient / depending on particle Stokes number, , for the three considered values of . A decreasing trend with increasing is observed, which is due to the growing of the time scale with particle inertia, reflecting the fact that more inertial particles are less affected by turbulent fluctuations and that the particle motion becomes more correlated. As a consequence, large inertia particles tend to experience a ballistic motion between consecutive collisions. Regarding the effect of particle volume fraction on , it is found that the increase in the number of collisions has only a marginal effect on the particle Lagrangian autocorrelation timescale, where, nevertheless, a small decrease of this quantity can be observed with increasing (i.e., increase in / ), as expected. This means that the particle motion becomes slightly more decorrelated as a consequence of the increasing collision rate and the resulting faster change of particle velocity. It can be seen that, for all particle volume fractions and for all inertias, the flow can be considered dilute, as the collision Stokes number is always above unity. A pronounced decreasing of St col is observed for the three considered volume fractions and the curves run pretty much parallel, suggesting a linear dependence of St col with α p for a fixed particle Stokes number. The collision Stokes number decreases continuously with St K , following a nearly straight line in the double logarithmic diagram for the three values of particle volume fraction considered. Only for the most inertial particles did the values of St col become near the dilute-dense flow conditions boundary; however, in such case, the increment of the number of collisions is mainly ascribed to the increase of the caustics contribution (or random uncorrelated motion) to inter-particle collision rate [28].
According to Figure 5, all the considered cases are in the dilute regime, where τ p < τ c ; therefore, the main mechanism affecting the particles' motion is the aerodynamic transport, and particles have time to adapt to its local fluid environment between inter-particle collisions. Increasing α p has the effect of reducing inter-particle distance, so collisions are more frequent. However, such collisions are not dominant, they must be viewed instead as a modulating mechanism. This fact was demonstrated in Reference [39], where the cases with and without inter-particle collisions were compared: the differences between both cases in the considered clustering measures were marginal for α p < 0.1%. As a consequence, the increase in collision frequency with α p is expected to have only a moderate effect on the statistical quantities computed for colliding particle events, such as relative velocity, two-particle velocity correlation function or strain rate at collision locations, among others.
The importance of collisional effects on the turbulence-induced particle transport can be illustrated by plotting the behaviour of the particle Lagrangian integral time scale, τ p L , versus particle inertia. τ p L is computed as the integral of the particle Lagrangian autocorrelation function and it measures the time scale in which particle velocities are correlated in time along their trajectories [38]. Figure 6a shows the curves of the quotient τ K /τ p L depending on particle Stokes number, St K , for the three considered values of α p . A decreasing trend with increasing St K is observed, which is due to the growing of the time scale τ p L with particle inertia, reflecting the fact that more inertial particles are less affected by turbulent fluctuations and that the particle motion becomes more correlated. As a consequence, large inertia particles tend to experience a ballistic motion between consecutive collisions. Regarding the effect of particle volume fraction on τ p L , it is found that the increase in the number of collisions has only a marginal effect on the particle Lagrangian autocorrelation timescale, where, nevertheless, a small decrease of this quantity can be observed with increasing α p (i.e., increase in τ K /τ p L ), as expected. This means that the particle motion becomes slightly more decorrelated as a consequence of the increasing collision rate and the resulting faster change of particle velocity. . Kolmogorov time scale over particle Lagrangian time scale (a) and inter-particle collision time over particle Lagrangian time scale (b) as a function of particle Stokes number for the three evaluated average volume fractions. Figure 6b presents the behaviour of / versus particle inertia. The same as in Figure 6a, a decreasing tendency is found for this variable, which is related to the fact that increases and decreases with growing . For all the computed cases, the quotient / is above unity, meaning that inter-particle collision time is larger than the particle Lagrangian time scale. Similar to what happened with , the curves for the different particle volume fractions decrease pretty fast with increasing inertia and run quasi parallel. Moreover, the values of / are always much larger than those of / , which means that the main mechanism of particle velocity change is the interaction between particle inertia and underlying turbulence. Inter-particle collisions have a marginal effect modulating such primary mechanism. Figure 7 illustrates the variation of the mean free path, , with particle inertia and with particle volume fraction. The free path is computed as the distance travelled by a particle among two consecutive collisions. Particles with Stokes numbers well below one mainly interact with their closest neighbours, with whom they have shared the same fluid environment during long periods. This fact is reflected in a low relative velocity and collision angles between them (see Figure 8 below). Besides, as such particles do not experience remarkable clustering: interactions between the same pair of particles are common and only a small percentage of them undergo collisions (for = 0.1 and = 10 , only 17% of particles are involved in collisions with a partner and the percentage is even lower for smaller values of , 9% for 5 • 10 and 2% for 10 ). This means that the majority of particles present a free path equal to the length of the trajectory described during the total measured time, which in this case was 25 eddy turnover times. As a result, small particles present high values of mean free paths. As increases, the percentage of particles involved in collisions augments (65% for = 0.5 and = 10 ) at the same time that the number of repeated collisions among the same pair of particles is reduced; as a consequence, the mean free path reduces. For around unity and larger, nearly all the particles are involved in collisions, a fact favoured by clustering processes and the formation of caustics [22]. Both phenomena increase collision frequency and also reduce the value of the mean free path. Finally, for large enough , particles become nearly irresponsive to fluid fluctuations and do not undergo appreciable clustering, so the main mechanism for inter-particle collisions is the so-called Random Uncorrelated Motion [19]. In this case, the number of collisions increases slightly, and the mean free path seems to reach an asymptotic value. Obviously, the effect of decreasing is reflected in a reduction of the inter-particle collision frequency, hence an increase of the mean free path, as shown in Figure 7.  Figure 6b presents the behaviour of τ c /τ p L versus particle inertia. The same as in Figure 6a, a decreasing tendency is found for this variable, which is related to the fact that τ p L increases and τ c decreases with growing St K . For all the computed cases, the quotient τ c /τ p L is above unity, meaning that inter-particle collision time is larger than the particle Lagrangian time scale. Similar to what happened with St col , the curves for the different particle volume fractions decrease pretty fast with increasing inertia and run quasi parallel. Moreover, the values of τ c /τ p L are always much larger than those of τ K /τ p L , which means that the main mechanism of particle velocity change is the interaction between particle inertia and underlying turbulence. Inter-particle collisions have a marginal effect modulating such primary mechanism. Figure 7 illustrates the variation of the mean free path, λ p , with particle inertia and with particle volume fraction. The free path is computed as the distance travelled by a particle among two consecutive collisions. Particles with Stokes numbers well below one mainly interact with their closest neighbours, with whom they have shared the same fluid environment during long periods. This fact is reflected in a low relative velocity and collision angles between them (see Figure 8 below). Besides, as such particles do not experience remarkable clustering: interactions between the same pair of particles are common and only a small percentage of them undergo collisions (for St K = 0.1 and α p = 10 −3 , only 17% of particles are involved in collisions with a partner and the percentage is even lower for smaller values of α p , 9% for 5 × 10 −4 and 2% for 10 −4 ). This means that the majority of particles present a free path equal to the length of the trajectory described during the total measured time, which in this case was 25 eddy turnover times. As a result, small particles present high values of mean free paths. As St K increases, the percentage of particles involved in collisions augments (65% for St K = 0.5 and α p = 10 −3 ) at the same time that the number of repeated collisions among the same pair of particles is reduced; as a consequence, the mean free path reduces. For St K around unity and larger, nearly all the particles are involved in collisions, a fact favoured by clustering processes and the formation of caustics [22]. Both phenomena increase collision frequency and also reduce the value of the mean free path. Finally, for large enough St K , particles become nearly irresponsive to fluid fluctuations and do not undergo appreciable clustering, so the main mechanism for inter-particle collisions is the so-called Random Uncorrelated Motion [19]. In this case, the number of collisions increases slightly, and the mean free path seems to reach an asymptotic value. Obviously, the effect of decreasing α p is reflected in a reduction of the inter-particle collision frequency, hence an increase of the mean free path, as shown in Figure 7. Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 22  Figure 8a shows that average velocity collision angles, defined as the angle that particle velocities form at collision, increase continuously with . The same behaviour is seen for the relative velocity at collision (Figure 8b), which is a fact related with the increased uncorrelated motion of colliding particles with growing inertia. Therefore, particles can collide coming from rather different directions, a phenomenon that is enhanced by trajectory deflection after previous collisions. As mentioned before, particles with small inertia respond very fast to the fluid velocity fluctuations and colliding pairs are usually moving with pretty similar velocities in nearly parallel trajectories [22] because they are transported by the same eddy during long time periods, as a result, their velocities are very correlated, which is reflected in both small values of relative velocity and nearly zero-velocity angle of collision. In this context, there is no preferential concentration and the increased collision frequency, i.e., by augmenting , has no effect in the particle transport and collision mechanisms which are driven by the turbulent velocity fluctuations.

Analysis of Kinematic Properties of Colliding Particles
With increasing particle inertia, preferential concentration and the sling effects come into play, with the result of increasing the collision frequency, relative velocity at contact as well as velocity angle of collision. In this regime, the effect of increasing the particle volume fraction is hardly seen as the inter-particle collision time is still much higher than the particle Lagrangian time scale (see Figure 6b). Interestingly, for = 1, the mean velocity collision angle obtained in the present work, around 0.3 radians or 17°, is very close to that obtained in Choi et al. [22]. For larger inertia, the relative velocity of colliding pairs overcomes the Kolmogorov velocity and reaches a maximum in the present computations for = 10, the same as the velocity collision angle (around 1.17 radians or 67°), which again is very similar to that reported in Reference [22]. When inertia is further  Figure 8a shows that average velocity collision angles, defined as the angle that particle velocities form at collision, increase continuously with . The same behaviour is seen for the relative velocity at collision (Figure 8b), which is a fact related with the increased uncorrelated motion of colliding particles with growing inertia. Therefore, particles can collide coming from rather different directions, a phenomenon that is enhanced by trajectory deflection after previous collisions. As mentioned before, particles with small inertia respond very fast to the fluid velocity fluctuations and colliding pairs are usually moving with pretty similar velocities in nearly parallel trajectories [22] because they are transported by the same eddy during long time periods, as a result, their velocities are very correlated, which is reflected in both small values of relative velocity and nearly zero-velocity angle of collision. In this context, there is no preferential concentration and the increased collision frequency, i.e., by augmenting , has no effect in the particle transport and collision mechanisms which are driven by the turbulent velocity fluctuations.
With increasing particle inertia, preferential concentration and the sling effects come into play, with the result of increasing the collision frequency, relative velocity at contact as well as velocity angle of collision. In this regime, the effect of increasing the particle volume fraction is hardly seen as the inter-particle collision time is still much higher than the particle Lagrangian time scale (see Figure 6b). Interestingly, for = 1, the mean velocity collision angle obtained in the present work, around 0.3 radians or 17°, is very close to that obtained in Choi et al. [22]. For larger inertia, the relative velocity of colliding pairs overcomes the Kolmogorov velocity and reaches a maximum in the present computations for = 10, the same as the velocity collision angle (around 1.17 radians or 67°), which again is very similar to that reported in Reference [22]. When inertia is further  Figure 8a shows that average velocity collision angles, defined as the angle that particle velocities form at collision, increase continuously with St K . The same behaviour is seen for the relative velocity at collision (Figure 8b), which is a fact related with the increased uncorrelated motion of colliding particles with growing inertia. Therefore, particles can collide coming from rather different directions, a phenomenon that is enhanced by trajectory deflection after previous collisions. As mentioned before, particles with small inertia respond very fast to the fluid velocity fluctuations and colliding pairs are usually moving with pretty similar velocities in nearly parallel trajectories [22] because they are transported by the same eddy during long time periods, as a result, their velocities are very correlated, which is reflected in both small values of relative velocity and nearly zero-velocity angle of collision. In this context, there is no preferential concentration and the increased collision frequency, i.e., by augmenting α p , has no effect in the particle transport and collision mechanisms which are driven by the turbulent velocity fluctuations.
With increasing particle inertia, preferential concentration and the sling effects come into play, with the result of increasing the collision frequency, relative velocity at contact as well as velocity angle of collision. In this regime, the effect of increasing the particle volume fraction is hardly seen as the inter-particle collision time is still much higher than the particle Lagrangian time scale (see Figure 6b). Interestingly, for St K = 1, the mean velocity collision angle obtained in the present work, around 0.3 radians or 17 • , is very close to that obtained in Choi et al. [22]. For larger inertia, the relative velocity of colliding pairs overcomes the Kolmogorov velocity and reaches a maximum in the present computations for St K = 10, the same as the velocity collision angle (around 1.17 radians or 67 • ), which again is very similar to that reported in Reference [22]. When inertia is further increased, the relative velocity of colliding pairs decreases because particles become less responsive to fluid fluctuations [37] and the velocity collision angle tends asymptotically to 90 • . For these large inertias, a marginal effect of α p starts to be appreciated as both relative velocity and collision angle are slightly reduced for increasing collision frequency. This fact was already noticed in Reference [38] and its explanation was ascribed to collective effects in clusters. Actually, the argument is that in situations where clustering is notable, particles spend most of the time in areas of high strain ratio and low vorticity so, besides some collisional effects, they have time to adapt its motion to that of the fluid, resulting finally in a reduced average relative particle velocity and collision angle. However, such explanation would not be appropriate for large inertia particles because in those conditions, particles are distributed fairly uniformly without appreciable clustering. This finding should be further investigated by using configurations with larger volume fractions and turbulent Reynolds numbers.
An interesting issue to be discussed is the distribution of the location of inter-particle collisions. As already mentioned, and illustrated in Figure 2a, particles with Stokes number around 1, experiencing appreciable clustering, tend to sample fluid locations with higher strain rate than particles with small and large Stokes numbers, where clustering is not observed. Figure 9a presents the ratio between the fluid strain rate at particle collision locations, averaged for all collisional events, S f @p p,coll , and the fluid strain rate at particle positions, averaged for all particles, S f @p p . Interestingly, inter-particle collisions at St K ∼ 1 clearly happen in fluid regions with noticeably higher strain rates than those averaged through all the particles. The maximum occurs for St K = 1, whereas for St K = 0.1, 10, the ratios are very close to 1, as in such cases no significant preferential concentration is observed. At this point, it is necessary to mention that Choi et al. [22] found that collision events in the case of particles with St K ∼ 1 occurred in locations where the fluid dissipation was higher than the corresponding mean value experienced by all the tracked particles. Both observations are related, as in the configuration of homogeneous isotropic turbulence, both quantities, fluid strain rate and dissipation, are strongly correlated (correlation coefficient of 0.97 in present simulations).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 22 increased, the relative velocity of colliding pairs decreases because particles become less responsive to fluid fluctuations [37] and the velocity collision angle tends asymptotically to 90°. For these large inertias, a marginal effect of starts to be appreciated as both relative velocity and collision angle are slightly reduced for increasing collision frequency. This fact was already noticed in Reference [38] and its explanation was ascribed to collective effects in clusters. Actually, the argument is that in situations where clustering is notable, particles spend most of the time in areas of high strain ratio and low vorticity so, besides some collisional effects, they have time to adapt its motion to that of the fluid, resulting finally in a reduced average relative particle velocity and collision angle. However, such explanation would not be appropriate for large inertia particles because in those conditions, particles are distributed fairly uniformly without appreciable clustering. This finding should be further investigated by using configurations with larger volume fractions and turbulent Reynolds numbers.
An interesting issue to be discussed is the distribution of the location of inter-particle collisions. As already mentioned, and illustrated in Figure 2a, particles with Stokes number around 1, experiencing appreciable clustering, tend to sample fluid locations with higher strain rate than particles with small and large Stokes numbers, where clustering is not observed. Figure 9a presents the ratio between the fluid strain rate at particle collision locations, averaged for all collisional events, 〈 @ 〉 , , and the fluid strain rate at particle positions, averaged for all particles, 〈 @ 〉 . Interestingly, inter-particle collisions at ~1 clearly happen in fluid regions with noticeably higher strain rates than those averaged through all the particles. The maximum occurs for = 1, whereas for = 0.1, 10, the ratios are very close to 1, as in such cases no significant preferential concentration is observed. At this point, it is necessary to mention that Choi et al. [22] found that collision events in the case of particles with ~1 occurred in locations where the fluid dissipation was higher than the corresponding mean value experienced by all the tracked particles. Both observations are related, as in the configuration of homogeneous isotropic turbulence, both quantities, fluid strain rate and dissipation, are strongly correlated (correlation coefficient of 0.97 in present simulations).  Figure 9b presents the ratio between the fluid vorticity at particle collision locations, averaged for all collisional events, 〈 @ 〉 , , and the fluid vorticity at particle positions, averaged for all particles, 〈 @ 〉 . As it was said previously, particles ~1 tend to sample fluid locations with lower than mean vorticity. In our simulations, this behaviour also happens for the smallest inertia particles, = 0.1, which otherwise approach the tracers' behaviour. They are located in areas with lower than mean fluid vorticity, as already noted in Reference [22]. Looking at the locations where inter-particle collisions happen, 〈 @ 〉 , , it is found that such locations are characterised by  Figure 9b presents the ratio between the fluid vorticity at particle collision locations, averaged for all collisional events, ω f @p p,coll , and the fluid vorticity at particle positions, averaged for all particles, ω f @p p . As it was said previously, particles St K ∼ 1 tend to sample fluid locations with lower than mean vorticity. In our simulations, this behaviour also happens for the smallest inertia particles, St K = 0.1, which otherwise approach the tracers' behaviour. They are located in areas with lower than mean fluid vorticity, as already noted in Reference [22]. Looking at the locations where inter-particle collisions happen, ω f @p p,coll , it is found that such locations are characterised by smaller values than ω f @p p (see Figure 9b). This observation is valid for all the Stokes numbers, although for particles with St K ≈ 10, the two values are very close. This fact is parallel to the tendency of particle collisions to occur in areas with fluid strain rate higher than the mean value computed at particle positions, S f @p p , and, therefore, could be somehow expected. However, this last observation is contrary to that presented in Reference [22], where the fluid enstrophy sampled at collision locations was higher than the average value sampled at particle positions. At the moment, the reason for such discrepancy is not clear.
Finally, Figure 10 shows the ratio of the average particle turbulent energy of colliding particles, k p,coll , and the average particle turbulent energy of all the tracked particles, k p , k p,coll /k p . In this case, colliding particles with St K ∼ 1 have smaller kinetic energy than the average, attaining a minimum of around 0.9 for St K = 1 and for all considered particle volume fractions. For particles with St K = 0.1, the ratio is closer to 1, although the value for α p = 10 −4 is higher than for the other two volume fractions. The fact that colliding particles with intermediate inertia, St K ∼ 1, possess on average less fluctuating energy than the mean value computed over all tracked particles could be understood as follows: as it has been shown before, inter-particle collisions tend to happen in locations with high values of fluid strain rate and, equivalently, high dissipation rate; in isotropic homogeneous turbulence, the decay of turbulent kinetic energy is larger at locations of higher dissipation, therefore, such regions possess, on average, lower values of turbulent energy. Consequently, it is expected that intermediate inertia particle collisions occur at fluid locations with lower than mean values of fluid kinetic energy, a fact that is found in the present simulations. Moreover, as it has been seen for particles of low and intermediate inertia, that St col values are much larger than 1, so those particles have time to accommodate to the local fluid conditions between collisions, meaning that their kinetic energy will also be lower than the mean kinetic energy computed for all tracked particles.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 22 smaller values than 〈 @ 〉 (see Figure 9b). This observation is valid for all the Stokes numbers, although for particles with ≈ 10, the two values are very close. This fact is parallel to the tendency of particle collisions to occur in areas with fluid strain rate higher than the mean value computed at particle positions, 〈 @ 〉 , and, therefore, could be somehow expected. However, this last observation is contrary to that presented in Reference [22], where the fluid enstrophy sampled at collision locations was higher than the average value sampled at particle positions. At the moment, the reason for such discrepancy is not clear. Finally, Figure 10 shows the ratio of the average particle turbulent energy of colliding particles, , , and the average particle turbulent energy of all the tracked particles, , , / . In this case, colliding particles with ~1 have smaller kinetic energy than the average, attaining a minimum of around 0.9 for = 1 and for all considered particle volume fractions. For particles with = 0.1, the ratio is closer to 1, although the value for = 10 is higher than for the other two volume fractions. The fact that colliding particles with intermediate inertia, ~1, possess on average less fluctuating energy than the mean value computed over all tracked particles could be understood as follows: as it has been shown before, inter-particle collisions tend to happen in locations with high values of fluid strain rate and, equivalently, high dissipation rate; in isotropic homogeneous turbulence, the decay of turbulent kinetic energy is larger at locations of higher dissipation, therefore, such regions possess, on average, lower values of turbulent energy. Consequently, it is expected that intermediate inertia particle collisions occur at fluid locations with lower than mean values of fluid kinetic energy, a fact that is found in the present simulations. Moreover, as it has been seen for particles of low and intermediate inertia, that values are much larger than 1, so those particles have time to accommodate to the local fluid conditions between collisions, meaning that their kinetic energy will also be lower than the mean kinetic energy computed for all tracked particles. Figure 10. Average of particle kinetic energy of colliding particles over mean particle kinetic energy over all particle locations, , / , as a function of particle Stokes number for the three evaluated average volume fractions.
In the case of large inertia particles, the situation is different as they are not preferentially sampling regions of high dissipation, and also, the collision Stokes number is only somewhat higher than 1. As a result, in the case of = 10, the , / ratio is clearly above 1 for all volume fractions, indicating in this case that colliding particle pairs possess higher fluctuating velocity than the mean. According to the random uncorrelated motion interpretation, such result is expected, as inter-particle collisions are favoured by a higher particle fluctuating velocity.

Collision Frequency
A way of characterising inter-particle collision probability is the collision kernel, which is equal to two times the collision frequency over the particle number density (particles per unit volume) Figure 10. Average of particle kinetic energy of colliding particles over mean particle kinetic energy over all particle locations, k p,coll /k p , as a function of particle Stokes number for the three evaluated average volume fractions.
In the case of large inertia particles, the situation is different as they are not preferentially sampling regions of high dissipation, and also, the collision Stokes number is only somewhat higher than 1. As a result, in the case of St K = 10, the k p,coll /k p ratio is clearly above 1 for all volume fractions, indicating in this case that colliding particle pairs possess higher fluctuating velocity than the mean. According to the random uncorrelated motion interpretation, such result is expected, as inter-particle collisions are favoured by a higher particle fluctuating velocity.

Collision Frequency
A way of characterising inter-particle collision probability is the collision kernel, which is equal to two times the collision frequency over the particle number density (particles per unit volume) [43]. In the literature, there are two limiting cases for the behaviour of such collision kernel with particle Stokes number.
Saffman and Turner [44] provided an approximation for the collision kernel of drops in atmospheric turbulence. In their analysis, the particles are small compared to the smallest length scales of the fluid flow. Hence, particles follow the turbulent structures completely. The resulting collision kernel for mono-disperse particles is given by: Another limiting case was studied by Abrahamson [45], in which the motion of heavy particles suspended in high-intensity turbulence neglecting external forces is analysed. According to the kinetic theory, the particle motion is completely uncorrelated concerning the underlying fluid motion. Under the assumption that all components of the particle fluctuating velocity are isotropic, the resulting collision kernel reads: where k p represents the average particle turbulent energy of all tracked particles.
In practice, particles have finite inertia and, therefore, their collision kernel is expected to lie between both limiting values. Figure 11 presents the non-dimensional collision kernel obtained in this work as a function of the particle Stokes number for the three considered particle volume fractions. In addition, the collision kernels of Saffman and Turner, Γ ST , as well as Abrahamson, Γ A , are also shown in Figure 11 for comparison.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 22 [43]. In the literature, there are two limiting cases for the behaviour of such collision kernel with particle Stokes number. Saffman and Turner [44] provided an approximation for the collision kernel of drops in atmospheric turbulence. In their analysis, the particles are small compared to the smallest length scales of the fluid flow. Hence, particles follow the turbulent structures completely. The resulting collision kernel for mono-disperse particles is given by: Another limiting case was studied by Abrahamson [45], in which the motion of heavy particles suspended in high-intensity turbulence neglecting external forces is analysed. According to the kinetic theory, the particle motion is completely uncorrelated concerning the underlying fluid motion. Under the assumption that all components of the particle fluctuating velocity are isotropic, the resulting collision kernel reads: where represents the average particle turbulent energy of all tracked particles. In practice, particles have finite inertia and, therefore, their collision kernel is expected to lie between both limiting values. Figure 11 presents the non-dimensional collision kernel obtained in this work as a function of the particle Stokes number for the three considered particle volume fractions. In addition, the collision kernels of Saffman and Turner, Γ , as well as Abrahamson, Γ , are also shown in Figure 11 for comparison. Figure 11. Collision kernel, normalised by the Kolmogorov velocity, and the square of particle diameter, obtained by the present simulations as a function of the particle Stokes number and for the three average particle volume fractions. The Saffman and Turner [44] and Abrahamson [45] kernels for colliding particles are also shown for comparison.
Let us observe that the non-dimensional Saffman and Turner kernel is actually constant but that of Abrahamson depends on particle fluctuating velocity, which varies with particle volume fraction. However, for the simulated cases, is similar for the three values of . As it can be readily seen in Figure 11, the computed collision kernel is between the Saffman and Turner (lower limit) and the Abrahamson kernels (upper limit), as it had to be expected. In the case of zero particle inertia, represented by the Saffman and Turner limit, particle collisions are only controlled by the smallest fluid scales. This regime is approached in the present simulations of the lowest inertia particles, where the inter-particle collision time is determined by the Kolmogorov time scale. On the other limit, the collision kernel is governed by the particle kinetic energy, which is Figure 11. Collision kernel, normalised by the Kolmogorov velocity, and the square of particle diameter, obtained by the present simulations as a function of the particle Stokes number and for the three average particle volume fractions. The Saffman and Turner [44] and Abrahamson [45] kernels for colliding particles are also shown for comparison.
Let us observe that the non-dimensional Saffman and Turner kernel is actually constant but that of Abrahamson depends on particle fluctuating velocity, which varies with particle volume fraction. However, for the simulated cases, k p is similar for the three values of α p .
As it can be readily seen in Figure 11, the computed collision kernel is between the Saffman and Turner (lower limit) and the Abrahamson kernels (upper limit), as it had to be expected. In the case of zero particle inertia, represented by the Saffman and Turner limit, particle collisions are only controlled by the smallest fluid scales. This regime is approached in the present simulations of the lowest inertia particles, where the inter-particle collision time is determined by the Kolmogorov time scale. On the other limit, the collision kernel is governed by the particle kinetic energy, which is approached in the present simulations for the highest considered Stokes number. Eventually, for larger inertia particles, the collision kernel will decrease because they tend to be less responsive to fluid turbulence. The trend of the computed kernel in present simulations agrees qualitatively with the results of Vosskuhle [46], showing a fast increase around Stokes numbers close to one and then growing much more moderately for higher Stokes numbers (here up to 10).
Regarding the dependence of the simulated non-dimensional collision kernel with particle volume fraction, Figure 11 shows that it does not have an appreciable effect: only a marginal increase of Γ/ u K d 2 p is observed with growing α p .

Two-Particle Velocity Correlation Function
As it was previously mentioned, it is necessary to account for the correlation of the motion of neighbouring particles and their preferential concentration when dealing with the description of inter-particle interactions in practical two-phase flow computations. As commented in the introduction section, in this context, Lagrangian-based stochastic collision models offer an alternative to the costly, and often impracticable, deterministic collision algorithms. In such stochastic models, along the trajectory of each particle, a fictitious collision partner is generated at each time step. The properties of the fictitious particle such as diameter or velocities are sampled randomly from previously constructed local distributions, and the probability of collision is computed according to kinetic theory. However, in this process, the fluctuating velocity correlation of both particles has to be taken into account [29,32] and such correlation depends on the particle response to the turbulent environment. In the original paper of Sommerfeld [31], it was modelled by the following decaying exponential: where St L is the Stokes number based on the fluid Lagrangian integral time scale. Coefficients appearing in Equation (6) were obtained after comparison with LES computations of Reference [34] by trial and error. Such correlation has been employed successfully for predicting the effect of inter-particle collisions in pneumatic conveying and particle separation processes [9][10][11]35], and it is shown in Figure 12. A more rigorous way to determine the particle-pair velocity correlation function between colliding particles, R(St K ), was developed in Lain et al. [36]. In such study, the model proposed in Zaichik and Alipchenkov [37] (or its refined version, [18]) developed for homogeneous isotropic turbulence was solved to build the Eulerian particle-pair velocity correlation function between colliding particles, R E (St K ). The conclusions obtained in Reference [36] are revised in this section in light of the new DNS results obtained. The information relevant to the estimation of the correlation function of particle-pair fluctuating velocities can be extracted from the solution of the Zaichik and Alipchenkov model. To this end, they define the Eulerian two-particle velocity correlation as: where k p is the kinetic energy of particles, S pij is the particle structure function of the relative velocity between two particles and r is the distance between the two particles [37]. Expression (7) is obtained assuming that there is a balance of net radial inward and outward fluxes implying a zero mean radial relative velocity. The longitudinal two-particle velocity correlation function is defined as R(r) = 3B pll (r)/2k p in terms of the longitudinal space correlation, B pll , which depends on the separation distance and Stokes number. As pointed out by Zaichik and Alipchenkov [37] and discussed in depth in Février et al. [19], in contrast to the well-known result for fluid turbulence, this correlation for inertial particles is not tending towards unity in the limit of r → 0 . This specific feature of the spatial correlation indicates that as a consequence of the inertial bias in their trajectories, the velocities of particles are not completely correlated, even at r = 0. In the case of colliding particles, the distance between them tends to zero (obviously it cannot be less than the sum of the radii of both particles), so the colliding particle-pair velocity correlation could be expected to be approximated by the behaviour near the origin of the longitudinal space correlation, R E = 3B pll (r → 0)/2k p , which clearly depends on particle inertia. Based on the results obtained by solving the Zaichik and Alipchenkov model, Lain et al. [36] proposed the following functional form for the particle velocity correlation of colliding particles: Equation (8) uses three parameters instead of the two originally proposed by Sommerfeld [30], Equation (6). The parameters take the values: α = 0.019, β = 1.725 and γ = 0.044.
Correlation (8) is shown in Figure 12 as a solid line together with our own results obtained solving the Zaichik and Alipchenkov model, represented by plus symbols. As it can be seen, R E (St K ) tends to one in the limit of very small Stokes numbers according to the Saffman and Turner limit, where velocities of neighbouring particles are spatially correlated due to the interaction with the same local fluid flow. As St K increases, particle velocity correlation decreases smoothly, and in the limit of very high Stokes numbers, it tends to zero, where statistics of heavy particles satisfy the assumption of molecular chaos according to Abrahamson [45].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 3 ( → 0)/2 , which clearly depends on particle inertia. Based on the results obtained by solving the Zaichik and Alipchenkov model, Lain et al. [36] proposed the following functional form for the particle velocity correlation of colliding particles: Equation (8) uses three parameters instead of the two originally proposed by Sommerfeld [30], Equation (6). The parameters take the values: α = 0.019, β = 1.725 and γ = 0.044.
Correlation (8) is shown in Figure 12 as a solid line together with our own results obtained solving the Zaichik and Alipchenkov model, represented by plus symbols. As it can be seen, ( ) tends to one in the limit of very small Stokes numbers according to the Saffman and Turner limit, where velocities of neighbouring particles are spatially correlated due to the interaction with the same local fluid flow. As increases, particle velocity correlation decreases smoothly, and in the limit of very high Stokes numbers, it tends to zero, where statistics of heavy particles satisfy the assumption of molecular chaos according to Abrahamson [45].  [18]. Cross (×) symbols show the results of Reference [28]. Correlation of Reference [30], Equation (6), corresponds to the purple dashed line. The rest of the symbols refer to the LBM-based DNS results obtained in this work for different average volume fractions. Closed symbols represent the Eulerian correlation evaluated at a distance = 3 /2, Equation (9), and open symbols represent the obtained colliding pair velocity correlation, Equation (10) (see text for details). The dot-dashed lines fit the colliding two-particle velocity correlation functions, ( ), of the different volume fraction cases studied in this work.
Results for the particle pair correlation function obtained in the present work within the frame of LBM-based DNS are also presented in Figure 12. The closed symbols correspond to the Eulerian two-particle velocity correlation function evaluated at the fixed distance = 3 /2, i.e., This corresponds to a situation where the longitudinal velocity correlation between neighbouring particles is evaluated at 1.5 times the particle diameter (i.e., around 0.3 ) at a fixed time. Correlation has been computed at different times after the initial transient has been completed, showing stable values in all the range of Stokes numbers studied.  [18]. Cross (×) symbols show the results of Reference [28]. Correlation of Reference [30], Equation (6), corresponds to the purple dashed line. The rest of the symbols refer to the LBM-based DNS results obtained in this work for different average volume fractions. Closed symbols represent the Eulerian correlation evaluated at a distance r = 3d p /2, Equation (9), and open symbols represent the obtained colliding pair velocity correlation, Equation (10) (see text for details). The dot-dashed lines fit the colliding two-particle velocity correlation functions, R c (St K ), of the different volume fraction cases studied in this work.
Results for the particle pair correlation function obtained in the present work within the frame of LBM-based DNS are also presented in Figure 12. The closed symbols correspond to the Eulerian two-particle velocity correlation function evaluated at the fixed distance r = 3d p /2, i.e., This corresponds to a situation where the longitudinal velocity correlation between neighbouring particles is evaluated at 1.5 times the particle diameter (i.e., around 0.3 η K ) at a fixed time. Correlation has been computed at different times after the initial transient has been completed, showing stable values in all the range of Stokes numbers studied.
On the other hand, open symbols represent the two-particle velocity correlation for particle pairs that have actually collided, i.e., This corresponds to a situation where the longitudinal velocity correlation between neighbouring particles is computed for particle pairs that have experienced a collision along the development of the simulation. Statistics are computed based on 10 6 collision events.
From Figure 12, it can be clearly seen that the R E values obtained from LBM-based DNS agree reasonably well with the correlation Equation (8). Therefore, the theoretical model of Zaichik and Alipchenkov is adequate for describing the Eulerian velocity correlation between neighbouring particles, as it could be expected. Moreover, the DNS results for different solid volume fractions show nearly no scattering of values for all the range of studied Stokes numbers.
As already commented in the introduction, the decreasing of particle velocity correlation of neighbouring particles, i.e., R E , with increasing particle inertia is usually explained [19] in terms of an increase of the relative velocity between them, consistent with the existence of caustics and a component of random uncorrelated motion of particle velocity field. In fact, the connection between the decreasing of particle velocity correlation and the increase of the modulus of relative velocity (or vice versa) can be easily cast in mathematical terms. The average of the modulus of relative velocity between two particles is computed as: Now, because the mean particle velocity in isotropic homogeneous turbulence is zero, where k p is the fluctuating kinetic energy of the particles. Finally, solving for the two-particle velocity correlation, it is obtained: where u p1 and u p2 are the velocities of particles 1 and 2 respectively, w = u p2 − u p1 is their relative velocity and denotes the averaging operator. As Equation (11) is a general relation, a reduction in the two-particle velocity correlation is accompanied by an increase in the modulus of relative velocity, and an increase of this quantity implies a decrease of the correlation. To support this affirmation, Figure 13 shows the modulus of radial relative velocity (RRV), defined as w r = w·r/|r|, between colliding particles, divided by the Kolmogorov velocity, u K , as a function of the particle Stokes number. It is clear that it increases with particle inertia, experiencing a rapid growth for Stokes numbers roughly between 0.5 and 2. The strong increase of the radial relative velocity for St K > 0.5 is connected with the strong enhancement of the inter-particle collision rate due to sling caustics occurrence. For larger Stokes numbers, of more than 10, a decrease of the RRV is expected because the particles become slower and less responsive to the fluid turbulence. Such results are in close agreement with those of Figure 2 in Reference [18].
The behaviour of the colliding particles velocity correlation function, R c , in Figure 12 deserves attention. At small values of particle inertia, here St K = 0.1, there is no difference with the Eulerian correlation as the collision events are governed by the fluid turbulent shear. As soon as inertia increases, R c adopts significantly lower values than R E ; for instance, in Figure 12, at St K = 0.5, the difference between the two correlations is already appreciable. Such differences are increasing with augmenting particle Stokes number. The reason for the smaller values of R c , regarding R E , is that in three-dimensions, not all particle pairs that are close enough to interact do really collide, it is necessary that they approximate each other with the right radial relative velocity to have the interaction. In fact, the R c curves cannot be obtained as R E r = d p because the RRV's are different, as already remarked in Reference [16]. The reason is that in the Eulerian case, the average of the RRV is zero, as it was implied in the derivation of Equation (7), while in the colliding particle case, this RRV is effectively negative. As a result, the modulus of the particle relative velocity is higher in the colliding case than in the Eulerian case which, due to Equation (10), has the consequence that the values of R c are smaller than those of R E . As the modulus of the RRV increases with Stokes number (Figure 13), the distance between both curves augments with particle inertia. Moreover, the results of the recent study of van Wachem et al. [29] for R c are also depicted in Figure 12 as crosses, showing a good agreement with those obtained in the present simulations. However, in Reference [29], the volume fraction values were different for the four particle types analysed and ranged from 0.3% to 2%. values of are smaller than those of . As the modulus of the RRV increases with Stokes number (Figure 13), the distance between both curves augments with particle inertia. Moreover, the results of the recent study of van Wachem et al. [29] for are also depicted in Figure 12 as crosses, showing a good agreement with those obtained in the present simulations. However, in Reference [29], the volume fraction values were different for the four particle types analysed and ranged from 0.3% to 2%. Another interesting fact that can be observed in the correlation functions in Figure 12 is that, for a fixed particle Stokes number, their values increase with increasing solid volume fraction. The counterpart, seen in Figure 13, is that the radial relative velocity between colliding particles decreases with increasing volume fraction. This fact was already shown by the authors of Reference [38], who also observed that the mean angle formed by the velocities of colliding particles was increasing with particle inertia. These authors ascribed such behaviour to a kind of collective motion of particles in clusters. Such phenomenon has also been found in the present work and it was illustrated in Figure 8a. However, in the present study, such trend is not as evident as in that research because they employed higher values of than in this work. In summary, it has been found that the colliding particle pair velocity correlation is smaller than that computed by the Eulerian particle pair velocity correlation as particle separation tends to zero. The practical consequence for the stochastic collision models in the Euler-Lagrange frame is that the velocity correlation between the real and fictitious particles cannot be drawn out from the Eulerian two-particle velocity correlation given by the model of Zaichik and Alipchenkov, but instead should be built from the colliding particle-pair correlation function, which depends on both particle inertia and volume fraction. As an example (which does not pretend to be universally valid), functional forms for such correlation can be devised from the presented DNS results and Table 3 shows the fitted coefficients for the curves to Equation (8) for the volume fraction cases studied in this work. The corresponding curves of such fits are also plotted in Figure 12 for comparison, which could be used for building the fluctuating component of fictitious collision partners in stochastic models of inter-particle collisions, e.g., Equation (23) in Reference [31]. Nevertheless, it would be necessary to study the dependence of such correlations on turbulent Reynolds number and in a larger range of particle volume fractions. Figure 13. Modulus of mean radial relative velocity (RRV) between particle pairs, w r = w·r/|r|, as a function of the particle Stokes number normalised by the Kolmogorov velocity, u K . Half-closed symbols represent colliding particles while open symbols denote the average of close particle pairs at a distance lower than r = 3d p /2 based on Eulerian evaluations of the radial distribution.
Another interesting fact that can be observed in the R c correlation functions in Figure 12 is that, for a fixed particle Stokes number, their values increase with increasing solid volume fraction. The counterpart, seen in Figure 13, is that the radial relative velocity between colliding particles decreases with increasing volume fraction. This fact was already shown by the authors of Reference [38], who also observed that the mean angle formed by the velocities of colliding particles was increasing with particle inertia. These authors ascribed such behaviour to a kind of collective motion of particles in clusters. Such phenomenon has also been found in the present work and it was illustrated in Figure 8a. However, in the present study, such trend is not as evident as in that research because they employed higher values of α p than in this work.
In summary, it has been found that the colliding particle pair velocity correlation is smaller than that computed by the Eulerian particle pair velocity correlation as particle separation tends to zero. The practical consequence for the stochastic collision models in the Euler-Lagrange frame is that the velocity correlation between the real and fictitious particles cannot be drawn out from the Eulerian two-particle velocity correlation given by the model of Zaichik and Alipchenkov, but instead should be built from the colliding particle-pair correlation function, which depends on both particle inertia and volume fraction. As an example (which does not pretend to be universally valid), functional forms for such correlation can be devised from the presented DNS results and Table 3 shows the fitted coefficients for the R c curves to Equation (8) for the volume fraction cases studied in this work. The corresponding curves of such fits are also plotted in Figure 12 for comparison, which could be used for building the fluctuating component of fictitious collision partners in stochastic models of inter-particle collisions, e.g., Equation (23) in Reference [31]. Nevertheless, it would be necessary to study the dependence of such correlations on turbulent Reynolds number and in a larger range of particle volume fractions.

Conclusions
In this paper, the effect of particle inertia and void fraction, and therefore inter-particle collisions, on the colliding particle-pair velocity correlation in homogeneous isotropic turbulence has been addressed by using direct numerical simulations based on the Lattice Boltzmann Method. Three particle volume fractions have been considered, α p = 0.01%, 0.05%, 0.1%, in the dilute regime.
In the range of void fractions studied, collisions have a marginal impact on particle transport, with the aerodynamic interaction with fluid turbulence being the main mechanism that drives it. This is reflected in the fact that collision Stokes number has been always larger than one, keeping the dilute flow conditions. A main result of the performed study is that collision locations of particles of intermediate inertia, St K ∼ 1, happens in areas where the fluid strain rate and dissipation are higher than the corresponding averaged values at particle positions. Connected with this fact, the average kinetic energy of colliding particles of St K ∼ 1 is lower than the value averaged over all particles. This finding was explained by the fact that the fluid kinetic energy at collision positions is lower than the mean value and, because St col 1, particles have enough time to accommodate to the fluid environment between consecutive collisions. In case of the highest inertia, St K = 10, particles do not experience noticeable clustering in areas of high strain rate and in that case, inter-particle collisions are favoured by high kinetic energy; as a result, the colliding particle kinetic energy is higher than the corresponding averaged value over all particles.
Later, the Eulerian and colliding two-particle velocity correlation functions were computed in the same homogeneous isotropic turbulence configuration. As a result, the Eulerian function was well reproduced by the Zaichik and Alipchenkov [37] model, but the colliding particles correlation showed remarkably smaller values due to the role played by the mean radial relative velocity. In particular, it was demonstrated that the colliding pair correlation cannot be obtained as a limiting value of the Eulerian correlation. Moreover, it has been confirmed that the colliding particles velocity correlation increases with particle volume fraction, a fact related with the decreasing of the average relative velocity between the colliding particles and attributed to the particle dynamics inside clusters. However, the decreasing of radial relative velocity with increasing α p was also found for large particle inertia, where no appreciable clustering was observed. Therefore, further studies are needed to confirm if the same phenomenon happens for a wider set of conditions than those presented here, such as larger volume fractions and turbulent Reynolds numbers.
The obtained results in the present study can be used to support inter-particle collision stochastic models (e.g., Reference [31]) for use in general-purpose two-phase flow Euler-Lagrange codes. However, it is necessary to study the dependence of such correlations on turbulent Reynolds number and in a larger range of particle volume fractions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In the present study, the modelling and detection of particle-particle collisions was computed using a deterministic collision model proposed in Reference [47]. This model assumes that the collisions are binary and quasi-instantaneous. Furthermore, the elastic contact between two hard spheres occurs in a single point. In order to determine whether a collision occurs, the time-sequenced collision detection algorithm described in detail in Reference [38] was applied. Such algorithm is briefly summarised here.
Firstly, all particles are tracked to the end of the present Lagrangian particle time step, ∆t. Inter-particle collisions have occurred if particles overlap after the actual time interval. In the case of a collision, the distance between their centers is equal to or less than the half sum of the two particle diameters. The determination of all existing collision pairs is realised under application of the cell-index method [47]. The potential collision partners of a given particle are found among the other particles in the cell plus the particles in all 26 neighbouring cells. The list of particles in each cell is set up and updated using linked lists. Assuming that the particles move at constant velocities during the particle time step, the collision time for each overlapping pair, ∆t ij,c , is calculated by solving the following quadratic equation: x pij (t + ∆t) + u pij (t + ∆t)∆t ij,c = 0.5 d pi + d pj (A1) Here, x pij and u pij are the relative position and velocity of particles i and j which have diameters d pi and d pj , respectively. Equation (A1) suffers from the limitation that it will not detect a particle-particle collision if complete interpenetration occurs during the given particle time step. Subsequently, the determined collision times (i.e., the time between the beginning of the actual time step and the time of particle collisions) are sorted in descending order. In the time-sequenced implementation, the first overlapping pair is moved backward in time to the position of collision. The determination of the post-collision velocities of the collision couple is based on the laws of impact for hard spheres provided in Reference [48]. A detailed derivation of the used collision operators can be found in Reference [38].
After the collision is computed, both particles are again moved forward in time to complete the actual particle tracking time step. New overlaps may be added to the list of colliding particles, i.e., in the case of multiple particle collisions within one tracking time step, while processed or non-existent ones are deleted. These collision computations are carried out until no overlaps are left.
The computational cost of the searching collision partners can be estimated as follows. In the employed Cell Index Method, the simulation box is divided in a regular array of M 3 cells; on average, there are N p /M 3 particles (N p is the total number of particles) within each cell, where 27 neighbouring cells have to be searched per particle. Therefore, avoiding double counting, the search algorithm requires 27 2 N p /M 3 operations [47].