This paper examines the combined influences of turbulence and gravity on droplet collision statistics in turbulent clouds by means of direct numerical simulation (DNS). The essential microphysical mechanisms that determine the geometric collision kernel are explored by studying how gravity affects droplet relative velocities and preferential concentration of both monodisperse and bidisperse droplet distributions. To this end, collision statistics of large amounts of droplets with radii ranging from 10 to 90 μm, driven by a turbulent flow field and gravity, are calculated. The flow is homogeneous and isotropic and has a dissipation rate of ε = 4.25 × 10−2 m2 s−3. The results show that in the calculation of collision statistics, the interplay between gravity and turbulence is an essential element and not merely an addition of separate phenomena. For example, the presence of gravity leads to clustering of large droplets interacting with the larger scales of turbulence in the DNS. The collision statistics of a bidisperse droplet distribution, even with a very small radius difference, shows profoundly different behavior than the monodisperse case.
A classical challenge in cloud physics is the adequate description of microphysical processes that can explain the fast initiation of warm rain. Specifically, the fast growth of cloud droplets across the size gap of 10 to 50 μm in radius is not well understood (e.g., Wang et al. 2006b). In the past, various ideas have been put forward to explain the discrepancies between theory and observations. For example, entrainment and mixing of dry air by large-scale turbulence (e.g., Blyth 1993), the existence of ultragiant cloud condensation nuclei (e.g., Yin et al. 2000), and nonuniformity of the supersaturation field (e.g., Shaw et al. 1998) may lead to enhanced droplet growth in the condensational phase of rain formation. The general hypothesis, however, is that the interaction between in-cloud turbulence and droplet dynamics can have a pronounced effect on the collision probability. As a result, turbulence may significantly enhance droplet growth in the collision–coalescence stage of warm rain formation.
Reviews on the interactions between turbulence and droplets in clouds (Jonas 1996; Pinsky and Khain 1997; Vaillancourt and Yau 2000; Shaw 2003; Khain et al. 2007) mention at least three mechanisms that significantly affect the collision and coalescence probability of cloud droplets. First, relative velocities between colliding droplets are increased because of the continuously varying drag forces that droplets experience, leading to large variations in droplet velocities. Moreover, relative velocities of sedimenting droplets can be altered by turbulence through an increase in average settling velocities, a phenomenon known as preferential sweeping (Wang and Maxey 1993; Yang and Lei 1998; Dávila and Hunt 2001). Second, at the fine scales of turbulence, droplets are centrifuged out of vortical structures through their finite inertia and hence accumulate in regions of the flow where there is low vortical activity. This local clustering of droplets is known as preferential concentration (Squires and Eaton 1991; Reade and Collins 2000). Third, the gravitational description of the collision efficiency (see, e.g., Wang et al. 2005) needs to be modified in the presence of a turbulent flow because of the more complex behavior of local hydrodynamic droplet–droplet interactions (see also Wang et al. 2007; Pinsky et al. 2007).
Over the last decade, increased computational power has made it possible to study these microphysical processes in more detail through the use of direct numerical simulation (DNS). The major advantage of DNS is that it explicitly solves the smallest turbulent scales, where most interaction with droplet dynamics is known to take place (e.g., Shaw 2003), considering the small droplet sizes. Franklin et al. (2005, 2007) employed DNS in the context of cloud microphysics and found turbulence-induced increases of the collision probability of droplets to be dependent on the dissipation rate. They also showed that these increases originate from preferential concentration and increased relative velocities. Corresponding results were obtained by Wang et al. (2005, 2007), who also included hydrodynamic effects in their DNS to study its effect on droplet clustering, relative velocities, and settling velocities. In a theoretical study, Ghosh et al. (2005) explored the effect of increased settling velocities in the context of fast rain initiation.
Apart from collision efficiencies due to droplet–droplet hydrodynamic interactions, the collision–coalescence process is to a large extent determined by the geometric collision kernel Γ12, which measures the collision probability from the geometric overlap between droplets. In the past, various attempts have been undertaken to model the geometric collision kernel in a turbulent flow. Saffman and Turner (1956) developed an expression for weakly inertial particles in turbulence in the presence of gravity, which was corrected and extended by Wang et al. (1998) using a spherical formulation. Dodin and Elperin (2002) corrected the Saffman and Turner expression to be valid for the limit cases of inertialess particles without gravity and sedimenting particles without turbulence. However, neither study included the influence of droplet clustering and the coupling between gravity and turbulence, which are of major importance [as pointed out by, e.g., Grabowski and Vaillancourt (1999)]. To this end, Franklin et al. (2007) developed parameterizations based on DNS results of colliding droplets with sizes ranging from 10 to 30 μm, which include these effects. Physically based models with an extended size range were developed by Wang et al. (2000) and Zhou et al. (2001) for turbulent collisions in the case of monodisperse and bidisperse droplet distributions, respectively. However, both models are based on DNS results that used frozen flow fields without gravity. For this reason, the use of these models in a recent cloud physics study by Riemer and Wexler (2005) has brought on discussion in the cloud microphysics community (Wang et al. 2006a; Riemer and Wexler 2006). For example, Franklin et al. (2005) showed that the use of a frozen flow field may lead to an overestimation of the collision kernel by up to 30% compared to an evolving flow field. In conclusion, there is still a need for a physically based model for the geometric collision kernel that includes both the effects of turbulence and gravity.
In this study, we try to contribute to the development of such a model by capturing the essential physical collision mechanisms that occur for droplets moving in a turbulent flow under the force of gravity. We focus on the interplay between turbulence and gravity and the differences between monodisperse and bidisperse droplet distributions. To this end, we employ DNS to study the collision statistics of equal- and unequal-sized droplets, whose motion is driven by both gravity and turbulence. Droplets of radius ranging from 10 to 90 μm are used to study the effectiveness of the different mechanisms in different size ranges. The results are compared to the nonsedimenting case (i.e., where only turbulence is present). The paper aims at a broad understanding of the physical mechanisms that play a role; as such, rather than trying to mimic reality in detail, we choose to simplify certain aspects in an attempt to disentangle the various interplaying phenomena. The effects, for example, of hydrodynamic interactions between droplets will not be accounted for and, to avoid additional complexities, in most situations we will resort to a rather simplistic form of the drag force on droplets; the effect of a more complex drag law is assessed separately.
A limitation of employing DNS is that it cannot resolve the wide range of turbulent scales that is found in convective clouds. As a result, the flow Reynolds number in our DNS is several orders of magnitude smaller than in real clouds. In higher Reynolds number flows, droplet collision statistics may be subject to the intermittent nature of turbulence and to large-scale effects, but in what way and to what extent these effects scale with the Reynolds number is presently an open question, which is not addressed in this study. Another limitation of employing DNS is that because of the limited computational resources available to us, the computational box size is somewhat small. Because of the periodicity of the domain, the somewhat small box size might introduce some artificiality in the results, particularly for the large-scale effects.
Apart from the above limitations, which are associated with the simulation of the continuous phase, the use of point-particle Eulerian–Lagrangian DNS with one-way coupling (i.e., the flow does not feel the droplets) also introduces some limitations. In this approach the droplets are considered as point-particles, and the interaction between a particle and the surrounding fluid is represented through a force located at the position of the center of the particle. Because the droplets are significantly smaller than the Kolmogorov length scale and the concentration is rather dilute, this approach seems appropriate. However, in regions of high droplet clustering, the droplet concentration may become quite high and hydrodynamic interactions between the droplets and the flow modification by the droplets (two-way coupling) might play a role (Portela and Oliemans 2006).
Even though we are interested in the collision statistics and the collision kernel, in this work the actual collision and coalescence processes are not enacted. This is justified when the particle relaxation time is significantly smaller than the time between collisions, which is the case for the dilute concentrations we consider. However, in regions of high droplet clustering the actual collision and coalescence processes might play a role.
In spite of the limitations mentioned above, we believe that this study provides a useful contribution toward a better understanding of the essential physical collision mechanisms that occur for droplets moving in a turbulent flow under the influence of gravity, and as such it may help the development of a physically based model for the geometric collision kernel.
The paper is organized as follows: In the following section, we will treat the theoretical background of colliding droplets in turbulence with and without gravity. Our numerical method is described in section 3. The results of this study are presented in section 4 and summarized and discussed in section 5.
2. Theoretical background
a. Droplet–turbulence interactions: Nonsedimenting droplets
It is well known that the dynamics of finite-sized particles1 such as cloud droplets is very different from the behavior of fluid tracers, when driven by a turbulent background flow field, because these particles possess inertia. The parameter that measures particle inertia is known as the particle relaxation time τp, which is the characteristic time it takes for a particle of radius r to react to changes in the flow. In the case of linear Stokes drag, which is applicable for small cloud droplets (see section 3b), it is given by τp = 2ρpr 2/9ρf ν, where the symbol ρp denotes the density of the particle and ρf and ν are respectively the density and the kinematic viscosity of the carrier fluid. With the use of τp, one can describe the interaction between a particle and the turbulent flow with the Stokes number:
where τf is the characteristic flow time at a scale of interest. In the limit of St → 0, particles act as tracers, following exactly the streamlines of the flow, whereas in the other Stokes limit (St → ∞) droplets are not influenced by changes in the flow. In the nonsedimenting situation (no gravity), maximum droplet–turbulence interaction takes place for intermediate values of the Stokes number, which may cause a distribution of droplets to be altered in its spatial and temporal behavior.
1) Preferential concentration
In a turbulent flow, particles have the tendency to cluster in regions of high strain and low vorticity, a phenomenon that is called preferential concentration (Squires and Eaton 1991; Reade and Collins 2000). This results from the fact that inertial particles are centrifuged out of vortex structures, thus residing in regions of low vorticity. In the absence of gravity, the Stokes number plays a key role in the amount of preferential concentration. From numerical simulations (e.g., Sundaram and Collins 1997; Wang et al. 2000), it was found that the appropriate characteristic flow time is the Kolmogorov time scale τη because it is at this scale that the vorticity is most pronounced. In the limits where the Stokes number is small or large, preferential concentration does not take place. When the particle relaxation time is on the order of τη, however, particles become resonant with the small-scale vortices and cluster in regions of high strain.
Preferential concentration is often measured by a statistical quantity called the radial distribution function g(R) (e.g., Sundaram and Collins 1997), which is defined as the probability of finding a droplet pair at contact, normalized by the corresponding probability for a uniform distribution. Hence, for a uniform distribution, the radial distribution function (RDF) is equal to 1, whereas the RDF reaches values up to 10–100 (e.g., Reade and Collins 2000) when turbulence distorts the uniformity of the concentration field.
For droplets of unequal sizes, the radial distribution function can be related to the monodisperse RDF through the correlation between the two concentration fields (Zhou et al. 2001):
In Eq. (2), the correlation coefficient ρ12 ranges between 0 and 1 and measures to what extent droplets of different sizes cluster in the same regions of the flow. Zhou et al. (2001) found that for two size groups, the correlation decreases with increasing difference in inertia because droplets of different inertia respond differently to the underlying flow field. Hence, the RDF of a bidisperse system is bounded by the values of the monodisperse distributions:
2) Relative velocities
A second consequence of droplet–turbulence interactions on a distribution of droplets is variability in its temporal behavior, expressed through increases in the relative velocities of colliding droplets. Wang et al. (1998) showed that in this context a spherical formulation, in which one measures the radial component wr of the relative velocities, is the correct form to use. This component is given by
where w is the relative velocity vector and r is the separation vector between the two droplets.
Sundaram and Collins (1997) and Wang et al. (2000) found from DNS that the increase in relative velocities from turbulence is highest when τp is on the order of the large eddy turnover time Te. For a bidisperse distribution, the relative velocities are always larger than the values of its monodisperse counterparts, as shown by Zhou et al. (2001):
where 〈·〉 denotes an ensemble average. To explain these observations, one needs to consider the variance of the relative velocity as a measure of the average radial relative velocity:
where and are the mean squared particle fluctuation velocities and 〈υ′1 υ′2〉 represents the covariance of the motion of particles 1 and 2. Let us first treat the relative velocities of a monodisperse distribution. At low Stokes numbers, velocities of equal-sized particles are highly correlated with the fluid and thus with each other, leading to low relative velocities. As the Stokes number increases, the correlation of particle velocities decreases, resulting in higher relative velocities. The increase is bounded, however, because at larger Stokes numbers (τp > Te) particles are less efficient at picking up the energy of the fluid, which ultimately lowers 〈|wr|〉. For unequal-sized droplets, the relative velocities are larger than the monodisperse values because of the low velocity correlation due to the difference in particle inertia. This increase is not unbounded either because or decreases for particles with high inertia.
b. Droplet–turbulence–gravity interactions: Sedimenting droplets in a turbulent airflow
The above considerations apply for finite-sized particles whose motion is only driven by a turbulent flow field. In clouds, however, a gravitational field is present, causing an average downward motion of droplets relative to the flow. Therefore, the cross-trajectory effect reduces the time of interaction between droplets and turbulent structures. The dimensionless parameter that describes this effect is the ratio of an eddy turnover time and the time τυ it takes for a droplet to sediment across that eddy:
where υf is an appropriate velocity scale of the fluid and υt = τpg is the droplet terminal velocity, where g denotes gravitational acceleration. If υt/υf is low, interaction times are scarcely influenced by the average downward motion of droplets. On the other hand, if υt/υf is large, droplets sediment rapidly through a vortical structure and hence exhibit little turbulence interaction. Preferential concentration and relative velocities are thus expected to monotonically decrease with increasing υt/υf , for a given value of St. Obviously, for constant gravity, υt/υf and St are not independent, and as the droplet size changes the influence of both is combined.
Besides decreased interaction times, another phenomenon takes place for sedimenting droplets in a turbulent flow. Because droplets preferentially move on the downflow side of a turbulent eddy, the settling velocity υs of droplets is increased on average. Yang and Lei (1998) showed that this phenomenon, known as preferential sweeping, deals with both small and large scales of turbulence. First, the interaction of droplets with small-scale turbulence, measured by τp/τη, leads to clustering of droplets in regions of low vorticity. In these regions, where the flow velocity is on the order of u′ Yang and Lei (1998), settling droplets are accelerated through drag forces. The increase in the settling velocity is maximum when the terminal velocity υt of droplets scales with the rms of the velocity fluctuations u′.
Dávila and Hunt (2001) suggested that the motion of particles near a vortex can be better understood by combining the inertial and gravitational parameters (i.e., the Stokes number St and nondimensional terminal velocity VT) into a single parameter, the particle Froude number Fp. It is defined as the ratio between the relaxation time of the particle τp and the time for the particle to move around a steady vortex (uf2τf/υt2):
The particle Froude number can be interpreted as a rescaled Stokes number and is a measure of the ratio of the inertial forces that a particle experiences to the buoyancy force.
c. The geometric collision kernel
In a population of droplets, the average number of collisions 12 between two size groups per unit of volume and time is given by
where n1 and n2 are number concentrations of the different size groups; Γ12 is known as the collision kernel,2 which is a measure of the probability of two droplets colliding from geometric considerations. For two droplets of radius r1 and r2, Γ12 can be expressed (e.g., Sundaram and Collins 1997) in terms of the radial relative velocity of the two droplets 〈|wr|〉 and the radial distribution function g(R):
where R = r1 + r2.
In the absence of turbulence, an analytical expression for the collision kernel exists. If gravity is the only inertial force present, droplets sediment with a velocity that is equal to the terminal velocity υt. Wang et al. (1998) derived that in this case, the average radial relative velocity is 〈|wr0|〉 = 1/2 |υt1 − υt2|. Because droplet clustering does not take place, the RDF equals unity and an expression for the gravitational collision kernel is found:
3. Numerical model
a. Flow field
To simulate a flow field that resembles the small-scale features of in-cloud turbulence, we solve the Navier–Stokes equations for an incompressible flow:
where ρf is the fluid density and ν the kinematic viscosity. Using a 3D uniform staggered grid with domain size L and periodic boundary conditions in all directions, Eqs. (12) and (13) are implemented in a finite volume code, which is the modification of a standard channel-flow code described by Portela and Oliemans (2003). Time-step sizes are dynamically restricted by the Courant–Friedrichs–Lewy (CFL) condition, using a Courant number of 0.25.
In Eq. (13), F is the forcing term that maintains the flow to be statistically stationary by continuously adding energy to the system at the lowest wavenumbers. We have chosen to implement a simplified version of the deterministic forcing scheme introduced by Guichard et al. (2004). In the scheme, the linear forcing term is the product of the velocity field in wavenumber space and a time- and wavenumber-dependent forcing function γ(κ, t):
with the wavenumber κ = ‖κ‖ and with ^ denoting the three-dimensional Fourier transform. The forcing function is chosen such that the energy levels of the simulated energy spectrum E(κ, t) are nudged toward a target spectrum Et(κ):
where τf = 0.5ν/εt is a characteristic relaxation delay of E(κ) toward Et(κ) and εt is the target energy dissipation of the system. With the use of the Heaviside function H(κf − κ) a maximum forced wavenumber κf is set.
The flow is initialized with a random flow field that obeys the prescribed target spectrum Et. In our DNS, we use a target spectrum of which the energy level of the first wave band corresponds to the kinetic energy of the large scales and is equal to (0.25εtL)2/3. For κ < κf , the target spectrum scales with κ−1 with the maximum forced wavenumber set at κf = 3(2π/L).
After initiation, a homogeneous and isotropic turbulent flow is generated in a cubic box with sides of L = 0.1 m using a 803 grid. The influence of droplets on the background turbulence is neglected (one-way coupling) because the mass loading and volume fraction of droplets in clouds is small, typically on the order of 10−3 and 10−6, respectively. The viscosity is fixed at a value of 1.5 × 10−5 m2 s−1 (air) and the target dissipation rate εt is set at a typical in-cloud value of 0.05 m2 s−3 (MacPherson and Isaac 1977; Weil et al. 1993; Heus and Jonker 2008). The impact of changing the dissipation rate on the collision kernel was studied by Franklin et al. (2005), among others, and will not be pursued here. After about three large eddy turnover times Te, several statistical quantities are calculated and averaged over a period of 27Te; these are presented in Table 1. The small-scale features of the turbulent flow are completely determined by the fluid viscosity and the dissipation rate through the Kolmogorov scales for length, time, and velocity:
The large scales of the flow, being limited by the domain size of the DNS, are characterized by the fluid velocity rms u′, the integral length scale Le, and the large eddy turnover time Te:
For homogeneous and isotropic turbulence, the Taylor-based Reynolds number Reλ is often used to characterize a turbulent flow. It is some measure of the separation between large and small scales and is based on the Taylor microscale λ:
In Fig. 1a, we show the distribution of energy over different wavenumbers in the DNS in comparison with a model energy spectrum developed by Pope (2000). The graph clearly indicates that energy added at the large scales of the system cascades down to the smallest scales, where it is dissipated. Qualitatively, the DNS spectrum corresponds very well to the model spectrum. In Fig. 1b, the energy spectrum is once again depicted in a log-linear plot, together with the dissipation spectrum. Both spectra are premultiplied by a factor κ for the area under the curve to equal the total energy or dissipation rate. The figure shows that there is a separation of scales between production and dissipation of turbulent kinetic energy in the DNS. In addition, the graphs show that the dissipation range is well resolved by the model, which is an important aspect in the present study.
b. Droplet motion
The equation of motion for a spherical droplet in an evolving turbulent flow is given by Maxey and Riley (1983). Many terms in this equation can be neglected because of the high density of cloud droplets compared to that of air and the small droplet sizes compared to the flow Kolmogorov scales, and the equation reduces to
where V(t) is the droplet velocity at time t, X(t) is the particle position at time t, and g is the gravitational force. The equation above assumes linear Stokes drag, which, strictly speaking, is only valid in the limit of Rep → 0, where Rep is the particle Reynolds number. Because we are dealing with a finite particle Reynolds number, the use of linear Stokes drag does not give an exact representation of the droplet motion, particularly for the larger droplets for which the value of Rep is larger. A standard practice is to introduce an ad hoc nonlinear drag formula in the particle equation of motion. However, these drag formulas are usually empirical fits based on experiments for large spheres under situations that do not necessarily apply to the motion of droplets immersed in a turbulent flow. In this paper, we will mainly confine ourselves to linear Stokes drag to focus on the effects of turbulence and gravity, without taking into account the extra complexity of the nonlinearity in the drag force; that is, any possible effects due to the nonlinearity of the drag force are not mixed in the analysis. The question of whether the results would be significantly altered for an alternative (nonlinear) drag-law formulation is addressed in the appendix. There it is shown that the results are qualitatively similar.
In our point-particle DNS, trajectories of large amounts of cloud droplets of various sizes are computed by numerically integrating Eqs. (20) and (21). Droplet position and velocity are updated using an explicit second-order Runge–Kutta method. To compute the fluid velocity U[X(t), t] at the position of the particle, we use trilinear interpolation in three directions. The time step of particle motion is adapted dynamically, such that a droplet can only travel a maximum distance of half a grid size within one time step. The code is a modification of a code that was used extensively for wall-bounded flows (e.g., Portela et al. 2002), where similar techniques were used. Recently, a benchmark test comparing this code with other codes using higher-order interpolation schemes and different numerical implementation details showed that the statistical results are not very sensitive to the particular numerical implementation details (Marchioli et al. 2008).
Initially, the droplets are placed randomly throughout the domain of the turbulent flow, which has already evolved for about 30Te. The initial velocity of the droplet is the sum of the fluid velocity at the position of the droplet and its terminal velocity. To ensure that the results do not depend on the initial conditions, the droplets are tracked for about 3Te, after which collision statistics are calculated over a period of 57Te.
c. Collision statistics
With the use of an efficient collision searching algorithm that searches only for possible collisions between neighboring droplets (Chen et al. 1999), we can detect droplets that are on a collision course. If the droplets overlap at the time of minimum separation within a certain time step, then the pair is marked to collide. Averaging the total number of collisions in a time interval leads to the average number of collisions per unit time Ṅc, from which the collision kernel can be obtained with the use of Eq. (9) and 12 = Ṅc/V:
where ΓD is the dynamic collision kernel, as it is based on the dynamic detection (i.e., direct counting) of collision events; V is the volume of the computational domain, and N1 and N2 are the number of droplets in a certain size category.
A second way to calculate the collision kernel in a DNS is through the use of Eq. (10):
where ΓK is the kinematic collision kernel. It is based on the computation of the kinematic properties of a collision: the average radial relative velocity 〈|wr|〉 and the radial distribution function g(R).
The radial distribution function g(R) is defined at contact of two droplets and is therefore theoretically only determined by droplets that are exactly a distance R apart. However, because we are dealing with a sparse system, this is almost never the case. To enable a greater number of samples, we identify droplet pairs that are separated by a distance of R ± δ. Wang et al. (2000) showed that a value of δ = R/100 gives accurate results. Then the radial distribution function is given by the probability of finding particle pairs at a distance R ± δ, normalized by the corresponding value for a uniform distribution:
where Npair is the total number of pairs detected and Vs is the volume of a spherical shell: Vs = 4/3π[(R + δ)3 − (R − δ)3]. Averaging over time leads to the average value of the radial distribution function.
The same pair identification procedure is used for the calculation of the relative velocities. Wang et al. (2000) showed that the proper input quantity for the calculation of the collision kernel is the relative velocity of neighboring droplet pairs, rather than the relative velocities of droplets that actually collide. Therefore, we determine radial relative velocities wr of droplets that are separated by a distance between R ± δ using Eq. (4). Further averaging over time gives the average radial relative velocity 〈|wr|〉, which together with the radial distribution function g(R) determines the kinematic collision kernel ΓK.
The algorithms for the calculation of the dynamic collision kernel, the relative velocities, and the radial distribution function have been extensively validated against the gravitational kernel [Eq. (11)], the Saffman and Turner expression for tracer particles (Saffman and Turner 1956), and the expression derived by Abrahamson (1975) for particles with infinite inertia. The results are given in Table 2 and show that the routines for the calculation of the collision statistics are accurate.3 Also, the dynamic and kinematic kernels have been compared for all droplet sizes, and the difference between the two kernels was always smaller than 0.5%, which is comparable to the numerical errors found by Zhou et al. (2001). In this paper, we make use of the dynamic collision kernel when treating collision results. Although we are using expressions such as “droplet collisions” and “collision kernel,” it must be stressed that the actual collision and coalescence processes are not enacted. Thus, droplets are allowed to overlap in space and time. A justification for this exists if we compare the particle relaxation time τp with the time between collisions τc. For one droplet of radius r1 colliding with droplets of radius r2, it is given by τc = (n2Γ12)−1 (Shaw 2003). From previous studies (e.g., Franklin et al. 2005), we know that the turbulent collision kernel is of the order of 10−10 m3 s−1 for a 20–30-μm collision. With the number of droplets we are using (section 4), we find that τp is four orders of magnitude smaller than τc, which leads us to conclude that actual collisions will not change the droplet motion or collision statistics, and we are allowed to employ overlapping droplets.
With the numerical model described in the previous section, a total of 45 numerical experiments was undertaken to study the various collision statistics of droplets with radii ranging from 10 to 90 μm with a step size of 10 μm. Bidisperse distributions were used in our computational domain, where the statistics of the corresponding monodisperse distributions are also calculated. For the simulations in which both droplets were smaller than 50 μm, the amount of droplets in the domain ranged from 333 000 for 10-μm droplets to 67 000 for 50-μm droplets. Simulations with droplets larger than 50 μm were performed with droplet amounts ranging from 114 000 to 23 000, and in all other simulations 32 000 droplets of each size were tracked. Considering the fact that the quantities Γ, 〈|wr|〉, and g(R) are independent of droplet concentrations [cf. Eqs. (22) and (24)], we have chosen these large droplet numbers to decrease the time of statistical averaging, taking into account the time between collisions τc (section 3c), which should be larger than the particle relaxation time τp.
The underlying flow fields were identical in all experiments. We compare our results of sedimenting droplets in a turbulent field with the case of nonsedimenting droplets, as described in section 2a. In Table 3, the various dimensionless parameters that describe droplet–turbulence–gravity interactions are listed for the different size groups.
a. The turbulent geometric collision kernel for sedimenting droplets
1) Monodisperse distributions
To explore the effect of gravity on collision statistics, we first study monodisperse distributions of droplets moving under the force of gravity in isotropic turbulence and compare the DNS results with the case in which only turbulence is present. The comparisons, nondimensionalized by u′, are shown in Fig. 2 for the collision kernel ΓD, the radial relative velocities 〈|wr|〉, the radial distribution function g(R), and the particle kinetic energy (the sum of the mean squared particle fluctuation velocities). The collision kernel has been divided by the factor 2πR2 to remove the triviality that larger droplets have a higher chance of colliding.
Overall, we observe that the presence of gravity decreases the turbulent collision kernel (Fig. 2a) because of the nonnegligible values of υt/υη and υt/u′ that imply decreased interaction times with turbulent structures. This is best shown in Fig. 2b, where we observe that the relative velocities of sedimenting droplets of equal size is much lower than that of nonsedimenting droplets moving in isotropic turbulence. The explanation for this is twofold. On the one hand, sedimenting droplets are less efficient in picking up turbulent energy from the flow, as shown in Fig. 2d, whereas on the other hand the motion of two equal-sized droplets becomes more correlated as a result of the sedimenting behavior. Compared to the nonsedimenting case, both effects lead to a decrease in relative velocities [see Eq. (6)], which is most pronounced for large droplets with a high settling velocity.
Preferential concentration reacts in a different way to the presence of gravity, depending on droplet size, as can be seen in Fig. 2c. For droplets smaller than 40 μm, the presence of gravity decreases the amount of preferential concentration, again because of decreased droplet–turbulence interaction times. This corresponds with the results of Wang and Maxey (1993). However, for droplets larger than 40 μm, the RDF remarkably increases with droplet size and peaks around 70 to 80 μm, whereas in the absence of gravity the RDF drops to values of uniform distributions where g(R) = 1. It is interesting to note that in the region where this peak takes place, the particle Froude number Fp based on large-scale flow characteristics (see Table 3) is on the order of 1. Furthermore, this region corresponds to droplet sizes that experience an extra decrease in relative velocities (Fig. 2b). The latter observation can be explained by the reasoning that the motion of clustered droplets is more correlated, leading to lower relative velocities. It seems that the presence of gravity allows for interaction between large sedimenting droplets and large energetic eddies (τp/Te ∼ 0.5), a phenomenon that does not occur in the absence of gravity. As a result, these droplets cluster in certain regions of the flow, which we will explore further in section 4c.
2) Bidisperse distributions
For nonsedimenting bidisperse droplet distributions, we know from previous studies (e.g., Zhou et al. 2001) that preferential concentration decreases as the difference in droplet size becomes larger, while relative velocities are increased. In Fig. 3, we show the turbulence-induced increases of the collision kernel, the radial relative velocities, and the radial distribution function for various droplet sizes, where we have subtracted the theoretical gravitational kernel [Eq. (11)] to allow for a fair comparison with the nonsedimenting case. In the contour plots, light regions denote large increases and dark regions denote small increases.
We now examine the graph of the relative velocities (Fig. 3b), which clearly shows a region of maximum increase for bidisperse distributions. From Eq. (6), we know that there are two effects that determine the increase in relative velocities Δwr = 〈|wr|〉 − 〈|wr0|〉. First, the velocity decorrelation, determined by the inertial difference of the two droplets, accounts for the fact that Δwr of bidisperse distributions is larger than the monodisperse values. The presence of gravity enhances the decorrelation due to the different sedimenting behavior, resulting in a fast increase around the monodisperse region. The second effect, however, counteracts this increase because the particle kinetic energy decreases with increasing droplet inertia. Again the effect is enhanced by gravity, as shown in Fig. 2d. The combination of both effects, being dominant in different size regions, explains the presence of a distinct region of maximum increase.
It is interesting to note that the effect of gravity on the correlation of the motion of two droplets 〈υ′1 υ′2〉 is different for the monodisperse and bidisperse cases. For equal-sized droplets, the sedimenting behavior increases the velocity correlation, whereas 〈υ′1 υ′2〉 is decreased for droplets of unequal size in the presence of gravity because of the difference in terminal velocities.
Figure 3c displays the amount of preferential concentration for bidisperse droplet distributions. Similar to the case in which only turbulence is present (as described in section 2a), we observe that the RDF drops to lower values for droplets of unequal size. Physically, these droplets reside in different regions of the flow because of different inertial responses to flow accelerations. As a result, the correlation ρ12 between the concentration fields is low [Eq. (2)]. For larger droplet pairs, the isocontour lines diverge because the values of the monodisperse distributions are higher as well (Fig. 2c).
One should bear in mind that the data behind the graphs in Fig. 3 is rather coarse, especially around the region of monodisperse distributions (the diagonals from bottom left to upper right) where smoothing of the raw data may lead to an inaccurate representation of the situation. To investigate this issue, we make a more detailed analysis of the region between monodisperse and bidisperse droplet distributions in the next section.
b. The region between monodisperse and bidisperse distributions
Both from theoretical considerations and from the observations in the previous section, it was found that in going from monodispersity to bidispersity, Δwr = 〈|wr|〉 − 〈|wr0|〉 increases, while the amount of preferential concentration g(R) decreases. As can be seen in Fig. 2, the collision statistics are quite pronounced around r ≈ 40 μm. For r ≈ 40 μm, τp/τη ≈ 1, and for nonsedimenting droplets it is known that the collision statistics are more pronounced for τp/τη ≈ 1. Therefore, to examine this region more thoroughly, we zoom in by studying collisions of droplet pairs such that the radius of one of the droplets ranges from 40 to 42 μm and the other is fixed at 40 μm. The results are shown in Fig. 4 for the collision kernel ΓD, the radial relative velocities 〈|wr|〉, the radial distribution function g(R), and the correlation coefficient between two concentration fields ρ12 [Eq. (2)]. The dashed lines represent the nonturbulent gravitational case [Eq. (11)]; hence, the difference between the two lines corresponds to Fig. 3.
In Fig. 4b, we see that in the small region of r2 = 40–42 μm, 〈|wr|〉 linearly increases with increasing droplet radius, mostly because of the increased difference in terminal velocities of the two droplets. The linearity continues when the graph is extended to r2 = 50 μm, which is not shown here. Considering the increase of 〈|wr|〉 compared to the nonturbulent gravitational case, we observe that Δwr increases as well with droplet radius. The increase in the difference between the two lines is quite large: a relative difference in r of 5% leads to an increase in Δwr that is 3.69 times as large as the monodisperse case.
The behavior for preferential concentration is even more explicit, as shown in Fig. 4c. For droplets of equal size, the amount of preferential concentration is around 20. However, a small difference in droplet size leads to a large drop in g12(R). For example, when r2 = 1.05r1 the RDF is reduced to a value of 2.56, which is 7.88 times smaller than the monodisperse value. Comparable results were obtained by Wang et al. (2006a), who showed a sudden discontinuous peak in the RDF at r1 = r2. Here we have taken the size resolution such that this peak could be entirely resolved, providing detailed information on the extent of the monodispersity effect (less than 1 μm for 40-μm droplets). An explanation for this striking effect is found in Fig. 4d, where we have plotted the correlation coefficient between the concentration fields of the two droplets. This coefficient measures to what extent droplets reside in the same regions of the flow. For r2 = r1 it is equal to 1. The graph shows that a small difference in droplet size leads to a large decrease in ρ12 and hence to regions of clustering of the two distributions that do not overlap. This leads us to conclude that in the presence of gravity, the regions where droplets of a certain size concentrate can be very sensitively determined by droplet size, and even a small difference in size makes a distribution of droplets cluster in a very different regions.
Because the behavior of sedimenting and nonsedimenting droplets is clearly very different, we study more closely the effect of the gravitational acceleration on the results by increasing g in steps from 0 to 10 m s−2. We do so, not because it is realistic for cloudy situations, but because it provides a continuous link between the “extreme” cases of nonsedimenting and normal sedimenting droplets. In Fig. 5 we show how the amount of gravitational acceleration influences the RDF for the droplet pair r1 = 40 μm and r2 = 41 μm. In correspondence with Fig. 2c, the RDF of equal-sized droplets around r = 40 μm (dashed lines in Fig. 5a) is not significantly influenced by changes in the gravitational acceleration. However, the RDF of the corresponding bidisperse distribution g40−41(R) shows a dependence on gravity as it decreases with increasing gravitational acceleration because of the decreasing correlation coefficient of the two concentration fields (Fig. 5b). Physically, one can explain the gravitational effect on the correlation with the notion of droplet–turbulence interaction times, which are determined by the droplet terminal velocity υt. For two droplets of unequal size, increasing the gravitational acceleration implies enhancing the difference in interaction times between the two droplets. As a result, the droplets interact in a dissimilar way with the turbulent structure, leading to an enhancement in the decorrelation of the two concentration fields.
The combination of the effect of relative velocities and preferential concentration is shown in Fig. 4a for the collision kernel. The result is a monotonically increasing collision kernel that increases at the same rate as the gravitational kernel. Thus, the probability of collision between two droplets does not experience any dramatic changes when going from monodisperse to bidisperse distributions. Coming back to the possibly incorrect representation of this region in Fig. 3, we conclude that Figs. 3a and 3b are presented correctly, whereas in Fig. 3c the peak in the RDF around the monodisperse region should be narrower. The width of the peak depends on the amount of gravitational acceleration.
c. Clustering of large droplets in the presence of gravity
In Fig. 2c, we noticed that in contrast to the turbulence-only case, the presence of gravity allows for clustering of droplets with radii r greater than 50 μm. Interestingly, the amount of preferential concentration peaks for droplets with a particle Froude number Fp of the order of 1, based on large-scale turbulence characteristics. This remarkable large-scale clustering in the presence of gravity will be explored further by examining collision statistics of 70- and 80-μm droplets (both mono- and bidisperse) as a function of gravitational acceleration g, the results of which are presented in Fig. 6.
First, the relative droplet velocities 〈|wr|〉 are shown in Fig. 6b. Both for the monodisperse and bidisperse distributions, 〈|wr|〉 shifts with increasing g toward the gravitational value 〈|wr0|〉, which is zero for equal-sized droplets. This results from the fact that sedimenting droplets are less efficient in picking up energy from the turbulent flow as their settling velocity increases, lowering the particle kinetic energy (Fig. 2d) and hence the relative velocities. Moreover, one may expect a higher correlation of motion of two droplets due to high levels of preferential concentration in the presence of gravity (cf. Figs. 2b,c for r > 50 μm) to further reduce 〈|wr|〉.
Figure 6c shows how the amount of preferential concentration varies with gravity for droplets of sizes 70 and 80 μm. The dashed lines (monodisperse) clearly show that it is the gravitational force that is responsible for the clustering of large droplets. The bidisperse distribution (solid line) shows a point of minimal clustering as the RDF slightly decreases from a value of 1.90 to 1.52 for g = 2.50 m s−2. This is due to the decrease in the correlation coefficient ρ12 as shown in Fig. 6d, comparable to what we observed earlier for the droplet pair r1 = 40 μm and r2 = 41 μm in Fig. 5b. Further increasing the gravitational acceleration leads to a small increase in the RDF to a value of g(R) = 3.31 for g = 9.81 m s−2. Because the correlation coefficient remains more or less constant, this increase can be attributed to the large increase of the amount of preferential concentration of the monodisperse distributions as shown in Fig. 6c. To gain a better understanding of the large-scale droplet clustering of monodisperse distributions in the presence of gravity, we show 2D and 3D snapshots of a droplet concentration field of size r = 80 μm in Fig. 7. The graphs clearly show that this clustering involves the larger scales of turbulence, considering the large depleted regions of droplets. Furthermore, the droplet interaction with large-scale vortices forms streaks of droplets in stream and convergence zones of the flow (Squires and Eaton 1991) in which the droplets move downward, comparable to the preferential sweeping mechanism.
In conclusion, it seems that the presence of gravity allows for interaction of large droplets (τp/Te ∼ 0.5) with the larger scales of turbulence, resulting in the formation of large “curtainlike” manifolds of sedimenting droplets waving in the turbulent field. Some care must be taken when interpreting these results, however, considering that the periodicity of the domain may introduce some artificiality. In Table 3, we have calculated the ratio of the time it takes a sedimenting droplet to cross the domain TL and the large eddy turnover time Te. A ratio of TL/Te < 1 implies that a droplet may encounter the same large-scale vortex multiple times, which is unrealistic. On the other hand, in reality vortices of larger length and time scales are present, which may actually facilitate the large-scale organization of droplets. This can be tested in future DNS studies by increasing the domain size.
5. Summary and discussion
Numerous studies have previously considered collisions of inertial particles in turbulent flows. However, as discussed in a recent review paper by Khain et al. (2007), in most studies the combined effect of turbulence and gravity is not considered, which can lead to unrealistic cloud models. In this paper, by means of direct numerical simulations, we have studied the collision statistics of equal- and unequal-sized droplets ranging in radius from 10 to 90 μm, where the droplet motion is determined by both gravity and drag in isotropic turbulence. Our focus was on the mechanisms of preferential concentration and enhanced relative velocities and how they are affected by the presence of gravity, compared to the turbulence-only (i.e., nonsedimenting) case. Furthermore, we have examined in more detail how the influence of gravity determines the region between monodisperse and bidisperse distributions and how gravity allows for interaction between large droplets and large scales of turbulence.
It was found that for monodisperse distributions, gravity decreases the collision probability compared to the nonsedimenting case because of decreased interaction times between droplets and turbulent structures. For bidisperse distributions, however, the collision kernel is increased by gravity as a result of the fact that droplets of unequal size possess different terminal falling velocities. In the region between monodispersity and bidispersity, the amount of preferential concentration, measured by the radial distribution function, drops to low values very quickly as the difference in droplet size becomes larger. This finding can be compared to the sudden peak observed by Wang et al. (2006a) in the radial distribution function for equal droplet radii. By varying droplet sizes in steps of 0.1 μm, we were able to fully resolve this peak and get information on the extent of this monodispersity effect (∼1 μm for 40-μm droplets). By gradually changing the amount of gravitational acceleration, it was shown that gravity plays a crucial role in this matter, causing droplets with slightly different sizes to concentrate preferentially in different regions. In addition we observed that the presence of gravity allows for preferential concentration of large droplets (τp/Te ∼ 0.5), which clearly does not take place in the nonsedimenting case. It appears that the preferential sweeping mechanism described by Wang and Maxey (1993) is also rather effective for relatively large droplets, which on account of the presence of gravity engage in an interaction with large-scale structures of the flow.
With regard to the geometric collision probability of droplets in turbulent clouds, two major conclusions can be drawn. 1) The combined effect of turbulence and gravity is not merely an addition of separate phenomena; the interplay between turbulence and gravity strongly affects the behavior of the collision statistics. Examples include the correlation between the concentration fields of a bidisperse distribution being determined by the amount of gravitational acceleration and the clustering of large droplets interacting with large-scale vortices in the presence of gravity. 2) A distinction must be made between monodisperse and bidisperse droplet distributions for the collision statistics of droplets moving under the forces of turbulence and gravity. Different behavior of, for instance, the velocity correlation of droplet pairs of equal and unequal size was observed, and we found a vast difference in the amount of preferential concentration for the two different cases.
Furthermore, our results suggest that not only small-scale turbulence influences collision statistics, but also that the large scales have an effect. Therefore, a next study should investigate the Reynolds number dependence of the collision statistics in the presence of both turbulence and gravity. In this matter, there are still two issues that are open for debate. First, the nonlinear droplet collision–coalescence process may be sensitive to local changes in the collision kernel because of the more intermittent nature of turbulence at high Reynolds numbers (Wijngaard 1992; Sreenivasan and Antonia 1997). Second, Wang and Maxey (1993) and Yang and Lei (1998) observed that settling velocities of sedimenting droplets display a Reynolds number dependence. The larger range of scales involved may affect the relative velocities of colliding droplets and thus the collision kernel. For the interaction between large droplets and large-scale vortices, some caution is needed as well, considering the artificiality that may be introduced by the large-scale forcing and the periodicity of the domain. Moreover, the assumption that the particle Reynolds number and the ratio between particle relaxation time and average time between collisions τp/τc are smaller than unity is questionable for droplets of 80- and 90-μm radii. Hence, the motion of these droplets does not correspond entirely to reality and the phenomenon of large-scale clustering should be interpreted with all the abovementioned limitations in mind.
As pointed out recently by Xue et al. (2008), further research is needed before a complete kernel model can be developed. The results of the present study may contribute to the development of a physically based model of the turbulent collision kernel that accounts for gravity. Such a model, which should include all droplet, flow, and gravitational parameters (ε, Reλ, r, g) as we have shown, can be used in the numerical solution of the stochastic collection equation that describes the temporal behavior of a droplet size distribution in a cloud. This solution, based on the static calculation of collision statistics, could be compared to the dynamic case in a DNS where colliding droplets actually coalesce and form a larger droplet.
This work was sponsored by the National Computing Facilities Foundation (NCF) for the use of supercomputer facilities, with financial support of NWO. The authors gratefully thank Dr. J. C. R. Hunt for discussions concerning the particle Froude number, Dr. W. Grabowksi for helpful comments on an earlier version, and two anonymous referees for their critical remarks and suggestions to improve the manuscript. The authors also thank Dylan Dussel and the Virtual Reality lab of TU Delft for the insightful visualizations of the large-scale clustering.
Nonlinear Drag Simulations
As discussed in section 3b, we opted to use linear Stokes drag for the equation of motion of the droplets in order to focus on the effects of turbulence and gravity without taking into account the extra complexity of the nonlinearity in the drag force. Here we will assess the effects of a nonlinear drag-law formulation on the results. Considering only drag and gravity, the equation of motion of the particles can be written as
where τp is the linear Stokes drag particle relaxation time, given by τp = 2ρpr 2/9ρfν, and f is the drag factor, defined as the ratio between the drag coefficient and the drag coefficient for Stokes drag; that is, f is given by
where CD is the drag coefficient and Rep is the particle Reynolds number, defined as
We repeated the simulations using Eq. (A1) instead of Eq. (20), with f given by Eq. (A4). The monodisperse results are shown in Fig. A1 and should be compared with Fig. 2 (repeated in Fig. A1 for convenience). The bidisperse distributions, to be compared with Fig. 3, are shown in Fig. A2. One can observe roughly the same behavior for the entire range of droplet radius considered. That is, the differences for the two drag-law formulations are relatively small when compared to the differences between the sedimenting and nonsedimenting situation.
The roughly similar behavior for the linear and nonlinear drag can be explained if one realizes that the major effect of the nonlinear drag is to decrease the droplet terminal velocity and that for larger droplets the main contributor to Rep is the droplet terminal velocity (for r = 10, 50, and 90 μm, the value of Rep associated with the droplet terminal velocity is, respectively, 0.016, 1.6, and 7.3 for nonlinear drag, and 0.016, 2.0, and 12, for linear Stokes drag), whereas for smaller droplets the nonlinear component of f is negligible. If the main contributor to Rep is the terminal velocity, once the transient disappears the value of Rep is roughly constant, and the droplet equation of motion can be written as
where V′ and U′ denote the deviations of V and U. The effective equation of motion for nonlinear drag is the same as the equation of motion for linear drag but corresponds to a smaller droplet radius (i.e., a larger value of f/τp). Hence, at least within a moderate range of droplet sizes, one would expect that the droplet velocity fluctuations for nonlinear drag would exhibit roughly the same qualitative behavior as for linear drag, with the values corresponding to a smaller radius; this can be seen, for example, in Fig. A1.
Corresponding author address: Dr. Harm Jonker, Dept. of Multi-Scale Physics, Delft University of Technology, P.O. Box 5046, 2600 GA, Delft, Netherlands. Email: email@example.com
In this section, we will often refer to cloud droplets as particles.
From here, when the term “collision kernel” is used, we mean to describe the geometric collision kernel Γ12.
The observed deviation of 9% in the zero-inertia case is not considered serious because the theoretical value is based on the simplified assumption of Gaussian velocity statistics.