Reynolds number scaling of particle clustering in turbulent aerosols

Particles with finite inertia in turbulent flow cluster in low-vorticity regions of the fluid due to the inertial imbalance between the denser particles and the lighter surrounding fluid. This effect, sometimes referred to as preferential concentration, has been investigated in several recent numerical studies. Sundaram and Collins (1997 J. Fluid Mech. 335 75) considered the effect of particle clustering on the interparticle collision rate and showed that the radial distribution function, evaluated at contact, precisely corrects the collision kernel for this effect. An open question is how preferential concentration scales with Reynolds number. We investigate this question using direct numerical simulations (DNS). Over the limited range accessible by DNS, it appears that the radial distribution function approaches a plateau with increasing Reynolds number. This contradicts earlier studies that predicted linear growth with Reynolds number. The implications of these findings for very high Reynolds number applications such as cloud formation in the upper atmosphere is briefly discussed.


Introduction
Aerosol particles with finite inertia are known to cluster in low-vorticity regions of a turbulent flow. Vortices act as 'centrifuges,' forcing particles out of vortex cores and into the strain fields that lie between the vortices. Figure 1 shows a snapshot of a three-dimensional direct numerical simulation (DNS) of particles in isotropic box turbulence. The green surfaces are isocontours of the magnitude of vorticity evaluated at a relatively high value (six times the standard deviation of the fluctuating vorticity) and the white surfaces are isocontours of the particle concentration field evaluated at 10 times the mean particle concentration. The tube-like structures of the vorticity field are a well-known feature of isotropic turbulence (see e.g. [1]); the remarkably similar geometrical shape of the particle concentration isocontours highlights the strong connection between particle clustering and fluid vorticity.
Particle inertia is required to allow particle trajectories to deviate from fluid trajectories, which are incompressible in this study. The most important measure of particle inertia is the Stokes number, defined as where τ p ≡ βσ 2 /18ν is the particle relaxation time, σ the particle diameter, β the particle-to-fluid density ratio, ν the fluid kinematic viscosity, τ η ≡ √ ν/ the Kolmogorov time scale, η ≡ ν 3/4 / 1/4 the Kolmogorov length scale and the turbulent energy dissipation rate. We normalize the particle parameters by the Kolmogorov scales (small scales) because these are the scales that are principally responsible for clustering. Based on this normalization the degree of clustering peaks when St ≈ 1 [2].
Particle clustering has a profound effect on a number of aerosol processes (e.g. settling, dispersion, turbulence modulation). The focus of this paper will be on the effect particle clustering has on the interparticle collisions. Collision is a binary process with a rate that is proportional to the local number density squared; consequently, fluctuations in the particle concentration field due to clustering will enhance the overall collision rate. This effect was analysed by Sundaram and Collins [3], who introduced the radial distribution function (RDF) into the collision kernel to account for particle clustering. The resulting collision kernel for a monodisperse population Three-dimensional isocontour of the vorticity magnitude evaluated at six times the standard deviation of the vorticity fluctuation (green surfaces) and isocontour of the particle concentration field evaluated at 10 times the mean concentration (white surfaces) taken from a 64 3 DNS. Notice the localized pockets of high particle concentration.
of particles of diameter σ, in the absence of Brownian diffusion or other molecular effects (i.e. infinite Peclet number) can be written as [4] where g(r) is the RDF, w r the radial component of the relative velocity between two particles and P(w r |r) the probability density function (PDF) of w r conditioned on the particle pair separation distance, r. The correction to the collision kernel for particle clustering is completely captured by the factor g(σ), which is normalized to be unity for a uniform distribution of particles and increases with increasing clustering, reaching values as large as 10-100 for St ≈ 1.0. An effect of this magnitude can profoundly alter the size distribution of an aerosol undergoing coalescence (liquid) or coagulation/agglomeration (solid). Indeed, this was observed in the DNS study by Reade and Collins [5], who found that particle clustering increased the growth rate of the particle population and significantly broadened the particle size distribution when compared with theories that do not take clustering into account.
Recently, the clustering phenomenon was postulated to explain a decades-old question in the cloud physics literature. Standard cloud models assume that droplets below 10 µm have low coalescence efficiencies and, therefore, grow principally by condensation alone. However, two observations are inconsistent with this assumption: (i) cloud droplet distributions are much broader than can be predicted by surface growth models (note that models of condensational growth are reasonably accurate and so the discrepancy is not likely due to modelling error); and (ii) initiation of warm rain occurs in much shorter times than can be predicted by these models. Bridging the so-called 'size gap' between the nucleated drops and coalescing drops remains an open question in the cloud physics literature [6,7]. Falkovich et al [8] hypothesized that cloud droplet clustering may be at least partially responsible for the accelerated growth rate. To put it succinctly, droplet clustering may compensate for the lower collision efficiency of the smaller droplets, allowing them to coalesce at a much earlier stage of their growth history. Kostinski and Shaw [9] found evidence of cloud droplet clustering in experimental measurements of droplet distributions made by Chaumat and Brenguier [10]. An alternative argument put forth by Shaw et al [11] is that droplet clustering due to turbulence will lead to fluctuations in the vapour supersaturation field. This could then allow droplets to experience an enriched vapour field, if for example a vortex that existed previously had removed all droplets in a local region of the fluid. Model simulations of this scenario yielded broader drop size distributions; however, the broadening was mainly in the direction of smaller drop sizes, due to the increase in secondary nucleation events in the vapour-enriched regions of the fluid. This explains the broader size distributions, but does not explain the accelerated growth rate observed in clouds.
While both arguments are appealing, there are issues concerning the parameter values in the atmosphere that must be addressed. Dimensional analysis shows that the RDF for a monodisperse population of particles depends on the particle Stokes number, the size parameter σ/η, the volume fraction of particles, φ and the turbulence Reynolds number (based on the Taylor microscale), R λ [12]. Typical values for these quantities in clouds are given in table 1 4 along with the range of values for the DNS 5 and controlled experiments. 6 Under the dilute conditions found in clouds, the dependence on φ is known to be weak and can be neglected [16]. Numerical simulations have shown that clustering peaks around St ≈ 1 and falls off rapidly on either side of this optimum value [12]. Conditions in the cloud correspond to St < 10 −2 for the size range of interest (i.e. micron-size droplets); thus, the degree of clustering should be weak, unless there is compensation by the much larger value of the Reynolds number found in the atmosphere. It is not possible to reach the atmospheric Reynolds number in a DNS; thus, we are forced to consider the scaling of clustering with Reynolds number over the range achievable in the DNS, and then extrapolate to atmospheric conditions. 4 It should be noted that fluid mechanical conditions in a cloud are difficult to measure and are not that well known. Thus, these estimates should be viewed with some skepticism. For more details on this issue see the discussion of experimental observations in [7] (section 7). 5 The largest DNS of a Newtonian fluid (without particles) to date by Kaneda et al [13] used a 4096 3 mesh and reached R λ = 1201. 6 See [14,15] for typical experimental conditions. It is useful to first consider what input theory might provide to this question. According to the classical Kolmogorov [17] argument, increasing the Reynolds number (while fixing the small scales) involves increasing the integral scale of the turbulence. As vorticity is mainly contained within the small scales, changes to the large scales will be inconsequential as long as the Reynolds number is large (a more precise statement is the integral scale will introduce a correction to the fluid vorticity that scales like (1 − R −2 λ ), where R λ is the Reynolds number based on the Taylor microscale). The logical conclusion, based on this argument, is that the effect of Reynolds number on the RDF should vanish in the limit R λ → ∞.
However, the more modern view of isotropic turbulence recognizes that the fluctuating signal, particularly for small-scale phenomena such as vorticity, is highly intermittent (see e.g. [18] (chapter 8) or [19] for a detailed discussion); moreover, the degree of intermittency increases with increasing Reynolds number. The argument in the previous paragraph does not take intermittency into account. One important consequence of intermittency is the PDF of vorticity has exponential tails that broaden with increasing Reynolds number. It seems reasonable to assume that these stronger vortical events (with increasing Reynolds number) may influence the RDF. One can speculate on the consequences of intermittency using the Kolmogorov [20] correction to his earlier theory. In the modified theory, is no longer fixed, but assumed to have a log-normal distribution, with a variance that scales like (see e.g. [21] (section 6.7)) (3) As the particle Stokes number is proportional to √ , we can consider the probability that a particle with a nominal Stokes number of 0.01 (typical value in a cloud) achieves a local value > 0.1 (corresponding to >100 times its nominal value). The answer for R λ = 100 (typical DNS value) is 0.0016%, while the answer for R λ = 10 000 (typical atmospheric value) is 0.033%. Thus, it is more than 20 times more likely to occur in the atmosphere than in the DNS. Moreover, this fraction increases indefinitely with increasing Reynolds number, and so it seems at least plausible, based on these arguments and observations, that the RDF may increase indefinitely as well.
Determining which of these two arguments is correct has important implications for the cloud physics problem and other high-Reynolds-number applications (see e.g. discussions of turbulence effects in protoplanetary nebula by Cuzzi and Hogan [22] and references therein). Numerical simulations to date have yielded mixed results, mainly due to the limited range of Reynolds number that can be explored. For example, Reade and Collins [12] observed a slight increasing trend with increasing Reynolds number, but the statistical error associated with the highest Reynolds number in that study made it impossible to draw a definitive conclusion. Wang et al [23] proposed a linear scaling of the RDF with R λ , while Hogan et al [24] correlated their particle concentration fluctuations using multi-fractals [19].
In this study, we reconsider this question using DNS over the Reynolds number range 57-152 (corresponding to 64 3 -192 3 points). To isolate the effect of the Reynolds number, we employ a forcing scheme that maintains the average dissipation rate nearly constant and so minimizes fluctuations in the Stokes number. Nevertheless, the statistical scatter in the RDF is still significant (relative to the trends we are hoping to identify); we use a confidence analysis to sort out the most likely trend with Reynolds number.

