## Abstract

Binary collisions of large raindrops moving with terminal fall velocity are numerically investigated using FS3D, a direct numerical simulation (DNS) code based on the “volume of fluid” method. The result of this process can be a permanent coalescence or a temporal coalescence followed by a breakup of the coalesced system into smaller-sized remnants of the original raindrops and a number of fragment droplets of different sizes. In total, 32 drop pairs are studied with sizes chosen to cover nearly completely the entire size parameter range relevant to breakup. This is an important extension of investigations performed in 1982 by Low and List, who studied 10 drop pairs only. Moreover, eccentricity has been introduced as an additional parameter controlling the collision outcome. Eccentricity is defined as the horizontal distance of the initial drops’ centers with values equal to approximately 0 for centric and 1 for grazing collisions.

The main results include numerically calculated data of coalescence efficiencies and fragment size distributions with emphasis on eccentricity effects. It is shown that eccentricity largely determines the appearance of specific breakup modes and consequently the respective fragment size distributions. Comparisons are made with the main findings of Low and List. Coalescence efficiency values larger than those derived by Low and List show up for very small Weber numbers. Additionally, the existence of their definite limit value of collision kinetic energy necessary for coalescence could not be confirmed. The fragment size distributions are in some cases similar to those measured by Low and List but there are also major differences for other cases. The presented results are used for parameterizations of coalescence efficiencies and fragment size distributions as well as for calculations of stationary drop spectra shown in Part II of this study.

## 1. Introduction

Number density size distributions of raindrops observed at ground frequently show a maximum at diameters of about 1–2 mm with a nearly exponential decrease for larger raindrops and a drop of the spectra for smaller raindrops (Marshall and Palmer 1948; Du Toit 1967; Zawadski and de Agostinho Antonio 1988). Such spectra are often assumed to be mainly the result of two processes counteracting each other: large raindrops are produced by collision/coalescence of smaller drops, whereas large raindrops are diminished in size by breakup whereby smaller drops are created. It should be noted, however, that evaporation, taking place in undersaturated subcloud layers and thus diminishing the number and mass of small raindrops, can also shape raindrop spectra (Mason and Ramanadham 1954; Li and Srivastava 2001). This effect is neglected here.

Whereas collision/coalescence is more or less well understood and quantified—except for turbulence effects influencing this mechanism—our knowledge of breakup is limited. With regard to breakup, two single processes are distinguished: hydrodynamic (or spontaneous) breakup of single very large raindrops and collision-induced breakup following from binary collisions of large raindrops. The former possibility is mostly ruled out since raindrops have to be very large to show this effect. Consequently, collision-induced breakup is a candidate to inhibit the formation of very large raindrops and is considered to be the principal mechanism limiting the size of raindrops (Pruppacher and Klett 1997). Despite the importance of collision-induced breakup, only a limited number of investigations exist to provide data and appropriate parameters to be used in meteorological applications.

Before reviewing past and recent efforts on better understanding collision-induced breakup, we sketch the process in more detail. A collision itself can lead, on the one hand, to a permanent coalescence of two initial raindrops such that one final larger raindrop is formed. On the other hand, a temporal coalescence may take place followed by breakup of the coalesced system such that smaller-sized remnants of the original raindrops and a number of fragment droplets of different size appear. Thus, the latter mechanism has—besides reducing the number and size of large raindrops—a side effect, namely that the fragment droplets become part of a pre-existing drop spectrum of a cloud. This feedback mechanism may considerably influence cloud and precipitation evolution (Seifert et al. 2005).

There are a large number of publications on the topic of drop–drop collisions since this process is—in addition to the meteorological context—important in several disciplines including technical applications in sprays, polymer blending, and other multiphase flows. In such studies the conditions for permanent coalescence are equally addressed as those for different types of breakup modes as reflexive and stretching separation. In this respect, one of the most cited papers is that of Ashgriz and Poo (1990), which deduces criteria as to whether coalescence or separation can be expected as the outcome of a collision. Based on their investigations, several other authors addressed this process (e.g., Qian and Law 1997; Brenn et al. 2001), concentrating mostly on the collisions of equal-sized drops of different fluids by means of laboratory experiments supplemented by theoretical considerations. Orme (1997) gave an outline of experiments and in Menchaca-Rocha et al. (1997) some theoretical considerations and approaches are compiled.

Related to cloud microphysics, Adam et al. (1968) sketched the limit between a stable and an unstable outcome of colliding equal-sized drops of 60 and 300 *μ*m based on evaluations of laboratory experiments in which the drop trajectories were horizontal. Some years later, Brazier-Smith et al. (1972) and Bradley and Stow (1978) presented thorough experimental and theoretical investigations on the interaction of water drops colliding in free fall. They primarily studied conditions critical to drop separation following a collision event. Later on, very comprehensive studies were done by McTaggart-Cowan and List (1975), Low and List (1982a, hereafter LL82a), Ochs et al. (1986), and Czys and Ochs (1988) performing laboratory measurements. They determined mainly coalescence efficiencies and fragment size distributions. However, these experiments suffered from a strong limitation of the number of dual drop size combinations. Specifically, LL82a only considered 10 pairs of large raindrops, and Ochs et al. (1986) and Czys and Ochs (1988) studied four pairs of small raindrops. Very recently, Barros et al. (2008) presented laboratory experiments on collision-induced breakup with a larger range of drop pairs than that considered by LL82a. They showed that their fragment size distributions were slightly different from those determined by LL82a; however, they confirmed that a certain limit for coalescence appears as given by LL82a.

From their sparse database, Low and List (1982b, hereafter LL82b) derived parameterizations of both coalescence efficiencies of large drops and fragment size distributions where the latter are the only ones applied in current cloud models simulating precipitation processes in greatest detail by a spectral (bin) formulation. Note that recently McFarquhar (2004) reinvestigated the data of LL82a and arrived at an extension and improvement of the original parameterization given by LL82b.

It should be clear from this discussion that the parameterizations of LL82b^{,}^{1} [and of Beard and Ochs (1995), which are not further discussed in this study] rely on a dataset that is very limited with respect to the size range of the two colliding drops. This fact was a motivation to create an extended database by considering an almost complete size range of colliding drops for the derivation of an improved parameterization describing collision-induced breakup.

The first steps to overcome the mentioned deficiencies have been done by Beheng et al. (2006), who computationally simulated binary collisions for the same drop pairs as LL82a and for some additional ones. The present study is a further extension of that work by again increasing the range of drop pair combinations in order to broaden the database. In contrast to LL82a, we also account for the influences of eccentricity on breakup.

This paper is organized as follows. At first the governing equations, the numerical method and setup, and the simulation procedure are presented. In the results section the collision outcomes as well as the budget of kinetic energies are analyzed. This is followed by a discussion of coalescence efficiencies and the number and size distribution of the fragment drops, with special emphasis on the influence of eccentricity.

In the companion paper (Straub et al. 2010, hereafter Part II), we concentrate on deriving parameterizations based on our results. To this end, we provide parameterizations for the coalescence efficiencies and fragment drop size spectra. That paper concludes with a demonstration of the impact of the new parameterizations on the shape of a stationary raindrop spectrum as compared to application of the LL82b parameterization.

## 2. Formulation and numerical method

The numerical code applied is the DNS-CFD code Free Surface 3D (FS3D), which was developed many years ago at the Institut für Thermodynamik der Luft- und Raumfahrt (ITLR). FS3D has been used for numerous investigations of drop dynamics and other basic multiphase flow simulations and is well validated. Of the abundant publications in which this code has been applied, we only refer to Rieber and Frohn (1999), who successfully simulated the process of drops splashing onto a liquid film, and to Gotaas et al. (2007), who studied the effect of viscosity on binary droplet collisions. In both cases the results were consistent with laboratory observations. In this section we give a short outline of the governing equations and numerical solvers used in FS3D. More detailed information can be found in Rieber (2004).

For the fluid system to be simulated, the governing equations are the balance equations for momentum and mass. The considered momentum balance reads

where ⊗ is a dyadic product, **u** is the velocity and *t* denotes time; *ρ* and *μ* are the local density and viscosity determined by the “volume of fluid” (VOF) variable *f* (see below), and *p* is the pressure; **k** represents an external body force, which in the present case is gravity. Surface tension is included by the volume force **f**_{γ}.

For an incompressible medium such as water the full continuity equation degenerates to the condition of a nondivergent flow:

Possible condensation–evaporation of the drops is omitted here. This is an admissible simplification as long as the saturation deficit between the liquid and the surrounding gaseous environment is small. So we assume very humid to saturated conditions valid inside clouds. Note that with FS3D condensation–evaporation can basically be calculated (Schlottke and Weigand 2008).

The representation of different phases in FS3D is based on the VOF method by Hirt and Nichols (1981). To distinguish between the liquid and the gaseous phase, the additional field variable *f*, expressing the volume fraction of the liquid in a region containing a mixture of both fluid and gas, is introduced. Here *f* is defined according to

and for its spatial and temporal change the following conservation equation is applied:

To maintain a sharp interface, the convective transport of *f* is performed based on the reconstructed interface using the piecewise linear interface calculation (PLIC) method of Rider and Kothe (1998), minimizing numerical diffusion. Figure 1 sketches the basic idea of the PLIC method.

The fluid properties are assumed to be constant in each phase but vary in regions containing the interface. They are obtained by applying the *one-field* formulation; that is, viscosity and density are calculated using

where the subscripts *l* and *g* denote the liquid and gaseous phase, respectively. Fluid properties could be treated as temperature dependent with FS3D, but for the present calculations this effect is negligible.

The equations are spatially discretized using a second-order accurate finite volume method on an equidistant staggered grid. The second-order temporal integration is performed using the Crank–Nicholson method. A fast and robust multigrid solver is used for the projection of the emerging preliminary velocity field onto one fulfilling Eq. (2).

## 3. Numerical setup and simulation procedure

This section illustrates our approach to simulate raindrop collisions by means of computational grids, initial placement, shape, and velocity of the drops and simulation duration.

The setup mimics that a large raindrop positioned initially some distance apart from a smaller one moves because of its greater terminal velocity toward the smaller one. After contact is encountered, the surfaces and masses of both drops form an intermediate system that may oscillate, giving rise to a complex flow field in the interior. As a consequence of the delicate interplay between kinetic and surface energies, as well as the internal and surrounding flow fields, protuberances may appear which thin out and eventually break up. The latter leads to the formation of fragment droplets departing from the parent drops. The computations are stopped after a stable configuration is established.