Overview of simulations
Isotropic simulations in a periodic cube of length 2π along each side have been performed. The incompressible fluid is governed by the continuity and Navier-Stokes equations, shown below: where u is the fluid velocity, t the time, ρ the fluid density, p the pressure, ν the fluid kinematic viscosity and F a large-scale forcing function added to maintain stationary turbulence. The equations are solved using a standard pseudospectral technique that is described elsewhere (see [25] for details). One unique element of the simulation is the forcing function, F. To reduce the magnitude of the fluctuations in the average dissipation rate, , which complicates the analysis of the results, we apply deterministic forcing over the wavenumber range 1 k √ 2 as shown belowF The factors f 1 and f √ 2 , defined below, are chosen to maintain the fluid energy and dissipation rates at preset target values, e T and T , respectively, where e ≡ e T − e * , ≡ T − * and the asterisk indicates the value after one iteration without forcing. E * (1) and E * ( √ 2) are the energy spectrum, at the indicated wavenumbers, after one iteration. The inequalities in equations (7) and (8) guarantee f 1 > 1 and f √ 2 > 1 so that the forcing cannot remove energy from the system. One consequence of this constraint is that the control of the dissipation rate is only from below, meaning excursions in above the target value cannot be controlled by this forcing scheme. Instead, forcing is turned off until the system returns to or falls below the setpoint.  The particles are simulated by following the trajectory of each particle in a Lagrangian frame of reference. The governing equations are as follows: where x i and v i are the position and velocity of the ith particle, respectively, and u(x i ) is the undisturbed fluid velocity at the point x i . Equation (10) neglects the added mass, Basset history and Faxen corrections that arise in the complete analysis of the forces acting on a particle in a timedependent flow field (see e.g. [26]). These are reasonable assumptions under the circumstances that the particle size is small (i.e. σ/η 1) and the density of the surrounding fluid is small (i.e. β 1), both of which apply to this study. In this regime, the fluid inertia associated with the local flow around a given sphere is negligible due to the small local Reynolds number, while particle inertia is still significant because of the extra factor of β in the definition of the Stokes number; β for most aerosol systems is O(1000).
We also neglect hydrodynamic forces resulting from the disturbance flow around each particle, and collisions. These effects are important for understanding the so-called collision efficiency, defined as the fraction of crossing trajectories that achieve molecular contact. This issue is a delicate one due to the competing hydrodynamic and molecular forces, all of which diverge as the particles come into contact. Consequently, the efficiency is sensitive to molecular parameters such as the Hamaker constant (van der Waals attractive force) and the Knudsen number (ratio of the mean free path to the particle dimension). The most comprehensive treatment to date has been done by Koch and co-workers [27,28]. In this study, we are interested in the collision rate in the absence of these considerations. Reade and Collins [12] showed that a reasonable approximation to the true RDF is obtained neglecting all of these effects-errors arise only for r ≈ σ. As we are concerned with σ η, there is a region σ r η where the non-interacting particle simulation is useful. We will restrict our attention to this region. This allows us to simulate a much larger population of particles and obtain accurate statistics below r/η ≈ 10 −2 . As will be shown, this degree of accuracy is essential for shedding light on the Reynolds number scaling.