All simulations were performed using a 3D computational domain with free-slip conditions at the lateral boundaries and periodic boundary conditions on the bottom and top of the domain. This setup permits a numerically justifiable effort for the large number of conducted simulations because we can observe the collision process for a rather long time period in a limited computational domain. The spatial resolution was chosen to be 100 *μ*m for all cases except investigations of grid dependency, where the spatial resolution was 50 *μ*m. According to the drops’ sizes, we used different grids in order to have a sufficient distance between the drops and the lateral boundaries. This resulted in approximately 16 × 10^{6} cells for the largest grid volume (2.56^{3} cm^{3}) and standard calculations as well as nearly 135 × 10^{6} cells for grid dependency investigations (see section 4g).

The drops are assumed to move in a vertical direction with terminal fall velocity *υ*, which is calculated following Beard (1976). This is slightly different from LL82a and LL82b, who used the relation of Best (1950). However, this leads to only small differing values of certain parameters, such as collision energies and consequently Weber numbers (see below).

Figure 2 sketches the initial setup on a plane through the drops’ centers. The initial vertical distance between the drops was set to Δ*y* = *υ _{r}*Δ

*t*′ (

*υ*is relative velocity; see below) with Δ

_{r}*t*′ = 0.5 ms, which allowed the drops to interact to some extent with their surrounding flow field in the early beginning of the simulation. The nonspherical shape of large drops was considered by initializing them as ellipsoids with the ratio of the equatorial radii to the polar radius set according to Beard and Chuang (1987). The dashed circle in Fig. 2 indicates the same volume as that of the ellipsoid. As already mentioned, eccentricity

*ε*is taken into account here. This parameter—defined as the ratio between the distance of the drops’ centers

*δ*to the arithmetic mean of the diameters of volume equivalent spheres

*d*(with subscripts

*L*and

*S*for the large and the small drop, respectively),

has been varied from 0.05 to 0.95 in six discrete steps (*ε* = 0.05, 0.2, 0.4, 0.6, 0.8, and 0.95). Bulk properties for water and air were chosen at 20°C with a surface tension of *σ* = 73 × 10^{−3} N m^{−1}. This corresponds to the conditions of LL82a.

It should be noted that because of aerodynamic effects of the oncoming flow on the deformed drops, collision can occur even for *ε* > 1. Moreover, occasionally the large deformed drops oscillate initially before contact, so *ε* may not have exactly its prescribed value. But this inaccuracy is assumed to be of minor importance according to LL82a.

The simulations were performed until a stable outcome of the collision process was reached (i.e., when the number of fragments stopped changing). This condition was additionally assessed by visual inspection of the results. The longest collision process took about 25 ms. The average time step Δ*t* was ≈3 × 10^{−6} s.

The number and sizes of the fragments are determined by a region-growing approach. Mass conservation was satisfied for all cases. To compare the present results to literature data, the obtained results for different eccentricities have been weighted according to their impact probability (see below).

## 4. Results

Overall, 32 drop pairs are investigated. A compilation of the corresponding initial drop diameters is given in Table 1, where the first 10 drop pairs are those of LL82a, drop pairs 11–18 are from Beheng et al. (2006), and the remaining drop pairs (19–32) are taken into account to cover the matrix of pairs as completely as possible. The choice was based on the results of Beheng et al. (2006), which already exhibited certain structures of the fragment size distribution, notably the occurrence of multipeaked spectra depending on drop sizes. A matrix of the considered drop pairs is shown in Fig. 3 with the LL82a drop pairs highlighted by boxed symbols. Note that the diameter of the smallest drop considered is 0.035 cm (i.e., it is evidently a raindrop). Fragment drops are usually smaller in diameter and may be termed cloud drops, although no definite limit between the sizes of cloud drops and raindrops is applied here.

After introducing some necessary definitions of the relevant energies and their connection to the drop pairs, we analyze and discuss the results of the collision simulations. To this end, images of some collision events are shown first. Then the time evolution of the energy budget and surface area during this process is presented as well as the coalescence efficiency, the number of fragments, and their size distribution after collision with the influence of eccentricity on the aforementioned issues. We finally assess the presented results in terms of grid dependency and sensitivity to the boundary conditions.

### a. Parameters

The results that follow are characterized by specific parameters adopted from LL82a.

If bouncing is neglected, the collisional interaction of two drops forming a single larger drop is, physically speaking, a completely inelastic event. This process is ruled by the conservation of momentum and kinetic energy before and after the collision. Based on these conservation principles, the collision kinetic energy (CKE, as denoted by LL82a) is given by

where *ρ _{l}* is the bulk density of water and

*υ*=

_{r}*υ*−

_{L}*υ*the relative velocity, with

_{S}*υ*and

_{L}*υ*being the velocities of the large and small drop, respectively. Similarly, surface energies are related to give the difference by

_{S}where *S _{T}* is the sum of the surface energies of the initial drops and

*S*the surface energy of the coalesced (spherical) system. Rotational kinetic energy is excluded here because it is at least one order of magnitude smaller than CKE and Δ

_{c}*S*(see Part II for an explicit formulation of the rotational kinetic energy).

Summing up both energies yields the total energy difference (excess energy):

A dimensionless parameter characterizing effects by kinetic relative to surface energies is the Weber number We. In the present context it is defined by

Values of CKE and We for all 32 drop pairs are compiled in Table 1.

Note that in engineering and fluid dynamics applications the Weber number is defined as

Both Weber numbers are related to each other by

with *q* = *d _{S}*/

*d*.

_{L}As already mentioned earlier, terminal velocities of the drops have been calculated following Beard (1976), which slightly differs from LL82a, who applied the formula given by Best (1950). In scrutinizing the differences coming from usage of the different velocity formulations, it turned out that only minor discrepancies in regard to the numerical values of CKE and We appear; in fact, they closely correspond to those of LL82a.

### b. Collision outcomes

Results of the numerical simulations in the form of snapshots are shown exemplarily in Figs. 4 and 5 for two different drop pairs considering five eccentricity values (note that outcomes for *ε* = 0.95 are not shown because of space limitations). One pair chosen for this comparison has a high collision energy (pair 10, CKE = 12.53 *μ*J) and the other has a low one (pair 16, CKE = 3.93 *μ*J). The breakup modes (disk, sheet, and filament) are named according to LL82a. Note that this classification can be rather subjective but is carried out in order to compare to the results of LL82a on the one hand and to improve the understanding of the relevant physics on the other hand.

For pair 10, shown in Fig. 4, the dependency of the breakup mode on eccentricity is evident. As expected for the considered high collision energy, there is disk breakup for *ε* = 0.05 and *ε* = 0.2 (i.e., nearly central collisions, resulting in numerous small fragment droplets of different sizes after the collision process). For intermediate eccentricities, *ε* = 0.4 and *ε* = 0.6, the breakup mode is sheet breakup and the size distribution of the fragment droplets is more homogeneous. For the strong off-center to grazing cases, *ε* = 0.8 and *ε* = 0.95 (not shown), finally there is filament breakup: the initial drops persist after the collision process with almost no change in drop size but few additional very small fragment droplets are produced. The results concerning this example are common to all cases with relatively high CKE values. (An overview of the dependency on eccentricity is also shown in Fig. 14 for the discussed droplet pair and two additional ones with lower CKE.)

The collision outcomes for pair 16 (an example of a pair with low CKE) are shown in Fig. 5. For *ε* = 0.05 and *ε* = 0.2, the two drops form a single larger permanent drop after collision (coalescence case). For larger eccentricities, there is filament breakup with no additional droplets for *ε* = 0.4 and *ε* = 0.8 but with some small fragment droplets for *ε* = 0.6. For *ε* = 0.95 a grazing collision takes place that ends up in two single drops of nearly the same sizes as those of the incident drops. Note that the relative number of collisions resulting in a specific breakup mode—which is identical to the mode-specific breakup efficiency (see the appendix) *E*_{b,i}, with *i* ∈ {*F*, *S*, *D*} indicating filament, sheet, and disk breakup, respectively—is given for all 32 drop pairs in Table 1.

### c. Budgets of kinetic energy and surface area

The temporal development of the kinetic energy of the fluid contained in the system and its surface area during a collision process is depicted exemplarily for pair 26, *ε* = 0.6, in Fig. 6. The respective images of the collision process are presented in Fig. 7 as a time series. For convenience, the initial kinetic energy is plotted starting from zero without loss of generality.

Looking at the surface area, there is at first a sudden drop as the drops coalesce. This surface reduction is caused by the contact of the two drops and the merging of their surfaces. Subsequently, because of the equilibration of momentum and mass a strong deformation of the generated larger drop continues, with the surface area increasing by about 50%. That is followed by the separation into one large drop and several smaller fragment droplets, where the drops attain a roughly spherical shape. Note the undulations in the surface area curve, which are a clear indication of drop oscillations.

The kinetic energy decreases at first during the merging phase of the drops. Then it increases again and reaches a maximum as the liquid disintegrates whereby it is transformed to surface energy. Subsequently, the kinetic energy is continuously reduced, presumably because of viscous dissipation due to the internal flow inside the drops and frictional interaction with the surrounding gas.

Next, for drop pair 9 only the variation of surface area as a function of time is examined for different eccentricities *ε* since the appearance of the collision strongly depends on this parameter. In Fig. 8 the surface area *A*, consisting of contributions of all drops occurring [*A* is calculated by *A* = |**∇***f*|, with *f* defined by Eq. (3); see Rieber 2004 for more details], divided by that of the spherical coalesced system *A*_{coalesced} = *π*(*d _{L}*

^{3}+

*d*

_{S}^{3})

^{2/3}is shown for the indicated eccentricities as a function of time. Initially all curves show slight minima followed by a strong increase, attaining distinct maxima and decreasing gently thereafter. It is clearly seen that the amplitude of this area ratio is strongly connected to the collision eccentricity. For

*ε*= 0.95 (grazing collision) there is almost no fluctuation. Going to more centric collisions, the minima—caused by the first contact of the drops—and the subsequent maxima, as collision energy gets transformed to surface energy, are more and more pronounced. Note that the area maxima reached for