Results and discussion
Numerical simulations of particles with Stokes numbers of 0.4, 0.7, 1.0 and 1.5 are performed at the four Reynolds numbers (R λ = 65, 81, 107 and 152). In each simulation, N = 2 million identical particles are initially placed at random locations in the fluid and assigned the fluid velocity. Simulations are performed for six eddy turnover times to allow the particles to equilibrate with the fluid and then statistics are taken over the next 8-10 eddy turnover times. The statistical measure of particle clustering that is most relevant for predicting the collision kernel is the RDF, g(r), defined as the ratio of the number of particle pairs found at a given separation distance to the expected number if the particles are uniformly distributed [29] ('uniformly distributed' means the probability of finding a given particle at a location in space is independent of the positions of the other particles). We compute this quantity by binning all of the particle pairs according to their separation distance and calculating [30] g(r i ) = where N p ≡ N(N − 1)/2 is the total number of pairs, V the total volume of the system (8π 3 ), N p (i) the number of particle pairs with separation distances that lie between r i − r/2 and r i + r/2 and V i ≡ 4 3 π[(r i + r/2) 3 − (r i − r/2) 3 ] the volume of the shell associated with the nominal separation distance r i . In this study, r/η = 0.005, based on statistical convergence considerations. Figure 3 shows the Reynolds number dependence of the pair correlation function, defined as h(r) ≡ g(r) − 1, at the four values of Stokes number. What is immediately apparent is the sensitivity of all of the curves to changes in the Reynolds number is relatively weak. A more careful examination of the curves shows that the sensitivity to Reynolds number is greater at the higher values of Stokes number than the lower values. Both of these observations are more consistent with the classical Kolmogorov scaling argument that predicts saturation with increasing Reynolds number. To better quantitfy the dependence on Reynolds number, we fit all of the curves at small values of r/η to a power-law expression of the form