*ε*= 0.2 and

*ε*= 0.05 are approximately the same. All area ratios subsequently decrease again and relax to a stable level as the oscillations of the droplet fragments decay.

### d. Coalescence efficiency and eccentricity effects

The coalescence efficiency *E _{c}* measures the probability for a binary collision process to end up in a single permanent drop. Notice that reference to a specific drop pair

*k*is omitted for

*E*and the breakup efficiency

_{c}*E*throughout the text. Here

_{b}*E*is obtained by dividing the number of collisions leading to a permanent coalescence,

_{c}*c*, by the total number of collisions

_{c}*c*:

In calculating *E _{c}* we have to note that

*c*depends on eccentricity, as has been discussed. However, most detailed cloud microphysical simulations do not (and may cannot) take into account such a dependency because of numerical limitations. Therefore, the eccentricity-dependent

_{c}*E*values are averaged here as follows: because of the model setup employed, each collision is representative for all collisions taking place in the according annulus of the projected collision cross section (circle with radius

_{c}*r*+

_{L}*r*). The result for a specific eccentricity

_{S}*ε*is weighted with the area of the corresponding annulus (

*δ*= separation distance, cf. Fig. 2):

Hence the coalescence efficiency is calculated via

where the subscript *c* indicates only cases where coalescence occurs. The variable *ε _{b}* describes

*ε*values for which breakup appears. Some relations concerning the sets of eccentricities and definitions of coalescence efficiencies are given in the appendix.

The outlined procedure of calculating *E _{c}* thus assumes an equipartition of the eccentricity-dependent

*E*values over the respective collision cross sections. This procedure is a first guess as no other information on preferred impact regions is available. Equation (16) is numerically computed using Eq. (2) of Part II considering the discrete

_{c}*ε*values mentioned in the text below Eq. (7).

Figure 9 shows the coalescence efficiencies as a function of the Weber number with the present results and those of LL82a. The following discussion is restricted to the data only. A comparison between the *E _{c}* parameterizations derived from the present results and those given by LL82b is left to our companion paper.

In a gross view of the figure, one recognizes some similarities and some discrepancies between the two datasets. The most remarkable difference appears for small We numbers. Whereas LL82a did not observe *E _{c}* values larger than about 0.65, the present data provide in some cases

*E*values larger than 0.8. Moreover, in two cases (pairs 11 and 30)

_{c}*E*values of 1.0 appear. The occurrence of

_{c}*E*= 1 can have two reasons. On the one hand, it may show up if a very small drop collides with a large drop (case 1) such that the energy imparted by the small drop is too small for breakup of the coalesced system. On the other hand,

_{c}*E*= 1 may result if the drops are of comparable size (case 2) so that only a smooth contact is made with low relative velocity, preventing breakup. But in the latter case it is not clear (and cannot be resolved by the present simulations) whether bouncing may take place.

_{c}Inspecting the data of Table 1, case 1 is true for pairs 2, 3, and 17. For these pairs the calculations deliver *E _{c}* values of about 0.8. That a value of exactly 1.0 is not attained can be explained by the fact that the collisions of these drop pairs do not lead to a permanent unification under all circumstances but rather result in filament breakup for the grazing case of

*ε*= 0.95 (cf. Table 1). Case 2 refers to pairs 5, 11, 16, 23, 30, and 32 having size ratios

*d*/

_{L}*d*smaller than about 1.8. Here we find the same situation as in case 1: whereas drop pairs 11 and 30 do not show any breakup, the remaining ones do, so such that

_{S}*E*becomes relatively small.

_{c}It should be clear already at this stage that an explicit consideration of eccentricity, which is not included in the parameters defined above as CKE, Δ*E _{T}*, and We, is important in determining coalescence efficiencies.

After having discussed *E _{c}* values for small Weber numbers, we now consider large Weber numbers where

*E*approaches 0. The condition of

_{c}*E*= 0 is equivalent to the situation that for a specific drop pair no collision ends up in a permanent unification of both drops or, conversely, that breakup takes places under all circumstances. The simulations show that for small eccentricities (i.e., nearly head-on collisions), often a single drop results. For larger eccentricities mostly breakup appears. Examining cases with small eccentricities,

_{c}*E*= 0 occurs if CKE is very large. This is valid for pairs 9, 10, and 31 (cf. the three circles in Fig. 9 for

_{c}*E*= 0). That this

_{c}*E*= 0 condition does not always follow from the reasoning just presented can be seen for drop pairs 18, 28, and 29, which also have large CKE values. For these drop pairs

_{c}*E*values slightly larger than 0 are calculated (cf. Table 1).

_{c}For intermediate We numbers there is a more or less good agreement between the present calculations and the LL82a data. Because of the many more drop pairs considered in this study as compared to LL82a, more data points arise. In this range of We numbers all collision outcomes depend on eccentricity.

Finally, it is emphasized that the present calculations do not confirm the cutoff value of *E _{c}* with regard to CKE as it was given by LL82a and recently attested to by Barros et al. (2008).

In summary, we conclude that eccentricity has a distinct effect on coalescence efficiencies. We emphasize that this finding is the result of the application of a numerical simulation in which all conditions relevant to collision-induced breakup can be fixed and the results can be reproduced. This is different from situations where laboratory experiments are performed and reproducible eccentricity effects may only be investigated with great difficulty.

To elucidate the eccentricity effects more closely, the eccentricity-specific number of outcomes (contours; for definition see the appendix) and the breakup type (different symbols) for all 32 investigated drop pairs and all considered eccentricities are presented in Fig. 10.

As already shown by means of two exemplary drop pairs (cf. Figs. 4 and 5), disk breakup with numerous fragments occurs for almost centric collisions (*ε* = 0.05 and *ε* = 0.2) and the highest collision energies, whereas for smaller CKE most drops coalesce. As *ε* shifts to higher values, the breakup zone expands at the expense of coalescences and the breakup type changes to sheet and subsequently to filament breakup. At the same time, the number of fragments decreases. For values of *ε* approaching unity there is breakup for almost all drop pairs and the number of fragments is further reduced, approaching the persistence of the two initial drops for almost grazing collisions.

As mentioned in the introduction, some efforts have been undertaken in the past to calculate a criterion separating the range where coalescence occurs from that for breakup. The criterion considered here is the one of Brazier-Smith et al. (1972) giving a critical eccentricity *ε*_{crit} as

with *q* = *d _{S}*/

*d*and We* defined by Eq. (13). This relation assumes that the rotational kinetic energy has to exceed the surface energy for breakup. Details are given in Part II. Isolines for

_{L}*ε*

_{crit}≥ 0.4 are additionally plotted in Fig. 10 for all considered eccentricities.

One can observe a reasonable agreement, particularly for *ε* = 0.4 − 0.8. For *ε* = 0.95 there are several simulations predicting separation for Weber numbers smaller than those given by Eq. (17). However, Brazier-Smith et al. (1972) state that this region is hard to capture and Eq. (17) only discriminates between separation with fragment droplet production and coalescence whereas there are no fragment droplets produced in the simulations for the grazing collisions at *ε* = 0.95.

In weighting the data in Fig. 10 by the eccentricity area [i.e., constructing *f*_{k} (cf. the appendix)], this variable is presented in Fig. 11. The symbols indicate the most probable breakup mode (cf. Table 1) and the color contours the average number of outcomes (coalesced and fragment drops) produced. One can see that the most prominent breakup mode is filament, with total fragment numbers ranging between two and about eight. The largest fragment numbers (about nine) result from sheet breakup occurring if *d _{L}* ≥ 0.4 cm and

*d*≥ 0.14 cm. Disk breakup has disappeared since its probability is highest only for small eccentricity values that contribute only little to the average. Note that the region with largest fragment numbers coincides with that of largest CKE values (not shown).

_{S}In Fig. 12, indication of coalescence and breakup events in a We–*ε* plane can be seen. Included in this figure are three curves. The solid one is an approximate fit given by *ε*_{crit} = exp(−0.65We) separating the coalescence region from the sheet and filament breakup region. The dashed–dotted one is a schematic separating the coalescence region from the disk breakup region for very small eccentricities. The dashed one indicates approximately the boundary between the filament and sheet breakup modes given by We = 46.36*ε*^{2} − 51.06*ε* + 15.4. In comparing this figure to the respective ones of Ashgriz and Poo (1990), for instance, a striking similarity appears. In their nomenclature the upper and lower breakup regions are signified as those for stretching and reflexive separation, respectively. However, there is a strong difference as well. Whereas in the figure the Weber number is calculated by Eq. (11), the Weber number We* used by Ashgriz and Poo (1990) is defined by Eq. (12); both Weber numbers are interrelated by Eq. (13). Plotting the present data in terms of We* shows an erratic picture (not presented here). Also, other limit *ε*_{crit} conditions [we refer to Menchaca-Rocha et al. (1997) for a compilation of those conditions] are not appropriate to explain the behavior exhibited in Fig. 12. About the reason for this discrepancy one can only speculate. Perhaps the discrepancy has its roots in assumptions underlying the derivation of the other limit *ε*_{crit} conditions, such as suggestions on neck dynamics, inclusion of rotational energy effects, and presuming that after a collision event only two fragment drops are created, which for the process here considered seem not to be applicable.

Consequently, it can be argued that We* is not the appropriate number to include all parameters relevant to the collision process under consideration. In particular, it is hard to accept that the diameter of the smaller drop alone, as contained in the definition of We*, is decisive for the outcome of a collision event with a larger drop.

### e. Number of outcomes and fragment size distributions

According to LL82b the total number of collision outcomes for a specific drop pair *k* is given by

Here *f*_{k,b} is the number of fragments due to breakup and *f*_{k,c} = 1 considers the formation of a single permanently coalesced drop.

The variables appearing here can be deduced from the results of the numerical simulations. The information used comprises the sets of eccentricities for permanent coalescence and for breakup as well as the respective spectral drop size distributions: *f*_{k,b}(*ε _{b}*,

*D*) for the fragments and

*f*

_{k,c}(

*ε*,

_{c}*D*) for the coalesced drops, with

*D*= diameter of coalesced or fragment drop.

In characterizing coalescence and breakup we present the total number of outcomes weighted by the eccentricity area, *f*_{k}, the eccentricity-specific number of outcomes, , and the spectral number of fragments, *f̃*_{k,b}(*D*). For details on the definition of these and other variables we refer to the appendix.