Results for the radial distribution function
using a least-squares regression to obtain the parameters c 0 and c 1 . Analysis of variance (ANOVA) tables [31] are generated for each case (see [32] for a complete listing of the tables). From the ANOVA tables, we are able to say with greater than 99% certainty that a power law is the best approximation of h(r) for r/η 1. To try and discern, via statistical analysis, whether the regressed parameters c 0 and c 1 are converging or diverging with respect to Reynolds number, we performed a second linear least-squares fit of both parameters to R λ and 1/R λ . The argument is that we are comparing two potential Taylor series If we obtain a better fit with equation (13), the data support the indefinite growth argument (or at least are consistent with it), whereas if we obtain a better fit with equation (14), the data support the saturation model. It must be emphasized that our goal is not to determine the functional dependence of c 0 and c 1 on Reynolds number, but rather to distinguish between these two possibilities. Table 3 shows the lack of fit statistics for all of the linear regressions performed. Lack of fit (or F statistic) is the ratio of the error associated with a chosen model to the intrinsic error in the data [31]. A lower value of F implies the model is doing a better job (statistically) of representing the data. The baseline value refers to the fit that allows neither c 0 nor c 1 to depend on Reynolds number. The remaining columns show the lack of fit statistics, allowing c 0 then c 1 to depend on R λ and 1/R λ , respectively. With the exception of St = 0.4, the fit to 1/R λ is superior in each case, supporting the saturation argument. The St = 0.4 case is strongly affected by the value of c 0 and c 1 obtained at R λ = 152, which is clearly high; however, this point also has the greatest statistical error of all of the curve fits performed. If you remove this point, the trend for St = 0.4 is in line with the others.
Our results contradict the predicted scaling of Wang et al [23], who proposed linear growth of the RDF with Reynolds number. However, they simulated particles that are relatively large (σ ≈ η), and as such did not consider the power-law region at small separations. Consequently, their results more closely correspond to the pre-exponential factor c 0 , for which we too see a stronger Reynolds number dependence. Nevertheless, our results do not support their predicted linear scaling. Table 3. Lack of fit (F statistic) for the linear regression of c 0 and c 1 with Reynolds number. Base refers to the baseline fit with no Reynolds number dependence, c 0 (R λ ) is the case with R λ dependence for the intercepts, c 1 (R λ ) is the case with R λ dependence for the slopes, c 0 (R −1 λ ) is the case with 1/R λ dependence for the intercepts and c 1 (R −1 λ ) is the case with 1/R λ dependence for the slopes.