Note that all variables depicted in the following figures are evaluated numerically from the respective relations; examples can be found in Part II.

The total number of outcomes *f*_{k} obtained for all investigated drop pairs is depicted in Fig. 13 as a function of the collision energy together with data measured by LL82a. The results are very similar in both cases, especially if one considers the standard deviation of the LL82a data. This means that on average (i.e., regarding all eccentricities used in the present study), LL82a observed nearly the same total number of fragments. An approximate fit relating *f*_{k} to CKE given by *f*_{k} = (1.5 − CKE^{0.135})^{−1} for 10^{−2}*μ*J ≤ CKE ≤ 10 *μ*J is additionally included in the figure.

### f. Further influence of eccentricity

The number of outcomes versus eccentricity is plotted exemplarily in Fig. 14 for three pairs of drops with different CKE.

As already stated, there is a strong influence of eccentricity on the number of fragment droplets produced. However, this statement is not unambiguous. Looking at the high collision energy case, a disk mode breakup occurs for low eccentricities and many small fragments are produced. Increasing eccentricity leads to a decreasing number of fragment droplets until the initial drops persist for a grazing collision. On the other hand, the behavior in case of the medium and low collision energy cases is contrary. While the two drops coalesce for almost central collisions and persist for grazing collisions, a number of fragment droplets are produced in between.

The numerically calculated fragment size distributions in spectral resolution *f̃*_{k,b}(*D*) are shown as histograms for three examples in Fig. 15 together with the parameterization curves of LL82b. The fragment size distributions of all investigated drop pairs can be found in Part II.

As can be seen in the figure, the LL82b parameterization fits very well for pair 5 (upper plate). For pair 10 (middle plate) the agreement is qualitatively satisfactory but fewer relatively small fragments appear in the numerical calculations. For pair 26 (lower plate) major discrepancies appear concerning the position and number of relatively small fragments. The agreement between the numerical and experimental results is good for emerging larger fragments, whereas large deviations are observed especially for very small fragments. On the one hand, this finding results from the ability of the numerical method applied here to resolve small fragments, in contrast to the LL82a parameterization employing a lower limit of small fragments of *D* = 0.05 cm. However, note that the size of smallest fragment drops calculated in this study is limited by the grid resolution, which in the standard calculations is *D* = 0.01 cm. On the other hand, the prediction of small droplet fragments shows some grid dependency, as presented in the following section. In addition, the experimental results may have been influenced by systematic errors, which is particularly important if one takes into account the strong dependency on eccentricity, as shown earlier.

### g. Assessment of results

Assessing the obtained numerical results, various factors have to be considered. The most relevant ones seem to be both the spatial resolution and the boundary conditions.

To elucidate the grid dependency of the results, particular cases have been simulated with doubled number of computational cells in each spatial direction (i.e., halved spatial resolution), resulting in ≈135 × 10^{6} cells. The results for the standard resolution (Δ*x* = 100 *μ*m) and the fine resolution (Δ*x* = 50 *μ*m) are shown in Fig. 16 for pair 9, *ε* = 0.5 and *ε* = 0.8. It is well known that the ability of the volume of fluid method, which is employed for tracking the interface, to correctly reproduce the physical behavior of small ligaments is quite sensitive to grid resolution. This fact is reflected by the presented results. For both eccentricities, there are more small fragments in the case of finer grid resolution. For larger fragments, however, it is important to point out that there is only a minor discrepancy between the results applying either a fine or a coarse grid. From a physical point of view it seems feasible to apply a lower radius limit of fragment droplets to be considered as very small droplets may under subsaturated conditions evaporate quickly and thus disappear. This situation is different in a saturated environment.

As to the boundary conditions, the drawback of the presented approach with periodic boundary conditions is that the drops are falling in their own wake, which is of major importance to large drops where the aerodynamic influence of the oncoming flow plays an important role for the deformation and the breakup process.

One may be inclined to refine the spatial resolution further. However, it should be mentioned that for the standard resolution the numerical effort varies between several hours and a couple of days for one single simulation (i.e., one eccentricity of one drop pair!) using up to four processors of an NEC SX-8 machine. The variation is due to different computational domains needed to provide a sufficient distance of the drops to the boundaries and different real time durations of the collision process. Thus, a further significant refinement of the numerical grid, though desirable, seems presently out of reach.

## 5. Summary

The collision and coalescence/breakup process of different pairs of raindrops has been investigated using FS3D, an advanced DNS-CFD code dealing with free surfaces. An important ingredient in the present investigation is the consideration of eccentricity for which six different values have been chosen. In total 32 drop pairs have been studied; among them are also those investigated by LL82a. By the choice of the drop pairs the parameter space, spanned by both individual drop diameters relevant to collision-induced breakup, is covered nearly completely.

From the numerical experiments coalescence efficiencies are derived which for small Weber numbers show a significant difference to those obtained by LL82a. An empirical fit for the boundary separating the coalescence and breakup regions is given. This criterion differs from criterions determined by other experimental and theoretical investigations. Moreover, a certain value of the collision kinetic energy limiting coalescence as found by LL82a could not be confirmed.