St
Base

Theoretical considerations
In a recent study, Chun et al [33] developed a theoretical expression for c 1 for the limiting case St 1. The approach they used is similar to the ones taken in two recent papers [34,35]. Chun et al [33] made extensive comparisons of the theory with DNS and other stochastic simulations and found quantitative agreement. A brief review of the results relevant to this study is given below.
The theory considers particles in a Lagrangian frame of reference, following a particle with finite Stokes number. Putting oneself in the frame of reference of, say, the ith particle, we can derive the governing equations for the relative positionr ≡ x j − x i and relative velocity w ≡ v j − v i from equations (9) and (10): where ≡ ∇u is the undisturbed fluid velocity gradient tensor along the trajectory of the ith particle. In equation (16), we make the assumption thatr < η and, hence, the fluid velocity is a linear function of position based on the velocity gradient (t). If we further assume St 1, we can perform a perturbation expansion ofr and w as follows: r =r [0] + Str [1] Substituting these expressions into equation (16) and equating terms of equal order in St yields w [1] = ·r [1] − τ η d dt  [33]. Colours indicate the value of Reynolds number. the dependence of the pre-exponential factor and the power on Reynolds number. Within the range of Reynolds number, we could achieve (65-152); there appears to be saturation of both coefficients with increasing Reynolds number. There is a somewhat stronger dependence of the pre-exponential factor than the power, but overall both coefficients appear to approach a finite value in the limit R λ → ∞. The results are supported by a new theory for clustering that predicts the exponent in the power-law expression to be proportional to [ S 2 p − R 2 p ], which in our simulations is bounded since R 2 p > 0 and S 2 p is observed to be a decreasing function of R λ . The results have important implications for cloud physics. In particular, they suggest clustering is not strongly enhanced by the much larger Reynolds number in the atmosphere. Droplet Stokes numbers must be reasonably large (i.e. >0.1) to feel the effects of clustering; hence the argument of Falkovich et al [8] applies to droplets with diameters 10 µm. While clustering will surely accelerate the coalescence rates of droplets in this size range (or larger), the question of how these droplets are generated from sub-micron nucleates still appears to be open.