Eccentricity plays a crucial role for the appearance of a specific breakup mode (disk, sheet, and filament) and for the number of fragment droplets created. As could be expected, the highest values of collision kinetic energy result in a relatively large number of fragments. The spectral size distributions of the fragment droplets show an uneven picture. Compared to the parameterizations derived by LL82a, the spectral fragment numbers obtained in this study exhibit some similarities but also some discrepancies. Because of the numerical method applied, very small fragment droplets, which are difficult to capture by laboratory experiments, could also be detected.

The results presented here are the basis for parameterizations of the coalescence efficiency and the fragment size distributions shown in Part II, together with computations of stationary drop spectra resulting from application of the new parameterizations.

## Acknowledgments

The simulations were performed on the national supercomputer NEC SX-8 at the High Performance Computing Center Stuttgart (HLRS) under Grant FS3D/11142. The authors gratefully acknowledge support by the German Research Foundation (DFG) under Grants BE 2081/7-1 and WE 2549/17-1 in the framework of the priority program Quantitative Precipitation Forecast. Several useful comments by the reviewers helped to considerably improve this study. The authors are also thankful to Dr. C. Ancun for a very stimulating discussion.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Coalescence Efficiencies and Specific and Total Numbers of Collision Outcomes: General Definitions and Relations

In this appendix some general definitions and relations are presented concerning the calculation of coalescence efficiencies as well as of specific and total numbers of collision outcomes.

We start with considerations on eccentricity values. The set of discrete eccentricities chosen in this study is

The cardinal number of this set is obviously *N _{ε}* = 6. Note that in general 0 ≤

*ε*≤ 1 with arbitrary intermediate values. In the set chosen here eccentricities denoting permanent coalescences and breakup events are indicated in the following by the subscript

*c*and

*b*, respectively. The numerical results provide the set of coalescences,

*S*, containing all

_{εc}*ε*values, and the set of breakup events,

_{c}*S*, containing all

_{εb}*ε*values, with

_{b}*S*∪

_{εc}*S*=

_{εb}*S*and

_{ε}*S*∩

_{εc}*S*= {}. The latter relation shows that both sets are disjoint. For the cardinal numbers it yields

_{εb}*N*= 6 =

_{ε}*N*+

_{εc}*N*.

_{εb}##### Coalescence efficiencies

The coalescence efficiency is defined as the ratio of the number of collisions leading to permanent coalescence, *c _{c}*, to the total number of collisions

*c*[cf. the identical Eq. (14)]:

The coalescence efficiency for breakup, which one may call breakup efficiency, is complementary to the number of collisions resulting in breakup, *c _{b}*

such that evidently *E _{c}* +

*E*= 1. In considering the specific breakup modes (filament, sheet, and disk; cf. section 4b), it yields

_{b}*E*=

_{b}*E*

_{b,F}+

*E*

_{b,S}+

*E*

_{b,D}, where the single contributions are given in Table 1.

Weighting *E _{c}* by the respective eccentricity area gives [cf. the identical Eq. (16)]

Equivalently, it follows that

so that finally

##### Specific and total numbers of collision outcomes

The numerical calculations deliver for all *ε _{b}* values the fragment size distributions referring to a specific drop pair

*k*,

*f*

_{k,b}(

*ε*,

_{b}*D*) with

*D*= fragment diameter, as well as for all

*ε*values the singular size distribution of the coalesced drop,

_{c}*f*

_{k,c}(

*ε*,

_{c}*D*) =

*δ*(

*D*−

*D*) with

_{c}*D*= (

_{c}*d*

_{S}^{3}+

*d*

_{L}^{3})

^{1/3}= diameter of coalesced drop and

*δ*= Dirac delta function. Clearly the joint size distribution is given by

*f*(

_{k}*ε*,

*D*) =

*f*

_{k,b}(

*ε*,

_{b}*D*) +

*f*

_{k,c}(

*ε*,

_{c}*D*). With this size distribution and its contributions by breakup and coalescence, respectively, the integral relations presented in Tables A1 and A2 can be calculated. In Table A2 the integrals are given with explicit eccentricity weighting as applied in this study whereas in Table A1 the equivalent integrals are listed irrespective of eccentricity weighting for completeness. Note that only some variables of Table A2 are considered in this study.

###### Specific and total numbers irrespective to eccentricity weighting

According to this compilation we have ; that is, the number of outcomes of a collision event *k* is the sum of the number of coalesced drops and that of the total number of fragment drops.

###### Specific and total numbers with explicit eccentricity weighting

According to this list we have

from which it follows that

which is the same relation as from LL82b for the total number of coalesced drops and drop fragments.

Equivalent to Eq. (A8), the spectral number of outcomes is

which is a relation referred to in Part II.

## Footnotes

* Current affiliation: North Rhine-Westphalia State Agency for Nature, Environment and Consumer Protection, Leibnizstr. 10, 45659 Recklinghausen, Germany

*Corresponding author address:* Dr. Klaus D. Beheng, Institut für Meteorologie und Klimaforschung, Karlsruher Institut für Technologie, 76131 Karlsruhe, Germany. Email: klaus.beheng@kit.edu

^{1}

Note that Barros et al. (2008) do not propose a parameterization based on their observations.