Keywords: suspensions, collisions, computational uid

is thesis is devoted to the mechanisms leading to strong collision rates of inertial particles in turbulent suspensions. Our work is based on simulating the motion of particles, using both direct numerical simulations of the Navier-Stokes equations, and a simpler model (kinematic simulations). is subject is important for many applications, in industrial as well as natural (astrophysical, geophysical) contexts. We revisit the ghost collision approximation (GCA), widely used to determine the rate of collisions in numerical simulations, which consists in counting how many times the centers of two particles come within a given distance. eoretical arguments suggested that this approximation leads to an overestimate of the real collision rate.

is work provides not only a quantitative description of this overestimate, but also a detailed understanding of the error made using the GCA. We nd that a given particle pair may undergo multiple collisions with a relatively high probability. is is related to the observation that in turbulent ows, particle pairs may stay close for a very long time. We have provided a full quantitative characterization of the time spent together by pairs of particles.

A second class of results obtained in this thesis concerns a quantitative understanding of the very strong collision rates o en observed. We demonstrate that when the particle inertia is not very small, the "sling/caustics " e ect, i.e., the ejection of particles from energetic vortices in the ow, is responsible for the high collision rates. e preferential concentration of particles in some regions of space plays in comparison a weaker role. turbulence, inertial particles, turbulent

collision. Cette thèse fournit non seulement une estimation quantitative de cette surestimation, mais également une compréhension détaillée des mécanismes des erreurs faites par l'ACF. Nous trouvons qu'une paire de particules peut subir des collisions répétées avec une grande probabilité. Ceci est relié à l' observation que, dans un écoulement turbulent, certaines paires de particules peuvent rester proches pendant très longtemps. Une deuxième classe de résultats obtenus dans cette thèse a permis une compréhension quantitative des très forts taux de collisions souvent observés. Nous montrons que lorsque l'inertie des particules n' est pas très petite, l' e et « fronde/caustiques », à savoir, l' éjection de particules par des tourbillons intenses, est responsable du taux de collision élevé. En comparaison, la concentration préférentielle de particules dans certaines régions de l' espace joue un rôle mineur.

turbulence, particules inertielles, suspensions turbulentes, collisions, simulations numériques d' écoulements turbulents, écoulements multi-phasiques viii PUBLICATIONS Some of the ideas and gures presented in this thesis have been published before in the following papers:

Voßkuhle M,

INTRODUCTION

Among the unsolved problems in physics, turbulence is probably the one that a ects our daily life the most. e search for the ultimate components of matter, or for the structure of the universe leads to fundamental questions. However, our lack of understanding about dark matter or dark energy comes in large part from the scarcity of experimental information, which ultimately originates from the fact that they are relevant at scales which are very remote from our everyday life. is situation makes the success of the standard model even more impressive-despite the fact that it may still be improved, or that it does not include a theory of gravitation. In comparison, turbulence is encountered in everyday life, yet, we do not have a fully satisfactory theory to describe this important phenomenon.

If we ride our bicycle or Vélo'v along the Rhône, we know how di cult it e Vélo'v is a rental bike in Lyon.

is to gain that extra km h, that would allow us to be in time. e reason is the air, which we bring into motion by moving through it-just like the car in Figure . . Simple dimensional arguments allow us to understand, why it is so hard to be faster. Riding our bike at speed v and for a time t, we leave behind us a wake of air in motion of volume Avt, where A is roughly the cross-sectional area of us facing the direction we are heading to. To estimate the energy contained in this corridor we need to estimate its mass, which is simply Avtρ, where ρ is the density of air. is mass is moving at speed v and therefore contains the energy ⁄ Avtρv . is means the energy we loose just to agitate the gas surrounding us grows with v . Going at km h instead of km h is two times harder. Furthermore the linear prefactor A Looking closely at the wake behind the car in Figure . , one notices that the motion of the air is not structureless. While the motion seems chaotic and unpredictable on the one hand, one can on the other hand clearly discern structures on many length scales. is is one of the intriguing features of turbulence. A detailed analysis shows that the above estimate for the energy contained in the wake is too simple. In fact the surface A needs to be replaced by a sort of e ective surface described by the so-called drag coe cient. is drag coe cient depends on the geometry of the object that moves through the uid. Tiny changes in the design of vehicles may easily lead to changes in e ciency of the order of . erefore car producers a very interested in the drag coe cient of new models they develop. Understanding the turbulence in the wake is crucial in estimating the drag coe cient without costly wind tunnel experiments.

e above is just one example where a complete and fully predictive theory of turbulence is missing. ere are many other domains including weather prediction or various industrial processes, which would bene t from such a theory. But still this phenomenon evades a conclusive description. In a review on turbulence John von Neumann expressed "hope to at review has been reprinted in von Neumann ( ).

well-known Millennium Simulation (Springel ; Springel et al.

) considered interactions between "only" × particles. And still, the simulation of Kaneda et al. (

) achieves a ow at a Reynolds number¹ Re λ = which is relatively small in comparison to Reynolds numbers observed in many natural phenomena (e.g., Re λ = in a typical cloud, Siebert et al.

). e reason why it is so hard to perform simulations of turbulent ows at large Reynolds numbers are again the structures, which are present on a wide range of length scales. An exact simulation must resolve all these length scales which may span over several decades. e present study is in a sense in the tradition of the above cited works.

Here also numerical simulations are performed to study turbulent ows. But instead of seeking insights into the fundamental structure of turbulence, this work concentrates on the phenomena that arise when particles are entrained by turbulent ows.

. In nature, more o en than not, uids do transport material in another state of matter (Balachandar & Eaton

). One may think of aerosol particles in the atmosphere, plankton in the sea, or pollutants in riverbeds. Evidently, whether the uid is turbulent or not has a big impact on the particulate phase. Stirring in the ocean enhances the mixing of phytoplankton and fertilizers (Abraham et al.

). Also, animals, that rely on olfactory perception, may need to adapt their strategies in a turbulent environment (Ferner & Weissburg ; Mafra-Neto & Carde ). As long as the entrained objects are small enough, they follow the movements of the ow exactly. For larger particles however, with a di erent density than the uid, new phenomena emerge. Clustering of particles in certain regions of the ow is a well-known example.

In this work we focus in particular on collisions between the entrained objects. e collisions occurring between particles are important in a wide range of applications. Many industrial production processes, like manufacturing of titania (Moody & Collins

), rely on turbulent mixing to bring the base products into contact. In all kinds of sprays, the coalescence rate of the droplets has an impact on their size distribution. is is of importance for example in diesel engines (Post & Abraham ). ere, further aspects, like the evaporation of droplets, provide motivation for interesting research (Chareyron et al. ). Also in nature formation of larger particles by collision is very frequent. Wells & Goldberg ( ) nd, that there is an abundance of colloidal particles in sea water. For instance a signi cant part of the carbon, that is considered dissolved in the oceans, "may in fact be present in the form of colloid particles" (Stumm & Morgan , p. )². Another example for the importance of collisions in the dispersed phase of a uid is provided by sandstorms, whose intensity can be enhanced by colliding grains (Carneiro et al. ). Even animate matter depends on turbulence induced collisions as can be seen in the works of Lewis & Pedley (

) and Rothschild & Osborn (

). Our main motivation (discussed in Chapter ) is provided by two further domains. On the one hand, collisions between droplets in clouds, contribute to their growth process and therefore in uence the onset of rain (Shaw ). ey also a ect the size distribution of drops, which has a strong in uence on the re ection of light from the sun, and thus, on the energy budget of the earth and on the climate. On the other hand we refer to collisions of dust grains in protoplanetary disks that play a crucial role in the formation of planets (Safranov ). e contribution of turbulence to the collision rates has been the subject of many studies, in relation to the two problems referred to here.

e aim of this work is to provide further quantitative understanding of the problem.

. is work is organized in three parts. e rst part includes this introduction as well as two chapters on fundamental concepts. Chapter introduces the basic equations and some of the typical reasoning that will be used throughout the text. e following chapter provides a short review of cloud physics in Section . . In Section . di erent aspects of collision rates of particles suspended in turbulent ows are discussed. is section constitutes the basics which are important to understand our scienti c work presented in the second part.

We split our own results in ve chapters. e rst two, Chapters and , describe di erent aspects of our numerical simulations. Furthermore, we validate our code, by comparing its results to published data and theoretical predictions. Additionally we present some simple results which we obtained theoretically for ideal gas-like particles in Sections . and . .

Especially the ndings from Section . will be helpful in interpreting some of the results presented in Chapter . ere we investigate multiple collisions between a same pair of particles, which astonishingly may take place in a turbulent ow (Section . ). ese multiple collisions are related to the ghost collision approximation (GCA) in Section . . at approximation is o en made in numerical simulations but is also apparent in theoretical results. Inspired by the observation of multiple collisions we proceed to investigate the time that particles stay in contact in Section . . Before discussing an alternative to the GCA in Section . , we present an interpretation of our results in the light of sling/caustics collisions in Section . . Section . concludes that chapter.

In Chapter we come back to the sling/caustics e ect. A subtle postprocessing of our numerical data allows us to conclude on the prevalence of the sling/caustics e ect over simple shear induced collisions. ese results are presented in Section . and con rmed by additional investigations in Section . . e analysis performed in this work is in fact fully consistent with evidence that was already present in previous numerical results, or anticipated theoretically using simpli ed models, as we explain in Section . .

All our results presented until Chapter were obtained in direct numerical simulations (DNS). But we started our investigations using a simple model ow known as kinematic simulations (KS). In Chapter we discuss how the lower computational cost of this model allows some further investigations which would be prohibitive in DNS. e qualitative results in KS are the same as presented before for DNS. But quantitatively we nd large di erences which are shortly discussed in that chapter.

Finally, in the Appendices, we will present two more technical aspects of our work. Some of the approaches developed for our investigations and discussed in Appendices A and B, may be interesting for future work.

FUNDAMENTAL RESULTS

A er an introduction of the subject from a broader perspective, and a general presentation of the reasons why it is worthwhile to understand the behavior of particles in turbulent ows, we now discuss the relevant results our work relies on. is presentation does not have the ambition to be exhaustive. Details not covered here can be found in the literature cited in the text. Expert readers may skip this chapter and only get back to it, when it is referenced later on in the text.

.

-¹

To derive the basic equations governing the motion of a uid, one starts from two balance equations. e rst one describes mass conservation in terms of the density eld ρ(x, t), and the other one conservation of momentum, in terms of ρ(x, t)u(x, t), where u(x, t) is the uid velocity at position x and time t. For an incompressible uid, the density is independent of time and position and the continuity equation derived from mass conservation simpli es to ∇ ⋅ u(x, t) = .

( . )

e equation for the momentum has to accurately account for all forces acting on the uid's surface or volume. In case of incompressible Newtonian uids, one obtains

∂ ∂t u(x, t) + u(x, t) ⋅ ∇u(x, t) = -∇ (x, t) + ν∆u(x, t) + Q(x, t), ( . )
where the momentum equation has been divided by the constant uid density ρ. erefore (x, t) is the pressure divided by the density, but usually referred to as just "pressure". e kinematic viscosity ν = µ ρ is the quotient of dynamic viscosity and density, and nally Q(x, t) contains all di erent external forces like gravitation or the Coriolis force. Equation ( . ) has rst been derived by Navier (

) and Stokes ( ) and hence is known as the Navier-Stokes equation.²

Even ignoring external forces, Equation ( . ) may seem to be unclosed at rst sight, because the pressure (x, t) appears in it. But actually the pressure is fully determined by Equations ( . ) and ( . ) as can be seen when one calculates the divergence of the Navier-Stokes equation. In case of divergence free forces Q(x, t) this yields the Poisson equation

∆ (x, t) = -∇ ⋅ [(u(x, t) ⋅ ∇) u(x, t)] ,
is section and parts of the next one loosely follow a former work, originally written in German (Voßkuhle ). Or more precisely, as they are actually three di erential equations for each component of u(x, t) and need to be combined with Equation ( . ): the Navier-Stokes equations.

which can be solved with the help of Green's function for the Laplacian (see Friedrich

) or by inversion of the Laplacian in Fourier space (see Frisch ). Furthermore the solution of Equations ( . ) and ( . ) requires the determination of initial and boundary conditions. In realistic problems the boundary conditions are usually given by impermeability and no-slip walls. In many computer simulations however, it is convenient to impose periodic boundary conditions in a cube. is method will also be employed throughout this work.

. . he vorticity

For the sake of completeness, we introduce the vorticity ω(x, t) which is simply the rotation of the velocity eld

ω(x, t) = ∇ × u(x, t). ( . )
In the case of a solid body motion, this de nition of ω(x, t) is, within a factor of , the de nition of the rate of rotation. e equation describing the evolution of the vorticity ω(x, t) is deduced from Equation ( . ) by applying the rotation operator.

A description of turbulent uids in terms of "eddies", or vortices, with a wide distribution of scales has emerged over the years, starting with the famous illustrations of Leonardo da Vinci (see Figure . ). In fact, the vorticity eld is as fundamental as the velocity eld, as the latter can be obtained from the former by simply inverting Poisson's equations. Mo att (

) and Sa man ( ) review essential results in this eld.

. . elf similarity and the eynolds number

An adimensional form of the Navier-Stokes equation can be handily derived by introducing the dimensionless units (see, e.g., Argyris et al.

)

ũ = u U , t = t T , x = x L , p = p L U , Q = Q T U ,
with some characteristic length, time, and velocity scales L, T, and U = L T.

Inserting in Equation ( . ) and immediately dropping the tildas, one obtains

∂ ∂t u(x, t)+u(x, t)⋅∇u(x, t) = -∇ (x, t)+ Re ∆u(x, t)+ Q(x, t), ( . )
where we introduced the dimensionless Reynolds number Re = U L ν. It is named a er Osborne Reynolds ( ), who noticed in experiments on pipe ows, that whether the ow is turbulent or laminar depends only on the quantity Re. Indeed, Equation ( . ) demonstrates, that ows with the same Reynolds number behave similarly. is allows engineers to test new developments on model vehicles with reduced size in wind tunnels.

. e observation Reynolds (

) made in his experiments was the following. He was pumping a given uid, at a xed viscosity ν, through a cylindrical tube of diameter L. By increasing the pumping velocity, hence U, he could increase the value of Re. At low values of Re, the uid motion is observed to be laminar. Tiny tracer particles would along straight lines. Above a critical value of Re, however, Reynolds observed that the tracer trajectories were no longer straight. e uid would move in an irregular and unpredictable way. is state of uid motion is known as turbulence. Its emergence may become more evident, when one interprets the Reynolds number as the ratio between the nonlinear convection term u(x, t) ⋅ ∇u(x, t) and the dissipation term of the Navier-Stokes equations. When the smoothing action of dissipation is outdone by convection, the ow becomes turbulent.

We explained in the introduction (Chapter ) why it is of such an importance to gain a deeper understanding of turbulence. And it was also discussed, that a satisfactory theory of the phenomenon has not yet been proposed. Nevertheless, during its long history, research in turbulence has produced many important results. Nowadays there is a vast number of textbooks on it, from the classic ones by Batchelor ( ) or Tennekes & Lumley (

) to more modern ones like those from Frisch ( ) or Pope (

). Each of these books focuses on some speci c aspects of the problem. Here we will only discuss some of the most basic notions, that can be found in almost all textbooks. e typical reasoning we introduce in the coming paragraphs will reappear throughout the text.

. . he cascade picture

e observation that turbulent ows involve many length scales is already obvious in some of Leonardo da Vinci's writings and drawings. Figure .  shows clearly how the large scale structures at which the turbulence is generated -in this case water owing in a basin-evolve to produce structures on many length scales. A phenomenological explanation of this fact is attributed to Richardson ( ), who in his famous adaption of a poem by Jonathan Swi , explains his view e realize thus that: big whirls have little whirls that feed on their velocity, and little whirls have lesser whirls and so on to viscosity-in the molecular sense.

What he means is that turbulence receives its energy from large scale vortices (or coherent structures) generated by some external forcing. ese vortices will then break up into smaller vortices, who themselves get unstable and break up. e energy gets, in this way, transported to ever smaller scales until it reaches the smallest scales where viscosity acts and transforms it into heat.

e above described picture is known as the Richardson cascade; Figure . shows a typical illustration. We shall call the length scale at which the turbulence is generated L and the one, where the energy gets dissipated by viscosity η. e coming sections will show, how these heuristically de ned scales can get quanti ed.

. . omogenous and isotropic turbulence

It became clear quite early, that the assumption of turbulence being statistically homogenous and isotropic leads to signi cant simpli cations. Already William omson ( ), the later Lord Kelvin, supposed that the velocity eld of a turbulent uid was homogenous and posed the question, in what situation the assumption of isotropy would be justi able. Subsequent experimental studies supported the validity of both assumptions. Here we follow later work by Taylor ( ) to clarify the concept. An interesting account of the development of the scienti c ideas leading to our present knowledge can be found in a recent review of the correspondence between Prandtl and Taylor (Bodenschatz & Eckert

). We start by introducing the spatial velocity autocorrelation function

R i j (r, x, t) = ⟨u i (x + r, t)u j (x, y)⟩ . ( . )
Here, as well as later in the text, the angular brackets ⟨⋅⟩ denote a suitably de ned average. e velocity correlation function is a very interesting quantity; it is relatively easy to measure in experiments and numerical simulations, and also convenient to analyze from a mathematical point of view. Metaphorically speaking, it tells how much the velocity at point x is in uenced by the velocity at position x + r.

If we assume homogeneity, all points in space become statistically equivalent and the velocity correlation function does not depend on x any more. A parabola, corresponding to the second term in a Taylor expansion of the transverse velocity autocorrelation function, intersects the axis at the Taylor length λ, Equation ( . ).

For stationary statistics, the correlation function does not depend on time.

Finally, if the statistics are assumed to be isotropic³, the second order tensor R i j (r) can only take a certain form (see, e.g., Robertson )

R i j (r) = R tt (rr)δ i j + [R l l (r) -R tt (r)] r i r j r ( . )
is is discussed in more detail by, e.g., Batchelor ( ) and Landau & Lifshitz (

). R l l (r) and R tt (r) are respectively the longitudinal and transverse autocorrelation functions. Taylor ( ) used them to de ne two distinct length scales. Both may be de ned for the longitudinal as well as for the transverse correlation, but we limit ourselves to the latter.

First, we de ne a short length scale. To this end R tt (r) is developed into a series. e rst term of this series, R ′ tt ( ), vanishes due to incompressibility, Equation ( . ). e second term can be evaluated, R ′′ tt ( ) = ⟨(∂u ∂x ) ⟩ (Pope ), and used to de ne the length scale

λ = u rms ⟨ ∂u ∂x ⟩ . ( . )
A geometrical interpretation of λ can be found in Figure . . Taylor ( ) interpreted λ as the length scale of the smallest eddies. is is not correct, actually λ is intermediate between the smallest and largest turbulent length scales, η and L respectively (see, e.g., Pope

, for details). Despite this fact, λ is widely used and its simple but exact de nition is appreciated. e correspondingly de ned Taylor scale Reynolds number

Re λ = u rms λ ν ( . )
is one of two quasi-standard forms to give the Reynolds number in DNS and in experiments.

Next, we introduce the integral length scale L t . e index t does not only indicate that it is based on the transverse velocity autocorrelation function, but also di erentiates it from the conceptual integral length scale L, that has at means invariant under rotation and usually also re exion. e variable urms denotes the root mean square velocity of one component, i.e., urms = ⟨u(x, t) ⋅ u(x, t) ⟩ . been introduced before. It is a large length scale given by the surface under the autocorrelation function

L t = u rms ∫ ∞ R tt (r) dr.
( . )

e velocity eld decorrelates at L t , which therefore characterizes the length scale of the largest eddies.

As mentioned before, corresponding length scales can be derived from the longitudinal velocity correlation function R l l (r). It can be shown that those are not independent from λ and L t respectively. Indeed, exact relations can be derived as shown, e.g., in the textbook by Pope ( ).

. . nergy ux and energy dissipation rate

Before coming to the famous phenomenology of Kolmogorov, we need to introduce one more quantity characterizing the cascade process-the energy dissipation rate ε.

An equation describing the temporal evolution of the average energy per unit mass, ⟨u ⟩, can easily be derived from the Navier-Stokes equation.

To this end, one multiplies Equation ( . ) by u(x, t) and then averages the complete equation. By standard techniques, taking into account incompressibility, Equation ( . ), and especially the boundary conditions, one obtains

d dt ⟨u(x, t) ⟩ = ν⟨u(x, t) ⋅ ∆u(x, t)⟩ + ⟨u(x, t) ⋅ Q(x, t)⟩. ( . )
Frisch ( , p. ) explains the derivation for periodic boundary conditions, McComb ( , p.

) does the calculation in an arbitrary volume with vanishing velocity at the boundaries. In the stationary case the le -hand side of Equation ( . ) is zero and obviously energy production by the external forces and energy dissipation by viscosity balance. We may therefore de ne the energy dissipation rate per unit mass

ε = -ν⟨u(x, t) ⋅ ∆u(x, t)⟩ = ν i, j ⟨ ∂ i u j + ∂ j u i ⟩ = ν⟨ ω ⟩, ( . )
where the latter two equalities can be derived by similar techniques as before and ω(x, t) is the vorticity as de ned in Equation ( . ).

It is interesting to remark, that the nonlinear and nonlocal terms of the Navier-Stokes equation do not contribute to the energy budget in Equation ( . ). Although they neither produce, nor dissipate energy, they do contribute to the transport of energy between scales. is becomes clear when one considers the energy balance in Fourier space. Namely, one observes that

⟨u i (x, t)u i (x, t)⟩ = ∫∫∫ ∞ -∞ Rii (k, t) dk = ∫ ∞ E(k, t) dk,
We introduce the notation ∂i ∶= ∂ ∂xi for spatial derivatives.

.

k i F(k) ∼ k E(k) E(k) k Π(k) = ε F .
Illustration of the di erent terms contributing to the energy spectrum balance, Equation ( . ).

e forcing F(k) injects energy in a band of small wavenumbers; the energy ux Π(k) transports the energy to larger modes, where viscosity gains in uence and dissipates it.

where Einstein summation applies to repeated indices in the rst two terms. Ri j (k, t) is the Fourier transform of the velocity autocorrelation tensor R i j (r, t), de ned in Equation ( . ), with respect to r. e energy spectrum E(k, t) is a de nition that follows by simply integrating over the angular components (see, e.g., Tennekes & Lumley ). It can be interpreted as giving the energy contained in each scale k.

An evolution equation for E(k) can be obtained from the Navier-Stokes equation in Fourier space (see, e.g., Orszag

), or directly by Fourier transformation of the evolution equation for the velocity correlation tensor

R i j (r, t) ∂ ∂t + νk E(k, t) = S(k, t) + F(k, t). ( . )
Here S(k, t) arises from the convection term in Equation ( . ), F(k, t) corresponds to the external forcing, and contributions from the pressure term vanish identically. It has been mentioned before, that the nonlinear term does neither produce, nor dissipate energy. When Equation ( . ) is integrated, it should give the Fourier transform of Equation ( . ). is implies that the integral of S(k, t) over k has to vanish ∫ ∞ S(k, t) dk = .

( . )

S(k, t) can be interpreted as the transfer of energy into mode k.

We shall now follow Lesieur ( ) in a simple but revealing reasoning, that will nally put Richardson's interpretation on a theoretical basis. To this end we consider a stationary ow, i.e., we drop all the t-dependencies in Equation ( . ) and we de ne the energy ux

Π(k, t) = ∫ ∞ k S(k ′ ) dk ′ = -∫ k S(k ′ ) dk ′ ,
where the second equality follows simply from Equation ( . ). e energy ux is the amount of energy owing through the Fourier mode k. If it is positive, Π(k) > , energy ows towards larger wavenumbers and vice versa.

Now consider the following situation corresponding to Richardson's picture (comp. Figure . ). e external forcing is pumping energy into the system at a small wavelength k i . As we are in a stationary situation, the energy injection rate necessarily equals the energy dissipation rate and we have

ε = ∫ ∞ F(k) dk.
en the balance equation for Π(k), obtained by integrating Equation ( . )

Π(k) = ∫ k F(k ′ ) dk ′ -ν ∫ k k ′ E(k ′ ) dk ′ ,
may be investigated for two distinct regions. In the range k < k i , the integral over the external forcing is zero and in the limiting case of in nite Reynolds number, i.e., ν → , we obtain

lim ν→ Π(k) = -lim ν→ ν ∫ k<k i k ′ E(k ′ ) dk ′ = .
For k > k i the integral over the forcing becomes equal to ε and for some nite

k i < k ≪ ∞, one has lim ν→ Π(k) = ε -lim ν→ ν ∫ k i <k≪∞ k ′ E(k ′ ) dk ′ = ε.
If however we admit k → ∞, the integral in above equation must vanish in order to secure the validity of Equation ( . ), even when the viscosity ν tends towards zero, i.e., lim

ν→ ν ∫ k→∞ k ′ E(k ′ ) dk ′ = ε.
Summing up, we have for the energy ux in the three discussed regions in the limit of vanishing viscosity

lim ν→ Π(k) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ , k < k i ε, k i < k ≪ ∞ , k → ∞
erefore the energy, generated at large length scales (i.e., small wavenumbers), does ow towards the small scales, just as Richardson's intuitive picture predicted.

. . olmogorov scales "Moreover, I soon understood that there was little hope of developing a pure, closed theory, and because of the absence of such a theory the investigation must be based on hypotheses obtained in processing experimental data. "

A. N. Kolmogorov (Tikhomirov , p. ) In his seminal paper, Kolmogorov ( ) introduced two simple and powerful hypotheses on the statistics of locally homogenous and isotropic turbulent ows, which represented a crucial step in our understanding of turbulent ows. Kolmogorov ( ) starts by rigorously de ning the terms "locally homogenous" and "locally isotropic". e key assumption is that at scales much smaller than the forcing scale, the ow properties become homogeneous and isotropic. We present here a simpli ed version of the discussion, and refer to Frisch ( ) for a more thorough presentation. e main idea is conveyed by Kolmogorov's (

) " rst hypothesis of similarity", which states for distribution functions F n of small scale velocity increments or the locally isotropic turbulence the distributions F n are uniquely determined by the quantities ν and ε.

is statement allows us to nd the small length scale η, where energy gets dissipated by viscosity. It is uniquely de ned by viscosity ν and energy dissipation rate є and therefore must be η = ν є .

Likewise, corresponding time and velocity scales can be obtained

τ K = (ν ε) , u K = η τ K = (εν) .
ese quantities correctly describe the smallest eddies of the ow and are known as the Kolmogorov length, Kolmogorov time, and Kolmogorov velocity respectively.

. . olmogorov spectrum

Kolmogorov's ( ) "second hypothesis of similarity" extends his rst hypothesis and can be rephrased or intermediate scales, larger than the dissipation scale η, the distribution laws F n are uniquely determined by the quantity ε and do not depend on ν.

is statement may be used to derive a generic form of the energy spectrum E(k). Like for the Kolmogorov scales above, one employs simple dimensional reasoning. When in the inertial range the spectrum E(k) does only depend on the wavenumber k and, according to Kolmogorov's hypothesis, on the energy dissipation rate ε, it can only take the form

E(k) = C K ε k -. ( . )
Here C K is the dimensionless Kolmogorov constant. Sreenivasan ( ) reviews a large amount of data from experiments as well as numerical simulations and concludes, that the Kolmogorov constant is universally C K ≈ . . Interestingly Equation ( . ), was not given by Kolmogorov himself, but has been derived independently by Obukhov ( ), Onsager ( , ), and in collaboration by Heisenberg (

) and von Weizsäcker ( ).

----

E(k)ε - η kη Re λ = .
Re λ = .

F . Data from an experiment by Comte-Bellot & Corrsin (

). e straight dotted line gives the expected Kolmogorov scaling ∼ k -. Typically for lower Reynolds numbers, this scaling is only valid in a limited range. at range extends for larger Reynolds numbers.

.

--

In his paper entitled "On the E ect of the Internal Friction of Fluids on the Motion of Pendulums" Stokes ( ) derives the drag force acting on a small sphere of radius a, which moves uniformly with low velocity V through a uid at rest. He obtains (Eq. ( ) in the paper of )

F drag = πaµV , ( .
)
where µ = νρ f is the dynamic uid viscosity. From this result we can obtain in a hand-waving manner the equation of motion of a (small) particle moving freely in a turbulent ow. We assume that the size of the particle is much smaller than the smallest length scale of the ow, i.e., a ≪ η. It therefore sees a laminar ow around itself and we may apply Equation ( . ). e relevant velocity is the di erence between the particle velocity V and the uid velocity u(X, t) at the particle position X. In addition, we have to account for the (reduced) gravitational force acting on the particle ⁄ πa (ρ pρ f )G.

Assuming ρ f ≪ ρ p and therefore neglecting terms proportional to ρ f ρ p , we obtain

d dt X = V , ( . ) d dt V = u(X, t) -V τ p + G, ( . )
where we introduced the time scale

τ p = πa ρ p πaνρ f = ρ p ρ f a ν . ( . )
e time τ p is obtained by comparing the mass of the particle, ⁄ πa ρ p and the friction force, as expressed in the rst equality in Equation ( . ).

is very much simpli ed approach misses several e ects, such as the "added mass" or the "Basset history force" (Basset ; Boussinesq ; Oseen ). When the particle becomes bigger, further e ects described by the Faxén ( ) corrections come into play. Maxey & Riley ( ) and Gatignol (

) independently derived the full set of equations for freely moving particles. Assuming ρ f ≪ ρ p and a ≪ η in their results, all additional terms vanish, leading again to Equations ( . ) and ( . ). ey are therefore o en referred to as the Maxey-Riley equations. Also, throughout the work at hand, we only consider small (a ≪ η) but heavy (ρ f ≪ ρ p ) particles. Consequently the above set of equations in combination with the knowledge of the uid velocity u(x, t) at arbitrary points in space x, fully determines the trajectories of particles in the studied cases. e work of Daitche & Tél ( ) is a reminder that the history term, neglected in Equations ( . ) and ( . ), becomes important when the particle radius becomes comparable to or larger than η. For further information and an overview of recent advances, we refer the reader to the review article by Toschi & Bodenschatz ( ).

. . he tokes number e particle relaxation time τ p , Equation ( . ), may be compared to the shortest (fastest) time scale of the ow τ η = (ν є) . eir ratio de nes the dimensionless Stokes number

St = τ p τ η = ρ p ρ f a η , ( . )
which, in a sense, measures the particles' inertia. For St → the particles behave like tracer particles and follow the streamlines of the ow exactly.

We underline again, that our de nition of τ p , Equation ( . ), is only valid for heavy particles with ρ p ≫ ρ f . For light particles its de nition changes and therefore also the right-hand side of Equation ( . ).

Another handy dimensionless quantity is the particle Reynolds number

Re p = a u(X) -V ν.

It tells whether the assumption of smooth ow on the particle scale is valid-Re p ≪ in that case. erefore it is good practice to observe Re p , when one is using Equations ( . ) and ( . ), e.g., in a numerical simulation.

. . he particle velocity gradient tensor

For reasons that will become clear only later, we also introduce the particle velocity gradient tensor σ i j = ∂ j V i . Taking the partial derivative on both sides of Equation ( . ) and rearranging, yields

d dt σ i j = A i j (X, t) -σ i j τ p -σ ik σ k j ( . )
where the velocity gradient tensor of the uid A i j (x, t) = ∂ j u i (x, t) was introduced. is equation has the mathematical form of a Matrix Riccati equation. e solutions of this type of di erential equation typically display singularities. is type of behavior may at rst seem surprising, but it has physical signi cance, that will be discussed in Section . . . Furthermore, in Appendix B we discuss two ways to numerically overcome the caveat of those singularities.

INERTIAL PARTICLE COLLISIONS IN TURBULENT FLOWS

In the introduction (Section . ), we presented several problems, whose description relies on a good understanding of collisions between particles in a turbulent environment. Here, we discuss one of them, rain formation in clouds, in more detail. Another domain, the formation of planets in protoplanetary disks, will be rapidly reviewed. It will serve as justi cation for our choice of a broader parameter range than is of interest in the formerly mentioned context. Most technical aspects discussed in the coming sections are similarly applicable in other domains as well.

. e enormous range of involved length scales is commonly described as one of the most fascinating aspects of cloud microphysics, but it is also one of the most di cult to understand in quantitative terms. When approaching Earth from space, clouds are one of the rst distinct features one can perceive. ey can expand over lengths of hundreds of kilometers, while individual droplets and aerosol particles have radii of only several tenths of micrometers (Bodenschatz et al.

). Depending on what effects are taken into consideration, this range may even be extended in both directions (Siebesma et al. ). Evidently, clouds have a huge impact on our climate system and a better understanding is not only necessary for improved short term weather predictions, but also for long term projections of climate change-notably when aspects like "climate engineering" are discussed (Rickels et al.

). Consequently the interest in cloud physics has intensi ed in recent years. Particularly the in uence of turbulence has been widely discussed and related progress is documented in many review articles.¹ Especially physicists started to adopt the topic, as can be seen for example in the focus issue, that Falkovich & Malinowski (

) edited for the New Journal of Physics. Further reasons for the increased interest in cloud turbulence are certainly the improved numerical and experimental possibilities (Devenish et al.

). Here we shall present only the basic notions necessary for the understanding of the discussion of our work in the coming chapters. For individual details, we refer the reader to the already mentioned review articles. ). en comes a stage, where the droplets, continue to grow from initially ∼ µm to a few tens of microns. is step is in the beginning still dominated by condensation, but the range between µm and µm, known as the "size gap", is not very well understood (Grabowski & Wang

). Larger droplets can continue to grow by gravitational collisions in which they capture smaller ones with whom they coalesce. Di erent mechanisms have been proposed, that could drive droplet growth through the bottleneck of the size gap. One of them is turbulence, which can enhance the collision rate between medium sized droplets.

e whole process described above is conveyed for the ensemble of drops in a cloud by the droplet size distribution f (a, x, t). It gives the number density of droplets of radius a at position x and time t, in the sense that ∫ ∞ f (a ′ ) da = n, where n is the number density of all droplets in the volume

(Shaw

). We will focus on warm clouds, that is clouds, whose top is still below the freezing level and assume a spatially homogenous distribution. In this case the droplet size distribution obeys the Boltzmann-type equation

∂ f (a) ∂t = J(a) - ∂ ∂a [ ȧ f (a)] + ∫ a a a ′′ Γ(a ′′ , a ′ ) f (a ′′ ) f (a ′ ) da ′ -∫ ∞ Γ(a, a ′ ) f (a) f (a ′ ) da ′ , ( . )
where we introduced the substitution a ′′ = (aa ′ ) . e terms on the right-hand side consist of a particle source in form of the activation rate J(a), a term representing growth by condensation and the two integral terms, which convey production and destruction of particles of radius a by collisions. We refer the interested reader to Shaw ( ) for a more detailed discussion.

In terms of the droplet size distribution, the aforementioned "size gap" problem corresponds to the question, how a distribution, that is initially sharply peaked around some smallish value of a, can develop into a broader distribution displaying a nite probability for larger droplet sizes. e gures of Berry & Reinhardt ( a,b) illustrate this evolution nicely.

. . he collision kernel e integrands in the two integrals on the right-hand side of Equation ( . ) depend on the collision kernel Γ(a, a ′ ). It is de ned as the proportionality between the collision rate per unit volume of particles of radius a and a ′ and their respective number densities per unit volume c = nn ′ Γ(a, a ′ ).

( . )

In the monodisperse case, i.e., a = a ′ , a factor of ⁄ has to be added to the right-hand side to avoid double counting.

In realistic situations various phenomena contribute to Γ(a, a ′ ). For example, larger particles falling faster in the gravitational eld collect smaller ones, thus enhancing their mutual collision kernel. On the other hand, not all particles that come into contact will necessary coalesce-a fact typically accounted for by introducing a coagulation e ciency (see Section . . ). In the remainder of this work, we will entirely concentrate on geometric collisions due to particle-turbulence interactions. All other e ects, although important in more sophisticated models, will be neglected. is is a well established approach² allowing us to concentrate on a limited number of e ects. We achieve a further simpli cation by limiting ourselves to the case of monodisperse solutions, thus writing simply Γ(a) in that case.

Note that describing the droplet size distribution with Equation ( . ) requires the knowledge of the collision kernel Γ(a, a ′ ) for di erent particle radii a and a ′ . Technically, this can be determined from numerical simulations of bidisperse solutions. Such simulations have been carried out, for example, by Bec et al. (

). Typically the collision kernel is larger for solutions with particles of di erent sizes than for monodisperse ones.

e review of various models for the collision kernel Γ(a, a ′ ) presented by Ayala et al. (

) reveals some of the di erences between monodisperse and bidisperse solutions.

. . ypical values

Having brie y introduced the basic principles of droplet growth in clouds, we now discuss typical values of the two dimensionless quantities introduced in the last chapter-the Reynolds number and the Stokes number. Both depend implicitly on the energy dissipation rate ε found in clouds and the kinematic viscosity ν of air. e latter is tabulated and varies with temperature (Haynes ). We assume an intermediate value of ν = . × -m s. Measuring the energy dissipation rate in clouds is a complicated task. Usually, rst the energy spectrum is determined from velocity measurements during ights by balloon (Kitchen & Caughey

), aircra (Smith & Jonas ), or more recently by helicopter (Siebert et al. b)³, through clouds. en, from the spectrum, using Equation ( . ) with some empirical Kolmogorov constant, the energy dissipation rate ε is determined.

e measured value di ers

= ( ν ε) (Pope , p.
), the energy dissipation rate ε and the root mean square velocity u rms su ce to determine the Taylor scale Reynolds number Re λ -see Equation ( . ). It is, again depending on the same various factors as ε, of the order of to (Siebert et al.

; Siebert et al. a). To estimate the Stokes number of droplets, we further need the densities of water and air, as well as the size of the droplets. e range of the latter has already been mentioned. e density of water at typical temperature and pressure is roughly ρ W ∼ kg m , the density of air can be estimated ρ A ∼ kg m (Haynes ). With these values we nd from Equation ( . ) an upper bound for the Stokes number of about St < .

Another important dimensionless number is the volume fraction Φ. Consider a given volume V containing N p particles of size a, then the volume fraction is Φ = N p πa V . Its exact value depends on the droplet size and the type of cloud, but is typically of order -(Grabowski & Wang

). For a more detailed discussion, we refer to Pruppacher & Klett ( , Sec. . . ).

A word of caution e simplifying approach o en employed by physicists has repeatedly been criticized by the atmospheric research community (e.g., Grabowski & Vaillancourt ; Khain et al.

). is criticism is certainly justi ed, when one attempts to transpose the results to the cloud physical context without any further discussion. We argue however, that these simpli cations are necessary to isolate di erent phenomena and to gain a deeper insight. In a second step, the improved understanding can be adopted in more complete models. erefore, when in the coming chapters we present our results, we will deliberately extent the investigated range of Stokes numbers and neglect, as has been mentioned, aspects, like gravitation, that are incontrovertibly of huge importance in the context of clouds.

Dust grains in protoplanetary disks

Further motivation for our choice of an extended parameter range is provided by the fact that larger Stokes numbers may be important in other domains. In the astrophysical context where turbulence induced particle collisions are discussed as a possibly important factor in the formation of planets (Johansen et al.

; Safranov ; Shari ), the experimental data is much more uncertain. But Stokes numbers of O( ) and larger are expected (Pan et al.

; Wilkinson et al.

). For a short overview of the domain, the interested reader is referred to Beckwith et al. (

). We . e collision kernel Γ has been de ned in Section . . . In this section, we present di erent limiting cases, where the collision kernel can be analytically derived, as well as other important results.

. . ollision kernel for an ideal gas

To get a feeling for the typical reasoning, we start by deriving the collision rate in an ideal gas. To this end we follow the textbook by Moore (

). e result obtained in this section will be used later, to verify our numerical collision detection scheme.

Imagine particles with a number density n, each moving with a constant normally distributed velocity V . e collision rate for any such particle can be calculated by considering the collision cylinder as shown in Figure . . If all other particles were at rest, the considered particle would collide with all particles in a cylinder with radius a around its trajectory. erefore, the corresponding collision rate per unit time would be n⟨ V ⟩π( a) . e other particles do however move and hence the correct velocity to determine the cylinder's height is the mean relative velocity ⟨w⟩ = ⟨ V ⟩. To obtain the total collision rate per unit volume, the above quantity needs to be multiplied by n , nally leading to

c = n √ π( a) ⟨ V ⟩. ( . )
Based on this result and according to Equation ( . ) we can then de ne the collision kernel for an ideal gas Γ ig = π( a) ⟨ V ⟩.

. . ero inertia

Owing to its importance, but also to the di culties involved in this task, much e ort has been made to determine the collision kernel for turbulent ows. e paper by Sa man & Turner ( ) introduces some of the seminal ideas. It describes a way to theoretically derive an estimate for the collision kernel of two particles of identical size, which follow the ow exactly, i.e., for the case St → . e main idea is illustrated in Figure . . For a sphere of two times the particle radius a around one central particle, the total rate of in ow is determined. Multiplied by n, the number density of particles in the ow, this quantity determines the rate of collision for the particle at the center. Mathematically this idea is described with the help of an integral over the surface of the sphere

Ω = n ∫ -w r ( a, Ω)Θ[-w r ( a, Ω)] dΩ. ( . )
is integral sums the radial component of the uid velocity w r (x, t) on the sphere, whenever it is negative-the Heaviside step function Θ(⋅) assures that condition.

For an incompressible uid the in ow and the out ow into a given volume cancel exactly. erefore the integral in Equation ( . ) simpli es to ∫ -w r ( a, Ω, t)Θ[-w r ( a, Ω)] dΩ = ∫ w r ( a, Ω, t) dΩ. ( . ) Up to this point, only one central particle has been considered. To obtain an estimate for the overall collision rate, it is necessary to average . is further simpli es the situation, because for isotropic statistics it is now su cient to consider ⟨ w x ( a) ⟩. For small particles with radius a ≪ η this can be identi ed with ⟨ w x ( a) ⟩ = a⟨ ∂u x ∂x ⟩. Taylor ( ) had found ⟨(∂u x ∂x) ⟩ = ε ν. With this result, assuming Gaussian statistics for the velocity gradients, Sa man & Turner (

) nally obtain

c = n Γ ST , with Γ ST = ( a) τ K π . ( . )
In their paper Sa man & Turner go on to present a derivation of the collision rate for inertial particles-see Section . . for a short comment on that part.

. . n nite inertia

Abrahamson ( ) noted that most previous studies aiming at modeling the collision kernel for particles in a turbulent ow dealt with the case of particles following the ow (St → ). He argued that this condition is however not ful lled in many situations and therefore derived the collision rate for a "vigorously turbulent uid", which corresponds to the case St → ∞. In such a situation, particle velocities, even for those passing close to each other, can be considered randomly distributed. e correct approach to derive the corresponding collision rate is the same as for the ideal gas described in Section . . . Hence the problem reduces to estimating the average relative velocity between the particles. Assuming the velocities are Gaussian distributed, the average velocity can be related to the root mean square of its components

⟨ V ⟩ = π V rms . ( . )
Based on Corrsin's hypothesis on the relation between Lagrangian and Eulerian correlation functions and following former work by others, Abrahamson ( ) obtains from a simpli ed version of the Maxey-Riley equations (Section . ), adequate for particles with large Stokes numbers,

V rms = u rms + τ p ε u rms = η τ K Re λ √ + St Re λ ,
where in the second step ε = η τ K and (u K u rms ) = Re λ (see, e.g.,

Pope

) have been used. Given the empirical reasoning involved in the derivation of the above result, the exact factors may vary, but it is important to note the functional form V rms = (η τ K )γ(St, Re λ ). Combining this result with Equations ( . ) and ( . ), Abrahamson ( ) obtains the collision kernel

Γ A = √ π( a) η τ K γ(St, Re λ ), ( . )
which di ers from Sa man & Turner's result by a factor ∼ η a.

Mehlig et al. (

) point out, that Abrahamson did not account for the "multiscale nature of the ow". ey provide a di erent approach based on Kolmogorov-type reasoning: In the inertial range, V rms can only depend on the dissipation rate ε and the particle relaxation time τ p . By dimensional analysis they nd γ(St, Re λ ) ∼ St with a possibly Reynolds number dependent prefactor.

. . referential concentration

In both turbulent cases presented above-zero and in nite inertia-the particles are distributed homogeneously. Maxey ( ) however noted, that for inertial particles with intermediate Stokes numbers, the particle velocity eld becomes compressible and particles tend to cluster in certain regions of the ow. is becomes obvious, when one writes the formal solution of Equation ( .)

V (t) = ∫ t e (t ′ -t) τp u(t ′ ) τ p + G dt ′ .
Assuming small inertia, the above equation can be developed in τ p by consecutive integration by parts

V = u + τ p G -τ p u + O(τ p ), ( . )
where u has to be interpreted as the material derivative ∂ t u + u ⋅ ∇u. Taking the divergence of Equation ( . ) yields

∇ ⋅ V = -τ p ∇ ⋅ (u ⋅ ∇u) = -τ p (∂ i u j )(∂ j u i ) = τ p [(∂ i u j + ∂ j u i ) -(∂ i u j -∂ j u i ) ].
e rst term in the square brackets can be related to the so called "rate of strain tensor", the second term to the vorticity (see Section . . ). Maxey ( ) concluded, that the particle velocity eld is divergent in regions of high vorticity, while it is convergent in regions of high strain. erefore the particles would get ejected from vortices and cluster in regions of low vorticity.

Later work by Wilkinson et al. ( ) and Gustavsson & Mehlig ( b) provides an alternative approach. eir analytical results, which they obtain for a model ow, are able to explain features present in numerical simulations of turbulent ows.

e clustering, or preferential concentration, evidently needs to be taken into account, when one determines the collision kernel. In the following, we present the approach introduced by Sundaram & Collins ( ), which is based on ideas from statistical mechanics.

Consider a system with N p particles. We de ne the probability that particle " " is in dX at X , particle " " is in dX at X , etc. P (Np) (X , . . . X Np ). en the probability that any two particles are in dX at X and in dX at X is (e.g., McQuarrie )

ϖ ( ) (X , X ) = N p (N p -)P ( ) (X , X ). ( . )
A further useful de nition is the radial distribution function (RDF)

g ( ) (X , X ) = n -ϖ ( ) (X , X ).
For isotropic statistics all above quantities only depend on the two particles' distance r = X -X . In this case one has P ( ) (X , X ) = P(r) V and can de ne

g(r) = N p (N p -) n V P(r), ( . )
where we dropped the upper index for simplicity.

Here we omit the argument X(t), which should in principle appear in V and u. See the original paper by Maxey ( ) for a more rigorous treatment.

To decide whether two particles are about to collide, one needs, aside from their distance, knowledge about the relative velocity w = V -V .

erefore, instead of the above probability P ( ) (X , X ), one has to consider the joint probability of position and velocity P ( ) (X , V ; X , V ). Here again, we can make use of isotropy, which means P(r, w) is enough to describe the statistics. e latter can be benecially rewritten with the help of the conditional probability P(w r)

P(r, w) = P(w r)P(r) = P(w r)g(r) n V N p (N p -) , ( . )
where in the second equality the RDF from Equation ( . ) was introduced.

It is convenient to introduce a function that, based on the values of r and w, recognizes whether a pair of particles is going to collide. Sundaram & Collins (

) rigorously construct a function ψ(r, w; τ), which returns if the particles collide within a time τ, and if not. We shall see, that the derivative with respect to τ is enough to determine the collision rate. Sundaram & Collins ( ) obtain

dψ(r, w; ) dτ = -w r Θ(-w r )δ( a -r), ( . )
where

w r = (V -V ) ⋅ (X -X )
X -X is the radial relative velocity (RRV) and Θ(⋅) and δ(⋅) are respectively the Heaviside and Dirac functions. Note that this is again simply Sa man & Turner's ( ) integral kernel-see Equation ( . ). Now we put together Equations ( . ) and ( . ), as well as the collision operator, to determine the number of collisions that happen in a time τ N c (τ) = N p (N p -) ∫∫ ψ(r, w; τ)P(w r)g(r) n V N p (N p -) dr dw.

From N c (τ) the collision rate per unit volume may be obtained by derivation with respect to τ and division by

V c = V lim τ→ N c (τ) τ = V dN c (τ) dτ .
Noting that ψ(r, w; τ) is the only τ-dependent term and introducing Equation ( . ), Sundaram & Collins ( ) nally obtain c = n π( a) g( a) ∫ -∞ -w r P(w r a) dw r .

Making the same argument as in Equation ( . ), one can write for the colli- which di ers from Sa man & Turner's collision kernel, Equation ( . ), only by the appearance of the RDF, that accounts for preferential concentration. But note also, that-as opposed to the former results-no explicit expression is given for the RRV ⟨ w r ⟩. ), and will be discussed shortly towards the end of this section. ) provide a di erent explanation for the observed enhancement of the collision rate.

sion kernel Γ SC = π( a) ⟨ w r ⟩g( a), ( . )
An initially single valued distribution of inertial particles in phase space may become multivalued due to the formation of so-called "caustics". is e ect is similar to the formation of shocks in Burgers' equation (Frisch & Bec ) or to the formation of the characteristic light pattern that can be observed on the ground of a swimming pool, from which this e ect borrows its name ( Berry). Consider the one-dimensional example depicted in

X V t = t t > t F .
e particle phase space is spanned by their position X and velocity V . e lines may be interpreted as manifolds on which the particles are distributed homogeneously. Although initially there is only one particle at each position, particles overtake each other leading to the formation of a fold caustic (comp. Wilkinson et al.

).

Figure . . Over the course of time, faster particles will overtake slower ones. ose with negative velocities will move to the le , those with positive velocities will move to the right. us a region, where the particle velocity eld is multivalued develops in coordinate space. is has two consequences, which can potentially lead to collisions. On the one hand, particles may nd themselves close to each other, but with largely di ering velocities. In this case Abrahamson's ( ) approach for uncorrelated particle motion is more appropriate to describe the collision rate, than Sa man & Turner's (

). On the other hand, the number density of particles diverges at the edges of the caustics. Wilkinson & Mehlig ( ) had investigated this in an earlier publication, but do consider its contribution to the collision rate negligible in systems, where nη ≪ .

In the presence of caustics, an obvious decomposition of the collision

kernel is Γ = Γ ST g( a) + Γ A h S (St, Re λ ), ( . )
where the RDF g( a) has been introduced according to ( . ) to account for the e ect of preferential concentration, which a ects only shear induced collisions. e function h S (St, Re λ ) can be interpreted as the fraction of position space, where the velocity has become multivalued due to the formation of caustics (Wilkinson et al.

). is quantity can be related to the rate of caustic formation, which Wilkinson & Mehlig ( ) had determined analytically for a model ow in a previous publication. From their results it is nally concluded, that h S (St, Re λ ) = exp(-C St) with some universal dimensionless constant C. Duncan et al. ( ) determined this constant for particles advected by a three-dimensional random vector eld. In a subsequent work Falkovich & Pumir ( ) could con rm the functional form of h S (St, Re λ ) in a DNS of turbulent ow. ey found in addition, that the constant C falls with growing Re λ .

Falkovich et al. (

) predicted the sharp increase in the collision rate using a di erent string of arguments. ey imagine particles getting ejected out of vortices, like stones from a sling. ese particles will have a large RRV when they collide, as well as a totally uncorrelated motion with the particles they collide with. e caustics and sling e ects turn out to be essentially the same, despite the two di erent approaches and presentations.

Consider Figure .  again. It is obvious that the formation of a caustic is accompanied by the divergence of the quantity dV dX on its edges. In three dimensions this translates to singularities of the particle velocity gradient tensor σ i j as introduced in Equation ( . ). erefore this quantity was used by Falkovich & Pumir ( , ) to study the sling e ect. If one has access to σ i j for each particle in every collision, it is even possible to estimate the ratio of sling collisions. ). ey note that the velocity eld describing the motion of inertial particles can be decomposed into two contributions. One describing the smooth spatially correlated movement of the particles and another one describing "random uncorrelated motion" (Reeks et al.

). e latter is responsible for caustics and sling collisions.

IJzermans et al. (

) present a method to study similar e ects based on the deformation of an in nitesimal volume around a particle along its trajectory.

is work was extended by Meneguz & Reeks (

). e relationship between the di erent interpretations described in this section is discussed and studied in two model ows by Gustavsson et al. (

).

. . host collision approximation

A typical simpli cation, when studying collisions in numerical simulations is the so-called "ghost collision approximation" (GCA). A collision is detected whenever the distance between two particles falls below the collision radius a, but the particles do not interact in any way and are allowed to overlap. ) su er from this inaccuracy. When in Equation ( . ) the amount of in owing uid is summed up by the integral, it is not asked, whether it is "fresh uid", which enters the sphere. In a turbulent ow it is however possible, that a parcel of uid traverses the collision sphere several times. In this case, only the rst passage should be counted in for the collision rate. Sa man & Turner considered a locally hyperbolic ow as shown in Figure . , which they assumed to be persistent. is is however incorrect: a turbulent ow is neither at every instant and space point hyperbolic (Bec ; Chong et al.

) nor is its structure persistent in time (Brunk et al.

b). To circumvent the problems associated with uid elements that pass the collision sphere multiple times, Andersson et al.

(

) propose to incorporate an indicator function in Equation ( . ), such that

= n ∫ -w r ( a, Ω)Θ[-w r ( a, Ω)]χ( a, Ω) dΩ. ( . )
e indicator function χ is unity only if the point reaching the surface of the sphere at (Ω, a) has not previously traversed it, elsewise it is zero.

is function is analytically hard to determine though, as in principle it necessitates knowledge about the full history of every uid element.

We nally mention the work by Wang et al. ( b) who study the combined e ect of the GCA and hydrodynamic interactions in terms of a collision e ciency (see next section).

eir work is interesting from an engineer's point of view, as it provides a simple way to compensate the bias introduced by the GCA. From a physical point of view it is however unsatisfactory as it seems to miss a conclusive explanation of the origins of that bias. One of the main objectives of this work has been to provide a better understanding of the biases introduced by the GCA. Our results are presented in Chapter .

. . ollision e ciency

In most of the work dealing with collisions in turbulent ows, the interactions between particles are e ectively neglected. In particular, hydrodynamic forces are not taken into account. ese interactions may however cause a reduction of the collision rate (Brunk et al. a,b, demonstrate this for coagulating particles). Usually these e ects are described with the help of an e ective collision radius, from which the collision e ciency ≤ c ≤ can be obtained (Pruppacher & Klett ). is e ective collision e ciency is then introduced in the above expressions for the collision kernel as a multiplicative factor. We refer the interested reader to Pruppacher & Klett (

) for an overview.

. . wo interesting remarks

In the late s and early s some interesting remarks concerning the derivation of collision kernels were made by L.-P. Wang and co-workers.

Here we brie y summarize two important results that will be helpful in the coming chapters.

Wang et al. ( b) noted that two di erent ways to model collisions were used in theoretical descriptions of the collision kernel. On the one hand, there is the "cylindrical formulation", similar to the one we presented in Section . . , on the other hand there is the "spherical formulation" as used by Sa man & Turner (Section . . ). e former makes the assumption of uniform relative velocity w, at least on the scale of the collision radius a. is assumption is correct in the case of Maxwellian particles, it fails however for turbulent ows at intermediate Stokes numbers. In addition to the one we presented in Section . . , Sa man & Turner ( ) give a second derivation of the collision kernel, this time for inertial particles. at result di ers from the one in Equation ( . ), but they explain the discrepancy with inexact assumptions made on the turbulent ow. ). It explains an e ect which may seem puzzling at rst: When one determines, for example in a numerical simulation, the average radial relative particle velocity ⟨ w r ⟩ c , based only on colliding particles, then the result di ers from the one obtained, when one takes into account all pairs at a distance a. Obviously, the correct average ⟨ w r ⟩ is the one that considers all pairs, but interestingly the two values can be related. To this end it is helpful to re ect the procedure, how ⟨ w r ⟩ c would be calculated. Consider again Sa man & Turner's (

) picture of one central particle and especially only one surface element dΩ of the collision sphere. For each such surface element, the radial relative velocity w r is registered and summed up at a rate d = -nw r ( a, Ω)Θ[-w r ( a, Ω)] dΩ (comp. Equation ( . )). By integrating over the sphere, averaging, and multiplying ⁄ n to account for all particles, one obtains the sum of w r from all collisions per volume and time, ( )

- collisions w r ≈ n π( a) ⟨ w r ⟩.
In this casual notation and simply mean the observation volume and time. From the above quantity the average radial relative velocity conditioned on collisions, ⟨w r ⟩ c , is obtained by normalizing with the number of collisions per volume and time c -Equation ( . ). One nally obtains

⟨ w r ⟩ c = ⟨w r ⟩ ⟨ w r ⟩ , ( . )
which is the result given by Wang et al. ( ). Although maybe confusing at rst, this result is numerically veri ed and will be of importance in Section . . Furthermore Ducasse (

) discusses interesting implications that arise from this result for the PDFs of collision angles.

Part II MAIN SCIENTIFIC WORK OF THE THESIS COLLISION DETECTION AND IDEAL GAS PARTICLES

e work carried out in this thesis is based on numerical simulations of turbulent ows, in which we follow particles. e detection of collisions between particles requires some careful treatment of the numerical data, which we explain in this chapter. Some more technical aspects will be discussed in the appendix.

In Section . . the collision kernel for an ideal gas was derived. is simple case will be used in the following to check our scheme. In addition some results, which are of interest in the coming chapters, will be derived for this simple example.

. To detect collisions in a system of N p particles, one may want to check all particle pairs, which requires of the order of O(N p ) operations. In some of our simulations we have up to O( ) particles. With such a number of particles, the simple-minded approach leads to a prohibitively expensive computational cost. It is however common to use a technique originally developed in molecular dynamics simulations called "cell linked-list" algorithm (Allen & Tildesley

) to reduce that cost. Sundaram & Collins (

) describe an implementation of this technique to study collisions of inertial particles in a turbulent ow.

. . ell linked-list algorithm e basic idea of the cell list algorithm is, that only particles in a certain region around one central particle can collide with that particle during the next time-step. To make use of this observation, the (cubic) computational domain is divided into N d b boxes, where d is the dimensionality of the system-in our case d = . In a rst step, a list of particles in each box is generated. en, each particle in every box is checked for a collision with any particle in the ( d -) surrounding boxes. In the ideal case, N b would be chosen such that each box contains only one particle (or even less). en the computational cost would reduce to a mere O(N p ). ere is however a minimal size for the boxes, which is given by the maximal distance, one particle can travel within one time-step. When the simulation domain has length L sys in all directions, each box has length L b = L sys N b . Now consider the situation illustrated in Figure . -two particles move with maximal velocity V max perpendicular to the cell. One particle is located at the far end of the cell under consideration, the other one on the near end of the next but one cell. If the collision between these two particles would take

V max ∆t V max ∆t L sys N b F .
Illustration of a limiting case, where a collision may take place, but would not be detected because the box size of the cell list algorithm is too large.

place within the next time-step, it would be missed by the cell list algorithm.

erefore the maximal number of boxes per direction is

N b = L sys [ (V max ∆t + a)],
where ∆t is the time-step and a is the particle radius. Typically this constraint is not too restrictive. But other e ects, like particle clustering can reduce the e ciency of this approach. ere are other algorithms, like k-d trees, which require of the order of N p log N p operations to determine nearest neighbors, and can therefore further reduce the computational cost of detecting collisions. But in systems, where the particles are su ciently homogeneously distributed and where the constraint on N b is not too strong, they are unlikely to perform better than the cell linked-list algorithm. is showed also in some comparative tests we did with a freely available implementation of the k-d tree algorithm (Kennel ).

. . rajectory interpolation

In Section . . it has been described, how the number of particle pairs that need to be checked for eventual collisions can be reduced. It was however not explained, how a collision could be detected for a given pair. Simply measuring the distance between the two particles is not enough, because one cannot tell whether they have been in contact before. A more elaborate approach consists in taking into account the two particles' positions and velocities at one time-step. en one can extrapolate their trajectories linearly and determine whether they approach closer than a within the next time-step. With a su ciently small time-step, this approach is correct. But usually one has interest to keep the time-step as large as possible. In this case problems like the one illustrated in Figure . may arise. erefore we used a more sophisticated approach in our simulations. Particle positions and velocities of two consecutive time-steps are stored. is allows a third-order interpolation of the trajectories, which is more accurate. Comparing the two approaches in a typical situation, we found that the linear extrapolation detected up to % collisions in excess and missed a few others. In Appendix A both approaches are discussed in more detail.

t (n) t (n+ ) t (n+ )
distance F . e dashed line shows the exact distance between a pair of particles over the course of two time-steps.

e straight line corresponds to a linear extrapolation at time-step t (n) .

is falsely predicts a collision. At the next timestep t (n+ ) a second collision would be counted.

. . umerical determination of the main variables e variables we are most interested in are the collision kernel Γ, the mean RRV ⟨ w r ⟩, and the RDF g( a) at the collision distance a. All these quantities have been introduced in Chapter . In the following paragraphs we explain how they are determined numerically.

e collision detection scheme allows to determine the number of collisions N c (τ) that take place in a time interval of length τ.

en the collision rate per time and volume in this interval is

c = N c (τ) τ sys , ( .
)
where sys is the volume investigated for collisions. From Equation ( . ) the instantaneous collision kernel can be derived

Γ τ = sys N c (τ) τN p . ( .
)
e nal value for the collision kernel Γ is obtained by averaging over consecutive intervals. is procedure also provides an estimate for the uncertainty in the determination of the collision kernel, namely the standard error of the mean (Sachs & Hedderich ). For the length τ of the interval, we usually choose a large time scale of the ow. In the DNS, which will be introduced in the next chapter, we use the large eddy turnover time T L .

To determine the RDF at contact g( a), the cell-linked list algorithm gets employed again, but this time to detect all particle pairs that have a separation a -∆r < r < a + ∆r , where ∆r is some very small distance. Di erent values of ∆r were tested. If the value is chosen too small, one will detect almost no pairs and therefore introduce a large statistical uncertainty. If on the other hand ∆r is too large, the variation of the function g(r) can not be neglected anymore in the range a ± ∆r . All results reported on in this work were obtained with ∆r = . a, which was found to give correct results. Given the number N a of particle pairs with a separation a at one time-step, the momentary RDF can be obtained from Equation ( . )

g( a) ∆t = N a sys ⊚ N p , ( .
)
where ⊚ = π[( a + ∆r) -( a + ∆r) ] is the volume of the shell around each particle, in which other particles were searched. Note that Equation ( . ) gives the RDF at one instant in time, i.e., at one time-step. It is not necessary to apply the interpolation used to detect collisions in between time-steps. Again, the nal result g( a) is obtained by averaging over all time-steps.

As before, all particle pairs in a distance a of each other are detected. For each such pair the RRV w r is calculated and by averaging over all pairs and all time-steps, one obtains ⟨ w r ⟩.

. In this section we consider a system of particles, whose velocity components are distributed according to a Gaussian

p(V i ) = ς √ π exp - V i ς , ( . )
where the notation ς = V rms was introduced for clarity of presentation. is situation corresponds to the classic model for an ideal gas. In Section . . , the collision kernel for this case was derived. We chose this special system for two reasons: First, it is useful to calibrate our collision detection algorithms. Second, it will allow a relatively simple access to a rst theoretical description of some not yet studied aspects.

. . eri cation of collision detection scheme

We applied our implementation of the techniques described in the above sections, to a set of particles with constant (in time) velocities according to Equation ( . ). e expected mean velocity is ⟨ V ⟩ = ςπ -. erefore the collision kernel can be expressed in terms of ς as

Γ ig = √ πς( a) .
We have determined the collision kernel in our numerical simulations by counting the number of collisions per time and volume. Table . shows a comparison between the results and the theoretical predictions for three di erent ratios of the particle radius a over the size of the simulation box L sys . e boundary conditions were, as always throughout this work, periodic. e results con rm the correctness of our implementation.

.

a L sys Γ 4 √ πσ(2a) 2
2.1 × 10 -4 0.9995 ± 0.0003 3.0 × 10 -4 1.0001 ± 0.0007 3.8 × 10 -4 1.0003 ± 0.0014

T .

e theoretical prediction and our measured collision kernel match very well, con rming our collision detection scheme. e uncertainty increases due to slightly poorer statistics for larger particles.

. In Section . we will be interested in the time τ c that two inertial particles spend closer than a certain distance d c in a turbulent ow. For the choice d c = a this time corresponds to the contact time of particles. In this section we will derive the PDF of this quantity for an ideal gas of particles with radius a. To this end we start by deriving the PDF of relative velocities only for collisions.

. . robability of relative velocity and statistics conditioned on collisions e PDF of the relative velocity w = V -V between two particles in a Maxwellian gas can be determined from Equation ( . ). Elementary algebraic manipulations lead to

p(w) = ( √ πς) exp - w ⋅ w ς , ( . )
where it was assumed that all particles and all velocity components are independent.

To determine the PDF of w, conditioned on the fact that the particles are about to collide, one proceeds as described in Section . . . One obtains in analogy to Equation ( . )

p c (w c ) = ∫ δ(w -w c ) -w r Θ(-w r ) ⟨ w r ⟩ p(w) d w, ( . )
where the Heaviside step function Θ(⋅) picks only negative RRVs. e index c indicates that the PDF is restricted to collision events. Equation ( . ) can be evaluated for an ideal gas, if one takes into account isotropy. Without loss of generality one may assume that w r = w z . Inserting Equation ( . ) yields

p c (w) = w ⋅ ẑ ς π Θ(w ⋅ ẑ) exp - w ⋅ w ς , ( .
)
where ẑ designates the unit vector in z-direction. Note that we dropped the index c on the relative velocity w, assuming that it su ces to indicate it at the probability distribution itself.

. . robability of contact times

A er the validation of our collision detection algorithm, we discuss now the PDF of the contact time τ c as de ned in the beginning of this chapter.

We restrict ourselves to d c = a, but obviously more general results may be obtained by simple substitution. In the case of particles that are moving with a constant velocity without changing direction, the collision time can be easily related to their relative velocity (see Figure . )

τ c = a cos θ w = aw r w ,
where θ is the collision angle.

So the contact time is entirely determined by the relative velocity w at the moment of collision. erefore we may obtain its PDF with the help of the Dirac delta function from Equation ( . ). We start by placing ourselves in a spherical coordinate system

p(τ c ) = ∫ ∞ ∫ π ∫ π δ τ c - a cos(θ) w × w cos(θ) ς π exp - w ς w sin(θ) dw dθ dϕ. ( . )
e Heaviside function from Equation ( . ) has already been evaluated, such that the θ-integration spans only the interval [ , π ]. Note also that the index c at the probability itself has been omitted. e integration in θ can be further simpli ed

∫ π δ τ c - a cos(θ) w cos(θ) sin(θ) dθ = ∫ δ(τ c -a w ξ)ξ dξ = w a ∫ δ wτ c a -ξ ξ dξ = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w a τ c if w ≤ a τ c , else.
Introducing this in Equation ( . ) leads to

p(τ c ) = τ c ( ς) ( a) ∫ a τc w exp - w ς dw,
which can be solved by standard techniques Again we performed numerical analysis. e results shown in Figure . con rm our derivation. e PDF consists basically of an increasing branch proportional to τ c for small τ c to and a rapidly decreasing branch proportional to τ - c for large contact times. It peaks at about τ c = a ς, a value that is small compared to the typical time-step in a numerical simulation.

p(τ c ) = τ c ς a -+ a ςτ c + a ςτ c exp - a ςτ c

DIRECT NUMERICAL SIMULATIONS

In the last chapter we described, how collisions are detected given the trajectories of particles. To actually determine these trajectories, we integrated the Maxey-Riley equations ( . ) and ( . ), which in turn necessitates knowledge of the uid velocity at arbitrary points in space. In this and the coming chapters we present results for which the uid velocity eld was obtained by direct numerical simulation (DNS). Later, in Chapter , ndings from a di erent approach, using a simple model for the uid eld, will be shown.

Here we start with a summary of our DNS.

. We obtained the velocity eld u(x, t) with the help of a standard pseudospectral code. is numerical method makes use of the fast Fourier transform (FFT) to integrate the Navier-Stokes equation ( . ) in Fourier space. We note that one such code has been made publicly available under the GPL by Chumakov (

) and refer the reader to the introductory book by Boyd ( ) for further information. e code we use was originally developed by Emmanuel Lévêque and Christophe Koudella at the ENS Lyon. It has been used in several publications (e.g., Calzavarini et al.

) and was optimized to be especially e cient on the cluster of the PSMN. A periodic cube of grid points was used in our simulations with a dealiasing following Orszag's (

) two-thirds rule, therefore leaving e ectively modes. e forcing term was chosen to continuously excite a small band of low wavenumbers such that the energy injection rate ε remains constant (Lamorgese et al.

). e energy injection rate and the viscosity can be freely chosen and were in our simulations ε = -and ν = × -respectively. ese values are given in dimensionless (or rather DNS) units. As explained in Section . . one may scale length and time in order to compare our results to di erent systems with the same Reynolds number. e Taylor scale Reynolds number of our homogenous isotropic turbulent ow is Re λ =

. In ) and is typical for this type of simulation. We nally note that time-stepping was achieved with a second order Adams-Bashforth scheme bearing in mind the Courant-Friedrichs-Lewy condition Co = u rms k max ∆t ≲ . (Courant et al. ) for the time-step ∆t.

-----

E(k)ε - η kη F .
Typical energy spectrum in our DNS (line) compared to an experimental spectrum (Comte-Bellot & Corrsin , Re λ = . , symbols). e inset shows the same spectra compensated according to the Kolmogorov scaling, see Equation ( . ), along with the value C K = . (comp. Sreenivasan

).

T . Characteristic values of the DNS. Given are the energy injection rate ε, the kinematic viscosity ν, the root mean square velocity u rms , the Taylor scale Reynolds number Re λ , the Kolmogorov length scale η, the Kolmogorov time τ K , the longitudinal integral length L l , and the integral (or large eddy turnover) time T L .

ε ν u rms Re λ η τ K L l T L . × - . × - . . . . .
Table . summarizes the values that characterize our DNS. Among these values the longitudinal integral length L l is given, which is calculated in the DNS according to the Fourier space equivalent of Equation ( . )

L l = π u rms ∫ k -E(k) dk ( . ) (e.

g., Batchelor

). e large eddy turnover time T L = L u rms is derived from this value and the root mean square velocity u rms .

. . article trajectories

Once the DNS had reached a statistically steady state, N p particles were introduced. We performed several runs with particles of di erent Stokes numbers, but in each such run, all particles had the same Stokes number. More details about the runs will be given in the next section. e particle trajectories were integrated by the velocity Verlet algorithm (Press et al.

; Swope et al.

)-a second order scheme-and the uid velocity at each particle's exact position was obtained by tri-cubic interpolation. A er a transient time of the order of T L , all particle positions and velocities were stored at a rate of one per . τ K , although the actual time-step used in the integration was a factor smaller. is data was then post-processed for all subsequent studies. At rst this procedure may seem cumbersome, but it allowed us to adapt our treatment of the data without the additional cost of integrating the trajectories anew.

e above mentioned transient time is necessary, because particles are inserted homogeneously distributed in the uid and with velocities equal to the uid velocity. It has however been described before, that inertial

. Γτ K ( a) t T L τ = ∆t τ = τ K τ = T L
F . e collision kernel Γ during the initial time period.

e result is shown for di erent sampling rates τ, according to Equation ( . ). e initial transient takes at least T L . In our simulations, we wait for T L before recording Γ with a sampling rate of τ = T L . particles in a turbulent ow are not homogeneously distributed, but tend to cluster in certain regions. Also, their velocities are di erent from the uid velocity. erefore the freshly introduced particles need a time of the order of a few large eddy turnover times to " nd their places" in the ow. is transient time is illustrated by Figure . , which shows the evolution of the instantaneous collision rate during this period.

Due to the very long integration time and a minor neglect in the setup, very few particle trajectories (less than % of the whole population) became de cient a er a certain simulation time. As we could deal with this on the level of the collision detection scheme, we decided not to repeat the integration of the trajectories. In the collision detection code, we check, whether any particle runs the risk of becoming de cient. If that is the case, the corresponding particle is replaced, by one additional particle that had not been under consideration for collisions before. e disquali ed particle gets then ignored for the rest of the simulation.

is approach was easy to implement, because we usually track more (ca. %) particles than are actually considered for collisions. e reasons for this procedure will become clear in Section . . . . As mentioned above, several runs were performed at di erent Stokes numbers (Section . . ). e di erent runs are summarized in Table . . e number of particles was determined by our choice of the volume fraction Φ = N ′ p πa L sys , where L sys = π in DNS. e volume fraction was typically Φ = . × -, except for three cases, where it was Φ = . × -. Both values can be considered dilute and e ects involving three particles and more can be a priori neglected. In Section . . we will discuss a case though, where such e ects have a measurable in uence. Instead of simply integrating the N ′ p particle trajectories necessary to achieve those volume fractions, we integrated an additional % of particles amounting to a total of N p = ⁄ N ′ p particles. e reason for this will be explained in Section . . . Usually it is assumed in our post-processing, that the density ratio is ρ p ρ f = , which corresponds to raindrops in air. But for every Stokes number, one may consider di erently sized particles by varying the density ratio. is fact will be used in Chapter . In all cases we assured that the T . Summary of the basic DNS runs. Given are the Stokes number St, the volume fraction occupied by the particles Φ , the total time t tot , for which the trajectories have been integrated a er the initial phase of T L , and nally the total number of collisions detected in the ghost collision approximation N GCA , when assuming a density ratio of ρ p ρ f = .

St 0.0 0.10 0.20 0.30 0.51 0.76 1.01 conditions necessary for the application of the simple form of the Maxey-Riley equations ( . ) and ( . ), namely a ≪ η and ρ f ≪ ρ p , were ful lled.

Φ 0 ×
Table . also lists the special case St = , which corresponds to Lagrangian tracer particles, that follow the ow exactly. In principle those particles do have no extent, but it was necessary to assume a nite size, to be able to detect collisions. erefore in terms of the collision detection algorithm, the radius of these particles was chosen to be the same as for particles with St = . . . e DNS code described above has been used in various scienti c publications. It has been validated many times, and its results have been demonstrated to be correct.

e post-processing code however, that does the collision detection and all particle statistics, has been written for this purpose and needs proper checking. To this end we will compare di erent exemplary results to published data. We start with the RDF.

. . adial distribution function

Since Sundaram & Collins (

) rst used it to study particle collisions, many studies have been devoted to a precise determination of the RDF. In one of the earlier extensive publications, Reade & Collins ( ) give the functional form ). We nd a good agreement, although the results of Woittiez et al. (

g(r) = c ′ g, (r η) -c ′ g, exp -c ′ g, (r η) + ,
) are smaller than ours.

. . adial relative velocity

Next, we want to check our results for the RRV, which presents the second ingredient of the kinematic formulation of the collision kernel. In this case we determine the mean value directly without passing by the PDF. e latter could again be t similarly to Equation ( . ) (see Rosa et al.

).

a) g( a)

St F .

Re λ = R + Re λ = R + Re λ = W + Re λ = b) . . . . . . ⟨ w r ⟩ u K St Re λ = R + Re λ = R + Re λ = W + Re λ = F . (a)
Comparison of the dynamic and kinematic collision kernels determined within the GCA. is provides a cross-check between our collision detection scheme and the evaluation of ⟨ w r ⟩ and g( a). Also shown is data from Rosa et al.

( , R + ). Γτ K ( a) St Γ kin Γ dyn R + Re λ = R + Re λ =
We content ourselves however with comparing the mean value ⟨ w r ⟩ to published data, again from Rosa et al. ) is very good.

. . ynamic and inematic ollision kernel

Finally, to verify our collision detection scheme, we compare in Figure .

the kinematic and the dynamic collision kernel. Up to a tiny uncertainty, the two ways of determining the collision kernel, give the same result. Also the comparison to the results of Rosa et al. ( ) displays no problems. We conclude thus, that our post-processing gives the correct results and we proceed with the main scienti c questions treated in this thesis.

MULTIPLE COLLISIONS AND SHORTCOMINGS OF THE GHOST COLLISION APPROXIMATION

Section . . introduced the ghost collision approximation (GCA), which consists in counting every collision between a given pair of particles as a new collision, even though the pair may have already collided before. In this chapter we quantify the error, which this simpli cation introduces in the estimation of the collision kernel. is error will be related to multiple collisions between a same pair of particles, which in turn is closely related to a tendency of particles to stay in proximity for long times. We will present the contact time, which has been introduced in Section . , for inertial particles in turbulent ows. Finally, di erent alternative algorithms will be discussed.

. We introduce Γ GCA , the collision kernel obtained by using the GCA. e rate Γ GCA can be compared to Γ , which counts only the rst collisions occurring between a given pair of particles. e de nition of Γ can be generalized in a straightforward manner to Γ Nc , which counts the number of times that a particle collides for the N c -th time. To determine Γ Nc , we store for every pair the number of collisions it has undergone with each other. If a new collision of that pair is detected, the corresponding number of collisions of this pair is incremented as well as the number of collisions used to determine the corresponding Γ Nc .

e rates Γ Nc in fact allow us to decompose Γ GCA systematically as

Γ GCA = ∞ Nc = Γ Nc .
We argue, that in a system, where particles react upon their rst contact, e.g., a cloud where droplets coalesce (assuming unity e ciency), the collision kernel of rst contacts, Γ , is a more appropriate estimate for the "real" collision kernel. Our arguments will be presented in more detail in Section . . With these de nitions, one obtains that the GCA. To further describe this e ect, we introduce the PDF for a pair to undergo N c collisions with each other, conditioned on the fact, that the pair collides at least once. is PDF P(N c N c ≥ ) is shown in Figure . a and interestingly obeys a law

Γ GCA = Γ + Γ m , Γ m = ∞ Nc = Γ Nc . ( . )
P(N c N c ≥ ) = βα Nc , ( . )
for N c ≥ . is result can be interpreted in a Markovian sense: A er a pair has collided at least two times, it has a probability (α) to separate and not collide again, and a probability α to undergo more collisions. Together with our data Figure . a shows ts according to Equation ( . ) and the corresponding t parameter α. e probability varies with the Stokes number, such that multiple collisions are less probable for particles with larger Stokes numbers. is agrees with the ndings of the last section, where it was found that the spurious e ect of the GCA diminishes with growing St.

Combining Equations ( . ) and ( . ), we may re-express the collision rate of multiple collisions

Γ m = Γ ∞ Nc = βα Nc = Γ βα -α . ( .
)
As a test for consistency, we present in Figure . b a comparison between Equation ( . ) and the numerical results for Γ m . e di erence between the two is less than %. F . e distance between two di erent pairs of particles over the course of time, both with St = . e distance for the one pair is given as full lines, the distance for the other as dashed lines. Panel (a) shows the entire runtime of the simulation and the distance in units of the box size π. e dotted line gives the expected average distance for a pair in a periodic cube. In (b) a zoom of the data is shown. Here the distance is given in units of the collision radius, which is also highlighted by a dotted line. In both cases the time is given in units of the large eddy turnover time.

only one collision, or many collisions. It shows the distance r(t) between two particles. e dashed lines correspond to a pair of particles that collides only once, and separates right a er this collision. e full lines correspond, on the contrary, to a pair of particles that stays in proximity a er an initial collision. e distance r(t) uctuates around the collision radius of a, thus causing several events which are interpreted as di erent collisions when using the GCA. Overall, this pair stays closer than a for a time of T L , whereas the contact time of the other pair is orders of magnitude smaller.

is observation motivates the investigation of the contact time presented in the next section.

e fact that particles in turbulent ows may stay close for relatively long times has been noted before, for example by Jullien et al. . e contact time τ c has been introduced in Section . for ideal gas particles with only one consecutive collision in mind. Figure . shows however, that the situation for particles in a turbulent ow is more complicated. We therefore re ne the de nition of the contact time as illustrated in Figure . . With these de nitions, we can for example express the contact time during an i-th collision as τ c,i = t s,t e, . Similarly the time between two consecutive encounters can be expressed as t e,i+t e,i . . Illustration of our de nition of the di erent "contact" times. Shown is the distance between a pair of particles over time. When it falls below some value d c for the i-th time, we term this an encounter and name the corresponding moment t e ,i . We proceed correspondingly with the times, when the particles separate, and name them t s,i . From these de nitions di erent time intervals, like the contact time τ c ,i = t s,t e , can be derived.

a) - - - - - - - T L P(τ c , ) τ c , T L d c = a d c = a d c = a d c = a b) - - - - - - - - - - T L P(τ c , ) τ c , T L d c = a d c = a d c = a d c = a
F . e PDF of the rst contact time for di erent critical distances in semi-logarithmic (a) and doublylogarithmic (b) scaling.

We have seen that the contact times may be relatively large, therefore we express them in terms of the large eddy turnover time. Another plausible choice could be the inverse of the dominant Lyapunov exponent of the particles (Bec et al.

). ese two time scales are very similar for the Reynolds number considered here. Deciding which is the proper time scale would require simulations at di erent Reynolds number, which were not performed in the realm of this work.

Figure . shows the PDF of the rst contact time τ c, , which we have determined for di erent values of the critical distance d c . In the investigated range with d c ≪ η, the results depend only very slightly on this parameter, which will be set to d c = a in the following.

e data shown in Figure . was obtained for particles with St = . (the Stokes number dependence will be discussed later). As it was the case for the ideal gas particles, see Figure . , the PDFs exhibit a power law behavior for a range of values of τ c, .

e constraint that the PDF has to be normalizable imposes that the scaling law cannot extend all the way to τ c, = when the exponent is larger than . We suspect that the functional form changes again for very small times, as was the case for the ideal gas

a) - - - - - τ c ,i T L T L P(τ c, N c > ) T L P(τ c, ) T L P(τ c, ) T L P(τ c, ) b) - - - - - - - - τ c ,i T L T L P(τ c, N c > ) T L P(τ c, ) T L P(τ c, ) T L P(τ c, )
F . e PDFs for the second, third, and fourth contact time, as well as the PDF for the rst contact time conditioned on the fact, that the same pair will collide at least one further time.

e data is shown in both semi-logarithmic (a) and doubly-logarithmic (b) scaling. e Stokes number is St = . . particles (comp. Section . . ). We do however not resolve such small times in our DNS.

Contrary to the results obtained for the ideal gas, the PDFs in Figure .  display exponential tails at long contact times, revealing a nite probability for particles to stay close for several large eddy turnover times. is feature seems to be related to the phenomenon of multiple collisions. Further evidence for this conclusion comes from Figure . , which shows the contact time probabilities for second and higher collisions. Here the power law behavior vanishes. is behavior was interpreted as being related to ideal gas particles-that is particles, which collide once and separate right a erward.

e fact that it disappears for multiple collisions con rms our interpretation. e exponential tails for long contact times however rest, con rming again the interpretation, that those are related to multiple collisions.

A further interesting conclusion from Figure . is, that also for the rst contact time PDF, the power law vanishes, if it is conditioned on having further collisions a er the initial one. Remarkably, this conditional PDF is very similar to the PDFs of higher contact times. Figure . further illustrates this phenomenon, by comparing the usual rst contact time PDF to two conditional PDFs: Once conditioned on having one unique contact, the other time conditioned on having multiple contacts.

We go on by presenting the PDF of time in between collisions in Figure . . Evidently, the typical timescales become larger, i.e., the PDFs get shi ed to the right with respect to the contact time PDFs presented up to this point. Even the cuto at small times gets resolved in this case. ere is however no power law for intermediate times, which is in accordance with our interpretation, that it is related to ideal gas like collisions. e exponential tails remain.

T L P(t e ,i+ - t e , i ) (t e ,i+ -t e ,i ) T L i = i = i = b) - - - - - - - - T L P(t e ,i+ - t e , i ) (t e ,i+ -t e ,i ) T L i = i = i = F .
T L P(τ c , T L ) τ c , T L St = . St = . St = . St = . b) - - - - - - - - - - - - T L P(τ c , T L ) τ c , T L St = . St = . St = . St = .
F . e rst contact time PDF for di erent Stokes numbers as indicated in the legend. Alongside the data (symbols) ts according to Equation ( . ) (lines) are shown.

e data is given in semi-logarithmic (a) and doublylogarithmic (b) scaling.

. . ependence on the tokes number

Up to here, all shown contact time PDFs were for the case St = . . Figure .

shows the rst contact time PDF for di erent Stokes numbers. e qualitative appearance does not change. But the overall probability for long contact times decreases with increasing St, whereas the range of the power law extends and its slope steepens. We nd that all these PDFs can be tted by

P(τ c, ) ∼ (τ c, T L ) -ξ exp (-κτ c, T L ) , ( .
)
where the omitted proportionality constant is determined by the correct normalization. e lines in Figure . show ts of this form. ey all work very well, but the quality of the ts degrades for St < . (not shown). e two t parameters κ and ξ are shown in Figure . as a function of the Stokes number. As noted before, the slope of the power law conveyed by the parameter ξ steepens and one may suspect that it reaches ξ = for St → ∞ according to our ndings for the ideal gas particles in Section . . . e exponential law seems to vanish.

. /

In Equation ( . ) a decomposition of the collision kernel into one contribution from shear induced collisions, and another one from sling/caustics collisions, was proposed. It is interesting to ask, whether the multiple collisions we report on in this chapter, belong mainly to one of these two groups.

An evident assumption is, that multiple collisions are mainly due to shearing motion. Sling/caustics collisions happen between particles originating from di erent regions of phase space, at large RRV. It is probable that those particles will continue their trajectories a er the collision and separate quickly without further collisions. is reasoning is con rmed by the results shown

u K P c ,i ( w r ) w r u K i = GCA i = - - -
F . e PDF of the RRV obtained for all colliding particles as in the GCA (full lines) and for only rst collisions (dashed line).

e inset shows the same PDFs in a wider range and in semi-logarithmic scaling. ese results are for St = . . in Figure . , which shows the PDF of the absolute value of the RRV. is PDF was obtained once, taking into account all collisions, as in the GCA, and another time, taking into account only rst collisions. Both PDFs are sharply peaked at low RRV and have exponential tails. A close comparison however reveals, that the PDF for rst collisions lies below the other PDF for small RRV and above it for larger RRV. As the di erence between the two PDFs must be attributed to multiple collisions, on may therefore conclude, that multiple collisions happen at small RRV. A strong evidence for this conclusion is provided by Figure . . ere the ratio of the mean RRV of rst collisions, ⟨ w r ⟩ c, , and the mean RRV of all collisions, ⟨ w r ⟩ c,GCA , is shown. At small values of the Stokes number, the mean RRV of rst collisions is larger, hinting again that multiple collisions have on average small RRV. At larger Stokes numbers the di erence between the rst collision kernel Γ and the ghost collision kernel Γ GCA vanishes and so does the di erence in the mean RRVs. e inset of Figure . shows the mean RRV for only multiple collisions, in comparison to the mean RRV of rst collisions. e value for multiple collisions is much smaller than the one for rst collisions. e assumption, that multiple collision are not due to sling/caustics collisions, which happen at large RRV, but rather due to shear induced collisions is therefore fully consistent with our numerical results.

⟨ w r ⟩ c, ⟨ w r ⟩ c ,GCA St ⟨ w r ⟩ c,i u K St i = i = m .
In the introduction to this chapter it has been argued, that the collision kernel, for which only rst collisions were taken into account, was a more realistic estimate than the collision kernel obtained in the GCA.

is argumentation was based on a system, where particles react upon their rst collision. To be more precise one should de ne, what kind of reaction the particles undergo. Krstulovic et al. ( ) investigate a system, where the unique species of particles Λ reacts with each other by disappearing, or in other terms Λ + Λ = ∅. In this case Γ is certainly a better estimate, but it is not clear, if it is the correct estimate. In such a system it is not only impossible that a pair collides with each other a second time. Actually each of the particles is not allowed to collide with any further particle a er an initial collision. To investigate such a system, we propose to mimic it with the help of the algorithm described in the following section.

. . eplacement algorithm

It has been mentioned before that we usually integrate a total of N ′ p particles in our simulations. Nevertheless, we restrict ourselves to a subset of N p particles to check for collisions. erefore we may, at each collision, replace the two colliding particles by two of the additional N ′ p -N p particles. Replacing the particles instead of just "deleting" them, has the advantage, that the system stays in a stationary state with respect to the number of particles N p . We will refer to this algorithm, that replaces both colliding particles, as R and to the corresponding collision kernel as Γ R .

A second similar algorithm consists in replacing only one arbitrarily chosen particle of each colliding pair. We will refer to this algorithm and the corresponding collision kernel as R and Γ R respectively. is algorithm could be seen as representing a system, where the particles react such that Λ + Λ = Λ. Of course one would need to take into account other e ects like momentum conservation in a more realistic model.

Γτ K ( a) St Γ GCA Γ Γ R Γ R F .
Here the same data as in Figure . a is shown. Additionally the collision kernels Γ R and Γ R that were obtained with the replacement algorithm described in Section . . are shown.

We nally note that for tracer particles, it is not necessary to integrate additional trajectories. In this case two new particles may be placed randomly in the ow and their velocities are fully determined by the uid velocity eld at that instant. For inertial particles however, the momentary velocity depends on their history and therefore it is necessary to keep the additional particles in the ow. ) investigate particles with nite inertia. ere, they use a further algorithm. ey simply remove all colliding particles and take the diminishing number of particles into account when calculating the collision kernel.

. . ynamic collision kernel with the alternative algorithms

We have implemented the alternative algorithms described in the last section and calculated the corresponding collision kernels. Our results are shown in Figure . . As expected, we nd that the results obtained with the new algorithms R and R lie below those obtained within the GCA. However, they lie even below the collision kernel for only rst collisions, which seems puzzling at rst. Before explaining the observed di erence between Γ and Γ Ri , i = , , we note that also Γ R and Γ R di er from each other. When we replace both colliding particles, fewer collisions take place.

To explain the di erence between the di erent algorithms, consider the following situation: Two groups ("jets") of particles, originating from di erent regions in phase space, are brought into proximity by the sling/caustics e ect. In the GCA, as the two jets cross each other, any one particle from one group, may collide with several other particles from the other group.

e same is true, of course, for any particle with which that particle collided. Now in this special case, Γ and Γ GCA would not di er at all, because none of the multiple collisions that any particle undergoes takes place with a particle it collided with before. e collision kernels obtained in the two new algorithms however will di er from Γ GCA . With the algorithm R only one of the colliding pairs has the chance to collide with further particles, explaining, why Γ R is smaller than Γ . e other algorithm, R , replaces both Γ R (Φ )

Γ GCA (Φ ) Γ (Φ ) Γ m (Φ ) Γ R (Φ ) Γ R (Φ )
F . e collision kernel obtained in the GCA, Γ GCA ; the collision kernel for only rst collisions, Γ ; the collision kernel for only multiple collisions, Γ m = Γ GCA -Γ ; and the collision kernels for the replacement algorithms introduced in Section . . , Γ R and Γ R . All these collision kernels are shown for two di erent volume fractions, once for Φ as given in Table colliding partners, therefore none of them will participate in any further collisions. at is why Γ R is even smaller than Γ R .

e above interpretation relies on some very particular hypothetical assumption, namely that there are events in a turbulent ow, where jets of particles collide with each other.

e assumption used here suggests in fact that the algorithms R and R a ect the number of collisions, insofar as groups of three particles (or more) are concerned; namely the two particles that initially collide and at least one other particle. If this is the case, the obtained collision kernel should depend on the number density n of the particles. Remember, that the rate of collisions is proportional to n . is proportionality is broken, if we suppress some collisions with a rate proportional to n . erefore the collision kernel, which was de ned as precisely the proportionality constant joining the rate of collisions and n , will depend on n.

To check our reasoning we performed the same analysis, that led to the results shown in Figure . , for half the number density as before. In other terms, we reduced the number of particles considered in the collision detection scheme by a factor of two, resulting in a halved volume fraction Φ .

e results are shown in Figure . . First of all, one notices, that as expected neither Γ GCA , nor Γ , nor Γ m depend on the volume fraction Φ. e collision

. . . . St Γ R ,R (Φ ) Γ Γ R ,R (Φ ) Γ Γ R ,R (Φ ) Γ F
. e ratios Γ Ri Γ , i = , for three di erent volume fractions and for di erent Stokes numbers are represented by symbols as indicated in the legend. We excluded the two points, where according to Table .  Φ is smaller than for the rest. e full lines are meant as a guide for the eye. ey are chosen, such that the distance between these curves and falls by a factor of two for the next higher line. kernels obtained with the alternative algorithms, Γ R and Γ R , however do change, when we vary Φ. is is in accordance with our reasoning.

Interestingly, looking closely one realizes, that the collision kernel obtained by replacing both particles coincides with the collision kernel obtained by replacing only one particle at double the volume fraction. is observation as well is in accordance with the above argumentation, that the discrepancy between the Γ Ri , i = , , and Γ stems from e ects involving three and more particles. In fact, according to our reasoning, the di erence Γ -Γ Ri , i = , , should be to leading order linear in n and therefore in Φ.

is is investigated in more detail in Figure . , which shows the ratios Γ Ri Γ , i = , for three di erent volume fractions all within a factor two of each other. We empirically nd the two rules

Γ R (Φ ) = Γ R (Φ) and Γ Ri -Γ ∼ Φ for i = , . ( . )
e latter seems to be only approximately valid, but both con rm our argumentation.

To make our explanation more clear, let us present it once more, but from another perspective. We try to explain, why Γ R and Γ R are smaller than Γ and why the former two vary with the volume fraction Φ. To this end, let us follow a particle along its trajectory. e rate of collision along this trajectory as been introduced in Section . . . It is proportional to n, and therefore to Φ, namely = Γn. is is exactly true for collisions in the GCA, but also when counting only rst collisions. e two alternative algorithms introduced in Section . however necessitate a correction. When the particle we follow collides, the algorithm R changes its position and its (former) surroundings, by replacing also the colliding partner.

e algorithm R changes only one of these, either the surroundings, or the position of the particle. Both approaches certainly prevent any further collisions between the two colliding particles. But they also suppress collisions between the test particle and a third particle, at least for the time, the ow needs to reorganize a er the "perturbation" of replacing one or two particles. During this time, which we shall call τ R , a part of the collisions, which would take place at a rate , are suppressed. With respect to the collision kernel, that counts a) only rst collisions, we therefore expect a correction of the order of τ R , or more precisely

Γτ K ( a) St Γ GCA,kin Γ GCA,dyn Γ R ,kin Γ R ,dyn Γ R ,kin Γ R ,dyn b) - - g(r) r a St = . St = . St = .
Γ Ri = Γ [ -O(i τ R )] . ( . )
As mentioned above, the rate is proportional to Φ and therefore Equation ( . ) ful lls the two empirically found laws ( . ). e above explanation of the discrepancy between the collision kernel Γ and the two collision kernels obtained with the alternative algorithm, Γ R and Γ R , implies that this di erence is proportional to the volume fraction Φ. As a consequence, the discrepancy must vanish in the very dilute limit, i.e., Φ → . erefore, in this limit the two algorithms R and R give the correct result.

. . inematic collision kernel with the alternative algorithms

One disadvantage of the replacement algorithm is illustrated by Figure . a. e collision kernel cannot be determined kinematically, that is with the help of Equation ( . ), anymore. For the GCA, the dynamic collision kernel, determined by counting collisions, and the kinematic kernel coincide. For both alternative algorithms, R and R , however, the kinematic collision kernel is too small. e reason for this becomes clear in It would be desirable to have a better understanding of how the replacement algorithm a ects the estimation of g( a) and ⟨ w r ⟩. At the moment we cannot provide a conclusive explanation. F . e RDF (a) as well as the mean RRV (b) at contact. Shown is data obtained in the GCA and in the alternative algorithms from Section . . . We have discussed two striking e ects, that can be observed for inertial particles in turbulent ows.

ose particles have a tendency to stay in proximity for very long times of the order of the large eddy turnover time.

is e ect had been reported for tracer particles, but we have shown, that it remains valid for particles with inertia. It slowly decays with growing Stokes number. We presented a thorough study, especially in terms of the contact time PDF. e latter shows some features, which are related to the motion of ideal gas-like particles. But additionally, the contact time PDF shows exponential tails for long contact times.

e fact, that the particles have a certain probability to stay close for long times, has as a consequence that multiple collisions between these particles take place. is e ect was studied as well and the corresponding PDF of the number of collisions per pair was discussed. It displays exponential tails leading to a simple modeling of the collision rate for multiple collisions.

ese multiple collisions lead to an overestimation of the collision kernel in the GCA. An improved estimate, namely the collision kernel for only rst collisions Γ was proposed. Determining the latter is however a tedious task. erefore an alternative algorithm was proposed, which delivers the correct result in the limit of very dilute suspensions.

An important aim for future work is to gain a better understanding of the multiple collision events. One can imagine that these events are related to coherent structures in the ow. For instance, the particles which stay together for a long time could be trapped in a strong vortex tube. Answering these questions would necessitate the knowledge of typical uid parameters, like the uid velocity or vorticity, along the particle trajectories. is has not yet been implemented, but is in principle feasible and can be envisaged for subsequent work.

In Section . we compared, as much as feasible, our results to those of Andersson et al. (

) and Gustavsson & Mehlig ( a). is theoretical work is based on simplifying assumptions (small Kubo numbers), and it has only been tested in Gaussian random ows. It would be interesting to test the validity of the theoretical approach, at small Kubo numbers, with the help of DNS.

Recently Gustavsson & Mehlig ( a, b) have developed a theory which describes the PDF of the RRV. At the moment, we do only determine this PDF for colliding particles, as presented in Figure . . Applying techniques along the lines of Section . . to obtain the PDF for colliding particles from the results of Gustavsson & Mehlig ( a, b) proves nontrivial. But it is feasible to calculate exactly those quantities for which their theory holds in our simulations. is necessitates the implementation of additional features, which can be envisaged as a future project.

Let us conclude this section by noticing that several signi cant approximations have been used here, and that it would be important to quantify the errors they induce. For cloud physics applications, the in uence of gravity on the e ects described in this chapter would need to be studied

(Grabowski & Vaillancourt
). Also, the role of polydispersity is known to be very important, so similar studies should be carried out for bidisperse solutions with a wide range of particle size di erences. e particles that stay together are those which simply follow the ow. We know that for Stokes numbers of moderate size, a fraction of the particle population essentially follows the ow. erefore we do expect to observe the e ects described in this chapter for polydisperse suspensions as well. But they remain to be quanti ed.

Finally the role of hydrodynamic interactions (see Section . . ) should be studied. is would demand the implementation of slightly more complicated models, like the one by Ayala et al. (

) and Wang et al. ( a). e most interesting question in this realm is probably, in how far the phenomenon of multiple collisions a ects the collision e ciency (introduced in Section . . ).

PREVALENCE OF THE SLING/CAUSTICS EFFECT

In Section . . we presented the decomposition of the collision kernel

Γ SC = π( a) ⟨ w r ⟩g( a), ( . )
which was rigorously derived by Sundaram & Collins ( ). It has been demonstrated to give the correct results (see also Figure . ). e aim of introducing this composition was to study the in uence of the preferential concentration e ect on the collision rate. But the interpretation of Equation ( . ) is actually more subtle, and the aim of this section is to reconsider it in light of our own numerical results.

We begin by noting that, of course, by the introducing of the RDF, the clustering e ect is correctly taken into account by Equation ( . ). Also, the RDF at contact g( a) grows with the Stokes number, as can be seen in, e.g., Figure . a. erefore the clustering of particles enhances the collision rate to some extent. But at the same time also the mean RRV ⟨ w r ⟩ increases (Figure . b) thereby additionally enhancing the collision rate. e increase of the mean RRV is due to the presence of a further e ect, which is-like the preferential concentration-not present at St = , namely the sling/caustics e ect as presented in Section . . . is e ect was rst described by Falkovich et al. ). Equation ( . ), while correctly describing the collision rate in the GCA, does not distinguish the two e ects, and simply rests on an "e ective theory". It is therefore di cult to evaluate precisely which of the two e ects accounts for the enhancement of the collision rate as a function of the Stokes number. In this chapter, we demonstrate that our numerical data show very clearly that the sling/caustics contribution provides the dominant mechanism enhancing the collision rate, even at relatively small Stokes numbers. We also show that our results can be best understood with a representation of the collision rate as introduced in Equation ( . ).

e decomposition ( . ) describes the collision rate in the GCA, hence all collision kernels presented in this chapter are obtained in this approximation and su er from the insu ciencies discussed in Chapter . is fact will however not a ect the conclusions, which are based on qualitative features of the results, rather than on exact quantitative values.

Γ = Γ ST g( a) + Γ A h S (St, Re λ ) ( . )
from Section . . [there printed as Equation ( . )]. It states that the total collision kernel in a turbulent ow is a combination of Sa man & Turner's ( ) and Abrahamson's ( ) collision kernels. It is interesting to note that the rst term, which is due to shearing motion and a ects particles that follow the ow, scales like Γ ST ∼ ( a) τ K -see Equation ( . ). e other term however, which describes particles that are decorrelated from the uid motion and may collide with large relative velocities, scales like Γ A ∼ ( a) η τ K -see Equation ( . ). Both terms have a prefactor, but only one of them depends on the particle radius a. It has been shown in Section . . that g(r) ∼ r -c g, .

erefore the entire rst part scales like Γ ST g( a) ∼ ( a) -c g, , which is always di erent from Γ A ∼ ( a) as c g, < (see Figure . b). Hence both terms on the right-hand side of Equation ( . ) behave di erently when the size of the particles is changed.

e dimensionless Stokes number, St = (ρ p ρ f )(a η) , introduced in Section . . , is the crucial parameter that speci es the dynamics of the particles. is means that the particle size a can be changed without a ecting the value of the Stokes number, if at the same time the density ratio ρ p ρ f is varied. We used this fact to investigate the scaling with regard to the particle radius a of the collision kernel in our DNS. To this end we post-processed our data as described in Chapters and , but here under the assumption of three di erent density ratios.

ese ratios were chosen such that the according particle radii are within a factor two of each other.

e collision kernel obtained for the three di erent particle radii is shown in Figure . . Both scalings, the one according to Γ ST and the one according to Γ A , are presented. Sa man & Turner ( ) obtained their result for St → and in this limit the scaling of their collision kernel, ∼ ( a) τ K , works perfectly. But already for relatively small Stokes numbers our results start to diverge, demonstrating that the scaling is not valid anymore. e scaling which one expects for sling/caustics collisions is correct for larger Stokes number, as can be seen from the good correspondence of our results in From Figure . it becomes clear that for St → the shear induced collisions describe the overall collision kernel correctly and that for larger Stokes numbers the sling/caustics e ect is more important. But the in uence of the preferential concentration e ect described by g( a) is not evident. erefore we plotted in Figure . the ratio of the collision kernels obtained for two di erent density ratios. In each case the particle radii di er by a factor of two. Hence one would expect a value of ⁄ , if ( a) τ K was the correct scaling and a value of ⁄ for ( a) η τ K . e entire rst term of Equation ( . ) including g( a) leads to a value of ( ) -c g, . Figure . shows these three

ρ p ρ f ) Γ(ρ p ρ f ) St Γ( ) Γ( ) Γ( ) Γ( ) ( ) -cg,

F

. Ratios of the collision kernels of particles with the same Stokes number, but whose radii di er by a factor two. If the collisions according to the description by Sa man & Turner ( ) dominate, one expects a value of . ; else, if most collisions correspond to the description by Abrahamson ( ), the expected value is . . We compare raindrop-like particles (ρ p ρ f =

) once to two times larger particles (ρ p ρ f =

) and once to two times smaller particles (ρ p ρ f =

). e dashed line represents the expected scaling of the rst term on the right-hand side of Equation ( . ).

values as well as the ratios of the collision kernels we obtained in our DNS. For St < . the numerical results follow the dashed line given by ( ) -c g, .

is means that in this range the enhancement of the collision kernel is almost entirely due to the e ect of preferential concentration. At larger Stokes numbers however, the dashed line and the numerical results start to diverge. For St ≳ . the curve ( ) -c g, has reached its maximum and starts to fall. Nevertheless the numerical results continue to grow and reach the value ⁄ , which is expected for sling/caustics collisions. is leads to the conclusion that at larger values of the Stokes number the majority of collisions is due to the sling/caustics e ect.

F . e cumulative PDF of the RRV F( w r ) as de ned in Equation ( . ) and the contribution of RRVs smaller than w r to the ux entering the collision sphere φ( w r ) de ned in Equation ( . ). Both are shown for St = . and in two di erent scalings: Once in units of the Kolmogorov velocity u K (bottom axis) and once in units of the shear rate at the particle size a τ K (top axis). . To provide further evidence for the prevalence of the sling/caustics e ect, we calculated the cumulative PDF of the RRV, de ned as

F( w r ) = ∫ wr p( w ′ r ) dw ′ r . ( . )
e cumulative PDF F( w r ) allows us to see, which values of w r are the most probable. For instance, one can read from Figure . that for St = . only about % of the particles have relative velocities larger than w r ≳ a τ K .

In addition to F( w r ), we calculated the ratio of particles entering the collision sphere with RRV smaller than w r

φ( w r ) = ∫ wr w ′ r ⟨ w r ⟩ p( w ′ r ) dw ′ r . ( . )
is can be interpreted as the contribution of velocities smaller than w r to the total collision rate. Along the lines of Sections . . and . . , φ( w r ) can also be seen as the cumulative PDF of the RRV conditioned on collisions. Figure . shows that at St = . , although rare, particles with w r ≳ a τ K contribute about % of the collisions. e collisions with larger RRV are likely to be sling/caustics collisions. erefore, the numerical values presented here provide a strong evidence that the sling/caustics e ect starts to become the leading mechanism enhancing the collision rates already at relatively small Stokes numbers.

Figure . shows F( w r ) and φ( w r ) for two further Stokes numbers, St = . and St = . . By comparing the two, the transition from small to large Stokes values becomes evident-note especially the change of scale in the horizontal axes. At small values of St shear induced collisions prevail. Almost all particles have very low RRV and almost all collisions take place at such low RRV. For increasing Stokes numbers however, more particles move at larger RRV and even more collisions are contributed by these particles. F . e same quantities as in Figure . , but here in (a) for St = . and in (b) for St = . .

g( a)

St

ρ p ρ f = ρ p ρ f = ρ p ρ f = F
. e RDF at contact g( a) for di erently sized particles (as tuned by the density ratio ρ p ρ f ). e smaller the particles, the larger g( a).

. In Section . the prevalence of the sling/caustics e ect for particle collisions in turbulent ows was deduced from an analysis of the scaling dependence of the collision rate as a function of the particle radius a. To this end numerical simulations of particles with the same Stokes number but varying radii were carried out. Looking again at past numerical and theoretical results, we explain in this section that our conclusions are in fact fully supported by numerical evidence obtained in previous works.

Equation ( . ) reduces the collision kernel to two quantities which have been the subject of many studies, namely g( a) and ⟨ w r ⟩. Both quantities depend on the particle radius a as can be seen in Figures . and. , but in di erent ways. While the RDF becomes larger with decreasing particle radius, the RRV grows with increasing radius. According to Equation ( . ) these opposing behaviors must cancel exactly for larger Stokes values, in order to give the correct scaling of the collision kernel ∼ a . Figure .  demonstrates that this is indeed the case.

A similar observation allows to con rm our results concerning the scaling of the collision kernel by comparison with published data. e RDF g(r)

ρ p ρ f = ρ p ρ f = ρ p ρ f = b) - - - ⟨ w r ⟩ u K St ρ p ρ f = ρ p ρ f = ρ p ρ f = F
. e RRV at contact for particles of di erent radii, shown in panel (a) in linear scaling and in panel (b) in logarithmic scaling. e radius was changed at constant Stokes number by varying the density ratio ρ p ρ f . For growing radii, the RRV increases.

ρ p ρ f = ρ p ρ f = ρ p ρ f = b) - - - . . . g( a)⟨ w r ⟩ u K St ρ p ρ f = ρ p ρ f = ρ p ρ f = F .
e RRV multiplied by the RDF at contact. Both scale di erently with the particle radius as can be seen in Figures . and . obtained in our DNS was shown in Figure . . For small distances r ≪ η it is described very well by a power-law g(r) ∼ r -c g, . erefore g( a) must scale like ( a) -c g, . Up to now, only the RRV at contact ⟨ w r ⟩ was discussed in this work. But it is a known fact, that for small separations the RRV scales with the distance at which it is calculated as ⟨ w r (r) ⟩ ∼ r +c w , (e.g., Bec et al.

). eories to predict this power law behavior have been proposed by Gustavsson & Mehlig ( a, b), who studied in fact a simpler random ow model. As a consequence the collision kernel must scale with the particle radius as Γ SC ∼ a -c g, +c w , .

Rosa et al. (

) present a thorough study of the collision kernel in DNS. Among others, they show data for c g, and c w, . From that data we obtained the scaling of the collision kernelc g, + c w, , which we present in .

We investigated here collision rates in turbulent ows. Two e ects are known to lead to a strong enhancement, with respect to the Sa man-Turner prediction, of the collision rate. On the one hand shear induces collisions of particles which are already in proximity. is e ect is enhanced by preferential concentration which augments the particle density in certain regions.

On the other hand the sling/caustics e ect describes how particles, which come from di erent regions of phase space, may collide with large RRV. We have shown that for small Stokes numbers an enhancement of the collision rate is due to preferential concentration. But as the Stokes number increases the sling/caustics e ect becomes more and more important, whereas the e ect of preferential concentration looses in uence. We could identify the region St ≳ . as the one where the majority of collisions is due to the sling/caustics e ect.

Our conclusions rely on two arguments. First, the scaling of the collision kernel with the particle radius is di erent for the two e ects. But for St ≳ . our numerical results reach the value imposed by the sling/caustics e ect. Secondly, the cumulative PDF of the RRV F( w r ) and a derived quantity φ( w r ) allowed us to show that for growing Stokes numbers an increasing ratio of collisions happens at large RRV. ese collisions are induced by the sling/caustics e ect.

As was the case in Chapter , it would be interesting to study the in uence of the di erent approximations that lead to the model investigated numerically. In particular, neglecting gravity tends to "arti cially" enhance the collision rate (see, e.g., Rosa et al. ). erefore an important question would be, which of the two di erent contributions is more a ected by the inclusion of gravity. Also the investigation of polydisperse solutions could provide interesting insights. Generally, the collision kernel Γ(a, a ′ ) is a complicated function of both variables, especially in the presence of gravity (see, e.g., Grabowski & Wang ; Woittiez et al.

). It is therefore di cult to predict, whether the sling/caustics e ect is of similar importance in bidisperse solutions-a thorough investigation would be necessary.

KINEMATIC SIMULATIONS

When we started our investigations of the collision rates of particles in turbulent ows, we did not yet use DNS. Instead we used a simple model ow introduced by Fung et al. (

) known as kinematic simulations (KS). We learned much from this model, as it allowed us to perform simulations with a very low numerical e ort. Also, essentially all the qualitative information we obtained with KS turned out to be in full agreement with the results we obtained with DNS. However, in quantitative terms, we nd a large di erence between the two approaches. It is the objective of this chapter to discuss the similarities and di erences in the KS and DNS problems. e essential conclusion of our comparison is that the sling/caustics e ect, as well as preferential concentration, play a much reduced role in KS as compared to DNS. In the light of these observations, studying collisions using KS instead of exact DNS leads to results which should be interpreted with care.

In this chapter we shortly introduce KS before presenting results analogous to those from Chapters and . We will focus especially on di erences in comparison to the results obtained with DNS. Although a complete understanding of the origin of this di erence is still missing, we will present several remarks, aimed at explaining the observed discrepancies.

. e basic idea of KS is to obtain the uid velocity eld u(x, t) as a sum of N k random Fourier modes

u(x, t) = N k n= a n cos(k n ⋅ x + ω n t) + b n sin(k n ⋅ x + ω n t).
( . )

ese modes k n and the corresponding coe cients a n and b n are chosen such that the resulting ow ful lls certain properties.

e wave vectors k n = k n kn are pointing in random directions kn . To assure incompressibility the directions of the coe cients need to ful ll ân ⋅ kn = bn ⋅ kn = .

e amplitudes of these coe cients govern the energy spectrum of the ow

E(k n )∆k n = a n = b n . We impose a spectrum E(k n ) = E k - n
and de ne the discrete di erences between the wave vectors

∆k n = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (k -k ) , n = (k n+ -k n-) , n ∈ [ , N k -] (k N k -k N k -) , n = N k .
Di erent choices of the energy spectrum, as discussed by Malik & Vassilicos ( ), are possible. For the amplitudes of the wave vectors there are two possibilities. We choose a geometric distribution

k n = k L η (n-) (N k -) , k = π L , k N k = π η . ( . )
e alternative would be an algebraic distribution as discussed by, e.g., Fung & Vassilicos (

). In Equation ( . ) L denotes a large and η a small length scale. ey may be interpreted as the integral and the Kolmogorov length respectively [note however, that L does not ful ll Equation ( . )]. We have L η = in all simulations discussed here. Malik & Vassilicos ( ) present a comparison between KS and DNS. ey obtain a good agreement for a total number of modes N k ≳

. erefore we chose N k = which assures modes in each band [k, k].

e frequencies ω n in Equation ( . ) are de ned to be proportional to the eddy turnover time on the corresponding length scale

ω n = λ k n E(k n ).
e unsteadiness parameter in our simulations was λ = . . is value is expected to give results similar to those obtained from DNS concerning Lagrangian dispersion (Malik & Vassilicos ; Nicolleau & ElMaihy ). All parameters are determined in the initial phase of each simulation. ey are not changed for the rest of the run time. erefore the "randomness" in each run is limited and results do vary with the initial state of the random number generator (RNG). We use the well-known Mersenne twister (Matsumoto & Nishimura ) as RNG. For each set of parameters we perform runs initializing the RNG with a di erent seed. We obtain our nal results by averaging over the individual runs. ). As a rule of thumb we use for the time-step ∆t = min(τ p , τ K ). We also integrate the particle velocity gradient tensor σ i j as de ned in Section . . with the same scheme. is quantity may diverge, therefore its integration demands some special care.

e technical details of this procedure are described in Appendix B. In KS the system size L sys can be an arbitrary multiple of L. In our simulations it is L sys = L. Specifying the system size and a volume fraction Φ determines the number of particles N p at a certain Stokes number. Due to numerical constraints, we divided the range of Stokes numbers we investigated in DNS into two overlapping regions. In the range St ≤ . we performed simulations with Φ = . × -. In the range St ≥ . we did another set of simulations with Φ = . × -. e smaller volume fraction is represented by lled symbols in all graphics in this chapter. Empty symbols, if they correspond to data obtained in KS, represent values from simulations with the higher volume fraction. We used again the alternative algorithm R introduced in Section . . . In the corresponding simulations % additional particle trajectories were integrated.

. We determined the di erent collision kernels de ned in Chapter for KS. e results are shown in Figure . a. First of all one notes that the general aspect is the same in KS as it was in DNS. Of all the collision kernels Γ GCA is the largest. e two corrected collision kernels Γ and Γ R are smaller than Γ GCA . In contrast to our DNS results, here Γ and Γ R coincide although the volume fraction is even higher. In Chapter we had shown that Γ R was smaller than Γ by a factor proportional to the volume fraction Φ. A further di erence becomes obvious, when one compares the quantitative values.

e collision kernels we found in DNS were larger by more than one order of magnitude for St ≳ . . For St → the prediction of Sa man & Turner (

) is correctly reproduced in both approaches. b, but here for KS. Again the qualitative features are similar, but here we nd a maximal error introduced by multiple collisions which is almost %. For comparison, in DNS the maximal error was around %. We conclude that while the overall number of collisions is reduced in KS, the ratio of multiple collisions is higher than in DNS.

e reasons for the large quantitative di erences between KS and DNS become a bit clearer from Figure . . Further evidence is provided to some extent by Figure . . It shows the ratio of the mean collision velocity for only rst collisions and for all collisions in the GCA. e corresponding data was shown for DNS in Figure . . Because multiple collisions are so rare for larger Stokes values and because the statistics is dominated by sling/caustics collisions, this ratio drops very quickly in the case of DNS. For KS however, this ratio stays signi cant even for larger values of the Stokes number. is underlines again the suppression of sling/caustics collisions in KS.

e reduced preferential concentration provides also an explanation why we nd in Figure . a that Γ R ≈ Γ .

e di erence between these two collisions kernels was shown to be proportional to the volume fraction in Chapter . More precisely it depends on the e ective volume fraction in the region of a particle, namely g(r)Φ. erefore, if g(r) is smaller in KS, as was shown in

. . . . u K P c ,i ( w r ) w r u K i = GCA i = - - - .
F . e PDF of the RRV conditioned on the fact that two particles collide. To determine this PDF once all collisions detected in the GCA were taken into account (full lines). Another time only rst collisions of a same pair were considered to determine the statistics (dashed lines). e data is for St = . and was obtained in KS.

is gure may be compared to Figure . where the same quantity is shown for DNS.

peaks at a much smaller velocity and the exponential tails are shorter by an order of magnitude. is provides further con rmation that there are fewer sling/caustics collisions in KS.

. e presence of multiple collisions in KS has become obvious in the previous section. Here we present the statistics of these multiple collisions. In Figure . we show the PDF of undergoing N c collisions a er an initial one in KS, P(N c N c ≥ ). e same PDF has been shown for DNS in Figure . a. We nd again, that this PDF can be described by a function P(N c N c ≥ ) = βα Nc , although there are larger deviations in the tails as was the case for DNS. For KS we nd in general a larger probability for multiple collisions. is is especially evident in the values of the t coe cient α which is shown in the inset of Figure . . It is roughly two times larger than in DNS. e high probability for multiple collisions underlines again the fact that sling/caustics collisions are less probable in KS. e sling/caustics collisions do not contribute to the multiple collisions and would therefore diminish their probability. e di erent lines show ts according to Equation ( . ). In the inset the t coe cient α is shown.

ere the full and empty symbols have again the meaning as explained in

T L P(τ c , T L ) τ c , T L St = . St = . St = . St = . b) - - - - - - - - - - - T L P(τ c , T L ) τ c , T L St = . St = . St = . St = .

. . ontact time

We also determined the contact time PDF as de ned in Section . for KS. In Section . it was explained that this feature of the PDF is mainly due to sling/caustics collisions. erefore the absence of the power law emphasizes again the lack of this e ect in KS. Furthermore it is interesting to note that in KS our temporal resolution su ces to see that the PDF reaches τ c, = with a negative slope.

It is shown in

.

/ - - - -. . . PDF r × rσ r w r St ↑ GCA R
F . e PDF of the ratio ar ⋅ σ ⋅ r w r for St ∈ { . , . , . , . } as obtained in our KS. e peak at corresponds to shear induced collisions, the peak at corresponds to sling/caustics collisions. We show the PDF obtained in the GCA as well as with the alternative algorithm R .

.

/ -

Typically the particle velocity gradient tensor σ i j can be used to obtain a good estimate for the RRV w r (r) ≈ r ⋅ σ ⋅ r as long as the separation r is small enough. If two particles however collide due to the sling/caustics e ect, they come from di erent regions in phase space. is has been explained in Section . . . eir velocity gradients are di erent and lead to an estimate that is much smaller than the real RRV between them. is provides a means to decide at every collision, whether it is a shear induced or a sling/caustics collision. One simply has to compare the estimate with the actually measured RRV, i.e., calculate ar ⋅ σ ⋅ r w r . If this ratio is close to one, the collision is induced by shear. If the ratio is close to zero, this means that the actual collision velocity is much larger than the estimate suggests and the collision is therefore a result of the sling/caustics e ect. Instead of doing this for every collision it is of course simpler to determine the PDF of the ratio ar ⋅ σ ⋅ r w r . Ducasse & Pumir ( ) present this PDF in a detailed numerical study of collision rates in KS. Here we use the same quantity to reason once again that the spurious e ect of multiple collisions is due to shear induced collisions.

In Figure . we show the PDF of the ratio ar ⋅ σ ⋅ r w r as in our KS for several Stokes numbers. We estimated the e ective gradient tensor σ as the average of the two individual particle velocity gradient tensors σ and σ , i.e., σ = (σ + σ ). e PDF displays two peaks: One at ar ⋅ σ ⋅ r w r = which corresponds to shear induced collisions; another broader one around zero which corresponds to sling/caustics collisions. e sling/caustics e ect becomes more important as the Stokes number grows. is has also been found by Ducasse & Pumir ( ). Furthermore Figure . displays for every Stokes number the ratio as determined in the GCA and as determined in the alternative algorithm R . e di erence between the two sheds light on the nature of multiple collisions which are not present in R . e peak at diminishes when the alternative algorithm is used, while the peak at increases. erefore the multiple collisions correspond mainly to shear induced collisions. is is in accordance with the results we found before in Chapter .

We did not determine the same quantity in DNS. In the light of our results from Section . , it would however be very interesting to compare the ratio ar ⋅ σ ⋅ r w r in DNS and KS. We concluded in Chapter that the sling/caustics e ect is the dominant contribution to the collision rate in DNS, even at relatively small Stokes numbers. In Section . it was demonstrated that in comparison to DNS there are very few sling/caustics collisions in KS.

erefore one should expect that the PDF from Figure . would show a very pronounced peak around zero in DNS even for relatively small Stokes numbers.

. In this section we presented results obtained in a simple, but established, model for turbulent ows, namely kinematic simulations. On a general perspective these results con rm those from DNS which we presented earlier. But we nd a striking discrepancy in quantitative terms. All aspects we investigated con rm an absence or at least a strong suppression of the sling/caustics e ect in KS. Also the preferential concentration is much less pronounced. e reason for this could be that the modes in Equation ( . ) are chosen randomly. erefore it is very unlikely that one nds in KS the coherent structures which are so typical for real turbulent ows. According to the sling interpretation these structures are responsible for the sling/caustics e ect. ey can also be seen as being at the origin of the preferential concentration e ect. erefore the absence of coherent structures in KS would ultimately result in the suppression of the two e ects. Future investigations may shed more light on these relationships.

Part III

APPENDICES

A

INTERPOLATION OF PARTICLE TRAJECTORIES

In Section . . two di erent ways to detect a collision based on the positions and velocities of two particles at consecutive time-steps were introduced. Here, we describe both techniques in more detail.

. How can one determine that two particles in a distance r = X -X moving in relative motion to each other with velocity w = V -V eventually collide? To answer this question consider Figure A. . ere, two such particles with radius a are shown. In the following we will consider this situation in the frame of reference where the rst particle is at rest. A necessary condition is is the usual approach also used in the well-known paper by Sa man & Turner ( )-see Section . . . for a collision to take place is that the minimal distance r min between the centers of the particles (given the relative velocity w is constant) is smaller than two times their radius, r min < a. e minimal distance is the distance between the center of the rst particle P and its orthogonal projection onto the hypothetical relative trajectory of the second particle M. Initially the distance P M is given by the dot product between the relative distance and the relative velocity P M = r ⋅ w w . So the minimal distance can be obtained by the help of the Pythagorean theorem . ) De ning δ = w [( a)r min ], Equation (A. ) becomes simply δ > .

r min = r - (r ⋅ w) w ! < ( a) . ( A
To introduce a su cient condition, we have to ask on which timescale τ our hypothesis of a constant relative velocity is satis ed. In numerical simulations this timescale can usually be identi ed with the time-step. A collision takes place, if the relative distance between the particles falls below a within a time τ < τ given the particles move constantly with relative velocity w. Taking a look at where the particle distance is smaller than a for the rst time. To reach this point the second particle has to move a distance y. is distance is equal to

y = P M -CM = (r ⋅ w) w -( a) -x = r ⋅ w w - √ δ w ,
where the Pythagorean theorem was used on the triangle (P , C, M). Finally the time to travel this distance is given by τ = y w , which has to be compared to the "minimal" timescale τ. If however the collision has already happened in the past, the distance y will be negative and so will be the time τ. erefore to avoid double counting, i.e., counting the same collision two times in succeeding time-steps, one should also check if τ > .

. e method discussed in the above section is rst order in the sense that the particle trajectories are linearly extrapolated in time to detect if they collide within the (up)coming time-step. It is convenient to use, because it only necessitates the knowledge of the particles' position and velocity at one time-step. Unfortunately, as has been explained in Section . . , this comes not without cost. at is why in this section we present a higher order method.

Using such an interpolation scheme is numerically however relatively costly, therefore it is bene cial to decide early, whether it is necessary to calculate it. To this end Wang et al. ( a) de ne three di erent types of collisions, that can happen in between two time-steps, depending on the initial and the nal particle distance, r(t (n) ) and r(t (n+ ) ) respectively, where r(t) = r(t) .

. e easiest case is the one, when the particle distance was initially greater than the particle diameter, but is smaller at the end of the time-step, i.e., r(t (n) ) > a and r(t (n+ ) ) ≤ a. In this case, which Wang et al. ( a) refer to as type I, (at least) one collision evidently took place.

. Type II collisions are those, when the initial and nal distance are both greater than the particle diameter, r(t (n) ) > a and r(t (n+ ) ) > a, but the distance becomes smaller than a at an intermediate time t * , i.e., r(t * ) ≤ a for one t * ⊂ [t (n) , t (n+ ) ].

. Type III collisions are, as de ned by Wang et al. ( a), just the opposite of type II ones, i.e., r(t (n) ) ≤ a and r(t (n+ ) ) ≤ a, but d(t * ) > a for one t * ⊂ (t (n) , t (n+ ) ).

To detect collisions of type II and type III, knowledge of the particles position in between time-steps is needed. In these cases it is necessary to perform the interpolation.

. In this section, the notation di ers slightly from the one in the main content.

Here the interpolation is derived.

e two particles' spatial distance in dimension i (i = , d) in between two time-steps will be denoted as f i (t) and we seek a representation matching f i (t) = a + bt + ct + dt . e coe cients a, b, c, d can be found, by solving the linear system

f i ( ) = X i ( ) -X i ( ) = a ∶= α f ′i ( ) = V i ( ) -V i ( ) = b ∶= β f i (∆t) = X i (∆t) -X i (∆t) = a + b∆t + c∆t + d∆t ∶= γ f ′i (∆t) = V i (∆t) -V i (∆t) = b∆t + c∆t + d∆t ∶= δ
Here ∆t denotes the length of a time-step and the X i j (t) and V i j (t) ( j = , ) are the i-components of the two particles' position and velocity at time t respectively. Finally one obtains

a = α b = β c = γ -α ∆t - β + δ ∆t d = - γ -α ∆t + β + δ ∆t
Putting this together and sorting for appearances of α, γ, β, and δ one obtains . ) In more mathematical terms, these equations can be derived from the triangle inequality using the uniform (or supremum) norm. For the present purpose the given simple form should su ce. With this re-de nition a lower bound for the particle distance is given by (using Equation (A. ) again)

f (t
r min = i min f i lwr f i lwr , f i upr f i upr .
is estimate can be used in a code to decide if a more thorough check for a type II collision (see Appendix A. ) is needed. Equally the analogically de ned estimate r max can be used for type III collisions.

is means calculating the derivative of the function, nding the time(s) t * for which this equals zero and, given t * ∈ [ , ∆t], deducing the function's value at this point.

B G I P INTEGRATORS

As discussed in Section . . the particle velocity gradient tensor σ i j can diverge when the nonlinear term in Equation ( . ) becomes dominant. Physically this is related to the formation of a caustic. Based on this interpretation Falkovich & Pumir (

) developed a technique that allows to integrate σ i j despite the singularities.

ey drew inspiration for their method from earlier work by Girimaji & Pope (

). Here we present their approach as well as an alternative approach and compare them shortly.

. e temporal evolution of the gradient tensor σ is given by Equation ( . ). Neglecting the time-dependance of A i j , and limiting oneself to the onedimensional case, this equation becomes σ = Aσ τ p σ .

If now σ becomes smaller than -τ - p this equation is governed by the nonlinear term and its solution is σ(t) = (tτ p ). e local evolution of the particle density n(t) is described by the equation ṅ(t) = -n(t)σ(t), which can now be solved as well. One obtains n(t) = -n τ p (τ pt).

So the particle velocity gradient as well as the density diverge both in nite time. is, of course, is physically not possible. e maximal allowed particle density is one particle per volume of a particle, i.e., a. So the equations need to be regularized at the time t * when this value is reached n(t * ) = a =τ p n τ pt * ⇒ t * = n a + τ p . is time also de nes the maximal velocity gradient σ(t * ) = (n aτ p ) -. Assuming as a realistic value for the initial density n = η one nally obtains the estimate from Falkovich & Pumir ( ), namely σ max = η (aτ p ). If this value is reached the sign of σ has to be inverted. In physical terms this corresponds to a situation where the observed particle is overtaken by another particle.

In a three-dimensional system the scalar σ is replaced by the norm of the matrix σ . e norm ⋅ is in principle the Euclidean norm ⋅ , but in practice the Frobenius norm σ = [∑ i, j σ i j ] is su cient. When this norm reaches the maximal value η (aτ p ) the signs of all components σ i j should be ipped. is plot shows all components of the velocity gradient tensor σ i j for one particle (thin shaded lines) as well as the norm of that tensor σ (thick line) obtained from a numerical simulation (KS). In this case σ was allowed to surpass the prescribed value of η (aτ p ) (straight dashed line) by one order of magnitude. We implemented this approach and show one short example trajectory in Figure B. . In that case we allowed the norm of the particle velocity gradient tensor to be at most σ = η (aτ p ). As soon as it surpassed that value the signs of all components σ i j were inverted. is example shows nicely that, as soon as the limiting value is exceeded, σ is bound to diverge. On the other hand the ipping of the signs of the components σ i j makes sure that it returns to physical values.

. . hortcomings of the approach e sign-ipping method works quite good in simulations of the particle gradient tensor as can be seen in Figure B. . It avoids getting in nite values and therefore assures successful numerical integration. However its results are not necessarily as precise as they could be. In his section we will explain some of the shortcomings of this technique and illustrate them in a simple example.

In Appendix B. a possible way to determine a physically sound maximal value at which the signs of the particle velocity gradient tensor σ should be inverted was derived. But in reality, no matter, how good the estimate of this maximal value, one will always ip the signs "too early". Consequently all following values of σ i j will su er from errors. is is impressively seen, if one tries this technique in a simple one-dimensional case with known solution, e.g.,

x ′ (t) = + x(t) , x( ) = ⇒ x(t) = tan(t) (B. )
Now, if one tries to reproduce this result employing the sign-ipping technique, one will nd that the determination of the rst singularity's t-value is relatively exact. But the location of all subsequent singularities will vary with the choice of the maximal x at which its sign is changed, as well as with the time-step. Two cases are presented in Figure B. . ere, Equation (B. ) was integrated by a simple Euler scheme, employing two di erent techniques: e sign-ipping and a so-called GIP integrator. e latter will be discussed in Appendix B. and here we will only be concerned with the former. As in this case the exact solution is known, it is possible to determine the "optimal" maximal value. It is given by the solution, two time-steps before a singularity is reached, i.e., in our case tan(π -∆t). A er the sign of x has been changed, it is as if the numerical integration would start anew with a new initial condition. e real solution however lags behind this initial condition by ∆t. e thusly generated o set increases at every singularity resulting in a delay between the numerical integration and the exact solution. is e ect persists, of course, also for small time-steps as can be seen from . ere are quite a few integrators that were developed for matrix Riccati equations like Equation ( . ). Most of them, as for example the Möbius scheme by Schi & Shnider ( ), present a full specialized time-stepping method. But there is one interesting publication by Garrett & Li (

), that provides a way to turn every ordinary Runge-Kutta method into an integrator that reliably integrates matrix Riccati equations and even steps over singularities in the solutions. Here we will discuss the basic ideas of this method on the basis of Equation ( . ). e interested reader is referred to Garrett & Li (

) and Li & Kahan ( ) for details and mathematical proofs.

By the help of the Bernoulli substitution σ = P P -, Equation ( . ) can be transformed into a linear homogeneous di erential equation d dt P = BP, (B. )

  in a spherical coordinate system φ( w r ) Ratio of the ux into the collision sphere with RRV smaller than w r , Equation ( the replacement algorithm, Section . . rms root mean square s separation

F.

  Detail of "Studies of Water passing Obstacles and falling", drawings by Leonardo da Vinci (c. -). Coherent structures of di erent sizes are well brought out. the Richardson cascade. External forcing generates turbulence at length scale L. e large scale vortices of this size break up to give ever smaller vortices. Finally, at length scale η, the energy gets dissipated by viscosity.

  Figure . shows the energy spectrum obtained in a famous wind-tunnel experiment by Comte-Bellot & Corrsin ( ). See also the energy spectrum obtained from our own numerical simulations in Figure . . Typically a scaling according to Equation ( . ) becomes only clearly visible for very large Reynolds numbers.

  of radius a moving on a straight line with average relative velocity ⟨w⟩ (with respect to the other particles). Such a sphere can collide per unit time with all particles in a cylinder of radius a and length ⟨w⟩ around its trajectory.also note with anticipation that novel experimental work is under way byCapelo et al. ( ).

  the integral in Equation( . ). Shown is a central particle surrounded by the collision sphere of radius a, as well as streamlines of the ow. Arrows pointing radially outward, represent the out ow. Inward pointing arrows represent the radial component of the velocity eld where it enters the shell-those are added up by the integral in Equation ( . ).

  kernel as determined in a numerical simulation by Rosa et al. ( ) at Re λ ∼ . Note the strong increase at about St ∼ . . For larger Stokes numbers the collision rate decreases only slowly.

  . . austic/sling collisions and random uncorrelated motion Falkovich et al. ( ) and Wilkinson et al. ( ) noticed independently that the enhancement of the collision rate may not be solely due to the preferential concentration e ect. In fact, they predicted another mechanism, which could explain the very strong increase in the concentration rate observed at St ∼ . (Figure . ). is e ect is due to an increase of the RRV, independently of the clustering of particles. It has recently been veri ed experimentally by Bewley et al. ( ). A third independent approach to describe the same phenomenon was presented by Février et al. ( ) and Simonin et al. ( ). It has become known as "random uncorrelated motion" (Reeks et al.

  is has been done by Ducasse & Pumir ( ) in a model ow (kinematic simulations). Yet another way of describing essentially the same phenomenon has been proposed by Février et al. ( ) and Simonin et al. (

F.

  Numerical veri cation of Equation ( . ) shown here as lines. e points represent the results of a numerical simulation. σ p c (w)At this point it is instructive to perform a simple numerical experiment.To this end we rst transform Equations ( . ) and ( . ) to spherical coordinatesw ⋅ x = w sin(θ) cos(ϕ), w ⋅ ŷ = w sin(θ) sin(ϕ), w ⋅ ẑ = w cos(θ).en we integrate the polar and azimuthal contributions to obtain the PDFs of the absolute value of the relative velocity. e nal results are pWe estimated both PDFs in a numerical simulation, taking into account all particle pairs in the one case, and only colliding pairs in the other. e results, shown in Figure . , support our reasoning.

  the calculation of the contact time τ c . Two particles moving with a relative velocity w are just about to collide. e thick line passing through the point where both particles touch, gives the distance D = a cos θ for which both particles will intersect. Dividing this distance by their relative velocity gives the contact time τ c .

Finally

  introducing the dimensionless quantity ζ = a (ςτ c ) we obtain p(τ c ) that Jørgensen et al. () derive a very similar PDF for a slightly di erent situation.

F

  . e contact time PDF from Equation ( . ), shown here as a thick continuous line. e symbols indicate results from two numerical simulations with di erent size particles. Furthermore two pointed lines are shown to guide the eye along the limiting behavior for small (∼ τ c ) and large (∼ τ - c ) contact times.

  Figure . the energy spectrum from one of our DNS is shown. It corresponds well to classical experimental data (Comte-Bellot & Corrsin

  B +a Re λ = B +b Re λ = R + Re λ = F . (a) e RDF for, from top to bottom, St = . , . , . , . , . , . , . . Smaller Stokes numbers were omitted, because the slopes do not behave monotonically in the whole range (comp. Figure . b). Additionally to the numerical results (symbols), power law ts according to Equation ( . ) (lines) are shown. is law holds only true for r < η; larger values are grayed out in the plot. (b) e power law coe cient c g, as de ned in Equation ( . ) from our DNS in comparison to other published results. Shown is data from Bec et al. ( , B +a ), Bec et al. ( , B +b ), and Rosa et al. ( , R + ). which gives the correct asymptotic behavior for large distances r. Most recent studies however content themselves with the simpler functional form, which holds for r η ≪ g(r) = c g, (r η) -c g, ( . ) and so do we. Results of our DNS and ts according to Equation ( . ) are shown in Figure . a. As pointed out by Bec et al. ( ), the latter formulation of the RDF in Equation ( . ) stands in close relation to the correlation dimension known from dynamical systems theory (Grassberger & Procaccia ). is allows us to compare our results to those in Bec et al. ( , ) by the simple transformation c g, = -. e comparison is shown in Figure . b, where we also show data from Rosa et al. ( ). Our results are in good correspondence with the published data. As explained in Section . . , sole the value of the RDF at the collision radius, g( a), enters in the kinematic determination of the collision rate. In Figure . a we compare our results again to data from Rosa et al. ( ) and, in addition, from Woittiez et al. (

  e RDF at the collision radius g( a). (b) e RRV ⟨ w r ⟩ normalized by the Kolmogorov velocity u K for di erent Stokes numbers. Alongside our data (symbols), we show data from Rosa et al. ( , R + ) and Woittiez et al. ( , W + ) for comparison.

  ( ) and Woittiez et al. ( ), in Figure . b. In this case the values from Woittiez et al. ( ) lie above ours, but the correspondence with the results from Rosa et al. (

Figure . a

 a Figure . a shows the ghost collision kernel, Γ GCA , and compares it with the value of Γ . Clearly, Γ GCA overestimates the real collision kernel, the di erence being due to the multiple collisions collected in Γ m . To quantify the error introduced by the GCA, Figure . b shows the ratio Γ m Γ for a range of Stokes numbers. It shows that the GCA overestimates the collision

  . . n example of a trajectory displaying multiple collisions Aside the striking functional form of P(N c N c ≥ ), Figure . a reveals another astonishing feature. Even for particles with a relatively strong inertia, e.g., St = , we nd cases, where a same pair collides eight times.Figure .illustrates the di erence between pairs of particles that experience either

  Scatamacchia et al. ( ). But up to now, this feature was mostly documented for Lagrangian tracer particles. e trajectories in Figure.are from particles with St = . us here we report on the novel fact, that also inertial particles may exhibit this feature.

.

  d

  (τ c, N c = ) T L P(τ c, N c > ) T L P(τ c, ) (τ c, N c = ) T L P(τ c, N c > ) T L P(τ c, )F . e rst contact time PDF for St = . , as well as the same PDF conditioned on two di erent facts. One time, that pairs will collide only once, P(τ c, N c = ). e other time, that pairs collide several times, P(τ c, N c > ). Again the data is shown in semi-logarithmic (a) as well as doubly-logarithmic (b) scaling.

  PDFs of the time between contacts for the rst and second, second and third, as well as for the third and fourth contacts. e data is for St = . and shown in semi-logarithmic (a) and doubly-logarithmic (b) scaling.

  cients according to Equation ( . ) for our data.

F

  . e ratio of the mean values of the RRV in the GCA and for only rst collisions as a function of the Stokes number. In the inset the mean values are shown for only rst collisions (circles) and for all multiple collisions (triangles) in units of the Kolmogorov velocity u K .

  Wang et al. ( a) investigated such an algorithm and their results are in accordance with ours. In a companion paper Zhou et al. (

  .

  and once for half that value. As in Figure . a, the dashed line represents the expected result for Γ GCA , when St → .

  e kinematic and dynamic collision kernel for the three di erent algorithms, GCA, R , and R . (b) e RDF for the two replacement algorithms R ( lled symbols) and R (open symbols).

  Figure . b, which shows the RDF in case of the replacement algorithm. By replacing the particles, an arti cial drop in the RDF at r = a is created. erefore the RDF at the collision radius g( a), shown in Figure . a, cannot be estimated correctly. e same is true for the RRV at contact, shown in Figure . b.

  ( ) and independently by Wilkinson & Mehlig ( ) and Wilkinson et al. (

.

  Remember the decomposition proposed by, e.g., Gustavsson & Mehlig ( a) or Ducasse & Pumir ( )

  Figure . b. e work of Abrahamson ( ) has been criticized as being inexact, because it takes not into account the multiscale structure of turbulent ows. erefore alternative versions of the function h S (St, Re) have been proposed, e.g., by Völk et al. ( ). Mehlig et al. ( ) derive by dimensional analysis h S (St, Re → ∞) ∼ St . Such a behavior is not present in Figure . b, possibly due to a too small Reynolds number in our simulations.

  kernel for three di erent particle radii (tuned with the help of the density ratio ρ p ρ f ). In panel (a) it is shown in the scaling according to Γ ST and only for a range of small Stokes numbers. e value of Γ ST is also shown as a straight dotted line. Panel (b) shows the full range of Stokes numbers in the scaling imposed by Γ A .

  r ), St = . φ( w r ), St = .

  r ), St = . φ( w r ), St = .

  . Nevertheless for larger Stokes numbers the three curves, which are for particles of di erent sizes, coincide. Only for small Stokes numbers a di erence is present. is can be seen from the zoom in panel (b) which gives the same data as panel (a) but in logarithmic scaling. of the scaling of the kinematic collision kernel with the particle radius a (circles). Shown is data from Rosa et al. ( ), who report on c g, and c w , . e straight dotted lines correspond to the scalings for shear collisions (∼ a ) and sling/caustics collisions (∼ a ). e dashed line is the scaling of shear induced collisions including the preferential concentration e ect.

  Figure . . eir results con rm ours from Section . . For very small Stokes numbers the scaling of the collision kernel approaches the value expected from Sa man & Turner ( ). But for growing values of the Stokes number, the scaling exponent quickly approaches . In this sense, the results of Rosa et al. ( ) provide evidence that for St ≳ . the collision rate is determined predominantly by sling/caustics collisions. We note that this phenomenon had in fact been predicted, in the case of simple ow models by Gustavsson & Mehlig ( a, b) in the spirit of earlier work by Wilkinson et al. ().

  shows four di erent collision kernels obtained in KS. Γ GCA has been obtained in the GCA, Γ takes into account only rst collisions of a same pair, Γ m represents multiple (i.e., at least two) collisions between a same pair, and Γ R was determined with the alternative algorithm R . e straight dotted line gives the Sa man-Turner estimate for St → (Section . . ). Panel (b) shows the ratio Γ m Γ . As always in this chapter, lled symbols represent data obtained with a volume fraction Φ = . × -and empty symbols such obtained with Φ = . × -.

  a) and the RRV (b) at contact as obtained in our KS in comparison to results from DNS. e di erence between lled and empty symbols is explained in the caption of Figure . .

Figure

  Figure . b shows the ratio Γ m Γ as in Figure .b, but here for KS. Again the qualitative features are similar, but here we nd a maximal error introduced by multiple collisions which is almost %. For comparison, in DNS the maximal error was around %. We conclude that while the overall number of collisions is reduced in KS, the ratio of multiple collisions is higher than in DNS.e reasons for the large quantitative di erences between KS and DNS become a bit clearer from Figure.. ere the two constituents of the formulation of the collision kernel as proposed by Sundaram & Collins ( ) [Equation ( . )] are shown. In comparison to DNS there is almost no clustering in KS, as can be seen from the RDF shown in Figure . a. Also the RRV is much smaller (Figure . b). Apparently both e ects leading to an enhancement of the collision rate-preferential concentration and the sling/caustics e ect-are strongly reduced in KS.Further evidence is provided to some extent by Figure.. It shows the ratio of the mean collision velocity for only rst collisions and for all collisions in the GCA. e corresponding data was shown for DNS in Figure.. Because multiple collisions are so rare for larger Stokes values and because the statistics is dominated by sling/caustics collisions, this ratio drops very quickly in the case of DNS. For KS however, this ratio stays signi cant even for larger values of the Stokes number. is underlines again the suppression of sling/caustics collisions in KS.e reduced preferential concentration provides also an explanation why we nd in Figure.a that Γ R ≈ Γ .e di erence between these two collisions kernels was shown to be proportional to the volume fraction in Chapter . More precisely it depends on the e ective volume fraction in the region of a particle, namely g(r)Φ. erefore, if g(r) is smaller in KS, as was shown in Figure.a, also the di erence between Γ and Γ R will be reduced.Finally, in Figure . , we show the PDF of the collision velocity for St = . . It is qualitatively similar to the one for DNS shown in Figure . . But it

  Figure . b shows the ratio Γ m Γ as in Figure .b, but here for KS. Again the qualitative features are similar, but here we nd a maximal error introduced by multiple collisions which is almost %. For comparison, in DNS the maximal error was around %. We conclude that while the overall number of collisions is reduced in KS, the ratio of multiple collisions is higher than in DNS.e reasons for the large quantitative di erences between KS and DNS become a bit clearer from Figure.. ere the two constituents of the formulation of the collision kernel as proposed by Sundaram & Collins ( ) [Equation ( . )] are shown. In comparison to DNS there is almost no clustering in KS, as can be seen from the RDF shown in Figure . a. Also the RRV is much smaller (Figure . b). Apparently both e ects leading to an enhancement of the collision rate-preferential concentration and the sling/caustics e ect-are strongly reduced in KS.Further evidence is provided to some extent by Figure.. It shows the ratio of the mean collision velocity for only rst collisions and for all collisions in the GCA. e corresponding data was shown for DNS in Figure.. Because multiple collisions are so rare for larger Stokes values and because the statistics is dominated by sling/caustics collisions, this ratio drops very quickly in the case of DNS. For KS however, this ratio stays signi cant even for larger values of the Stokes number. is underlines again the suppression of sling/caustics collisions in KS.e reduced preferential concentration provides also an explanation why we nd in Figure.a that Γ R ≈ Γ .e di erence between these two collisions kernels was shown to be proportional to the volume fraction in Chapter . More precisely it depends on the e ective volume fraction in the region of a particle, namely g(r)Φ. erefore, if g(r) is smaller in KS, as was shown in Figure.a, also the di erence between Γ and Γ R will be reduced.Finally, in Figure . , we show the PDF of the collision velocity for St = . . It is qualitatively similar to the one for DNS shown in Figure . . But it

  Figure . b shows the ratio Γ m Γ as in Figure .b, but here for KS. Again the qualitative features are similar, but here we nd a maximal error introduced by multiple collisions which is almost %. For comparison, in DNS the maximal error was around %. We conclude that while the overall number of collisions is reduced in KS, the ratio of multiple collisions is higher than in DNS.e reasons for the large quantitative di erences between KS and DNS become a bit clearer from Figure.. ere the two constituents of the formulation of the collision kernel as proposed by Sundaram & Collins ( ) [Equation ( . )] are shown. In comparison to DNS there is almost no clustering in KS, as can be seen from the RDF shown in Figure . a. Also the RRV is much smaller (Figure . b). Apparently both e ects leading to an enhancement of the collision rate-preferential concentration and the sling/caustics e ect-are strongly reduced in KS.Further evidence is provided to some extent by Figure.. It shows the ratio of the mean collision velocity for only rst collisions and for all collisions in the GCA. e corresponding data was shown for DNS in Figure.. Because multiple collisions are so rare for larger Stokes values and because the statistics is dominated by sling/caustics collisions, this ratio drops very quickly in the case of DNS. For KS however, this ratio stays signi cant even for larger values of the Stokes number. is underlines again the suppression of sling/caustics collisions in KS.e reduced preferential concentration provides also an explanation why we nd in Figure.a that Γ R ≈ Γ .e di erence between these two collisions kernels was shown to be proportional to the volume fraction in Chapter . More precisely it depends on the e ective volume fraction in the region of a particle, namely g(r)Φ. erefore, if g(r) is smaller in KS, as was shown in Figure.a, also the di erence between Γ and Γ R will be reduced.Finally, in Figure . , we show the PDF of the collision velocity for St = . . It is qualitatively similar to the one for DNS shown in Figure . . But it

.

F

  . e probability for a pair of particles to undergo an N c -th collision conditioned on the fact that it collided at least once. e data shown here was obtained in KS.

  Figure . . An analogous graph was shown for DNS in Figure . a.

F

  . e rst contact time PDF obtained in KS for di erent Stokes numbers as indicated in the legend. e data is given in semi-logarithmic (a) and doubly-logarithmic (b) scaling. e same quantity is shown in Figure . for DNS.

  Figure . for the rst contact and for di erent values of the Stokes number. As in DNS, we nd exponential tails for long contact times (Figure . a). e dependence on the Stokes number seems to be a bit less pronounced. e most striking di erence however is the absence of the power law behavior for intermediate contact times, even at large Stokes numbers (Figure . b).

  Figure A. again, one can determine the point C, radius a before an eventual collision. eir vectorial distance is r = X -X and they move relative to each other with velocity w = V -V . More details are given in the text.

  In this plot the third order polynomial interpolation for some arbitrary α, β, γ, and δ is shown. Also shown are the di erent versions of the upper and lower bounds, given by Equations (A. ) and (A. )Hence we can estimate the bounds by calculating the extrema of each term in Equation (A. ). Doing the math² and taking into account, that even if there is no extremum within the interval ( , ∆t) the function might be extremal at the interval's borders, one nally obtains the following two values as lower and upper boundsflwr = min(α, γ) + min , β ∆t + min -(β + δ)∆t, fupr = max(α, γ) + max , β ∆t + max -(β + δ)∆t, (A. ) or another version (with di erent, but equally correct values) flwr = min(α, γ) + min , β ∆t + minδ ∆t, fupr = max(α, γ) + max , β ∆t + maxδ ∆t, (A. ) e two di erent versions stem from the two di erent representations of Equation (A. ). e above results are illustrated in Figure A. . A lower bound for the real particle distance in three dimensions can now be deduced from the above values. erefore we want to use the relation min f (t) = [min ( f (t))] , which is valid if f (t) ≥ or f (t) ≤ in the whole domain. To circumvent the resulting problems, one needs to re-de ne flwr and fupr in the following way if flwr fupr < fupr = max flwr , fupr flwr = else fupr = fupr flwr = flwr .

  . ) has been integrated numerically with a simple Euler scheme. Here we compare the signipping technique with the GIP integrator described in Appendix B. . Full lines represent the exact solution x(t) = tan(t), dotted and dashed lines give the results from signipping and the GIP scheme respectively. e straight dotted line gives the maximal value at which the sign of x was ipped [out of range in (b)].

  Figure B. b.

  Pumir A, Lévêque E ( ). "Estimating the Collision Rate of Inertial Particles in a Turbulent Flow: Limitations of the 'Ghost Collision' Approximation." In: J. Phys.: Conf. Ser. , p. Voßkuhle M, Lévêque E, Wilkinson M, Pumir A ( a). "Multiple collisions in turbulent ows." In: ArXiv e-prints (Oct.). Accepted for publication in Phys. Rev. E. arXiv:.[physics. u-dyn].he tone ge did not come to an end because we had a lack of stones [. . .] Numerical determination of the main variables . Particles with Gaussian velocity distribution . . Veri cation of collision detection scheme .

			xiv				xv
	F		.		. . ird order method cumulative probability density function Trajectory interpolation
	Ahmed Zaki Yamani ( Saudi Arabian Minister of Oil Contact times . . . . Probability of relative velocity and statistics condi-) -e third order interpolation and its bounds CONTENTS Acronyms xvi Notation xvi . g(r) radial distribution function, Equation ( . ) δ(⋅) Dirac delta function ACRONYMS G gravitation vector ∆t time-step in numerical simulations . Physically inspired approach . . Shortcomings of the approach . integrators DNS Direct numerical simulation rate of caustic formation, Section . . ε energy dissipation rate per unit mass, Equation ( . ) h S (St, Re λ ) J(a) term representing growth due to condensation in ζ dimensionless quantity in Equation ( . ) tioned on collisions . . Literature cited GCA Ghost collision approximation (Section . . ) Equation ( . ) η Kolmogorov length, Section . . Probability of contact times k wavenumber θ polar angle in a spherical coordinate system
	. b). "Prevalence of : t coe cient in Section . . A short introduction to our integral length scale . Θ(⋅) Heaviside step function / -/ / / . Voßkuhle M, Pumir A, Lévêque E, Wilkinson M ( i L . . Particles in turbulence and collisions . Outline . . Particle trajectories . n number density (of particles) κ Description of the di erent runs . N number of subscript λ Taylor scale, Equation ( . ) Veri cation of results the sling e ect for enhancing collision rates in turbulent suspensions." In: ArXiv e-prints (July ). arXiv: . . . Radial distribution function µ dynamic viscosity [physics. u-dyn]. . Navier-Stokes equations . . e vorticity . . . . Radial relative velocity . . ν kinematic viscosity Dynamic and Kinematic Collision kernel ξ t coe cient in Section . . Self similarity and the Reynolds number . Turbulence ϖ ( ) distribution function, Equation ( . )
					. . Π(k, t)	e cascade picture energy ux, Section . .
					. .	Homogenous and isotropic turbulence
					. .	Energy ux and energy dissipation rate
					. .	Kolmogorov scales
					. .	Kolmogorov spectrum
			. Motion of spheres in a turbulent ow
					. .		e Stokes number
					. .		e particle velocity gradient tensor
			.		Droplets in turbulent clouds-An overview
	St . t T L u	≤	t .	,	h a n Comparison to published results k W u b i x = s i n ( t ) c o s ( t ) l n t e droplet size distribution , Stokes number, Equation ( . ) t e collision kernel time Typical values Analytical results on collision kernels c o s ( t ) . . . . . . . . Collision kernel for an ideal gas a particle radius integral time scale Cumulative . A, A i j velocity gradient tensor uid velocity eld	t	I
	. . . . B, B i j Zero inertia . Conclusion and Perspectives matrix de ned in Equation (B. ) particle velocity vector In nite inertia . . Preferential concentration . . Caustic/sling collisions and random uncorrelated motion . . Ghost collision approximation . . Collision e ciency . . Two interesting remarks . Description of the approach . . Integration of particle trajectories . Collision kernels c g parameter in Equation ( . ) c w volume relative velocity between two particles t parameter similar to c g but for ⟨ w r ⟩ C constant coordinate, typically as argument in a eld . Multiple collisions . . Co Courant number, Section . particle position vector Contact time d dimensionality . Detecting sling/caustics collisions correlation dimension, Section . . Greek letters V w x X . Conclusion and perspectives ii . Collision detection scheme E(k, t) energy spectrum α t coe cient in Equation ( . ) iii β t coe cient in Equation ( . )
					. .	Cell linked-list algorithm . Linear extrapolation
					xvi			xi xiii

. Ghost collision and rst collision kernel . Multiple collision . . An example of a trajectory displaying multiple collisions . Contact time statistics . . Dependence on the Stokes number . Relation to sling/caustics collisions . Alternative algorithms . . Replacement algorithm . . Dynamic collision kernel . . Kinematic collision kernel . Conclusion and Perspectives / . Di erent scaling of the collision kernels .

GIP Generalized inverse property (Appendix B) KS Kinematic simulation (Chapter ) PDF Probability density function R Replacement approach replacing one particle (Section . . ) R Replacement approach replacing both particles (Section . . ) RDF Radial distribution function (Section . . ) RNG Random number generator RRV Radial relative velocity (Section . . ) NOTATION It has been attempted to keep the meaning of symbols consistent throughout the whole document. Here is a list of the most important de nitions used in this text. c collision e ciency, Section . . f (a) droplet size distribution, Section . . c overall rate of collision p probability density function (PDF) P probability or probability density function (PDF) pressure eld divided by uid density ρ, Section . Q forcing eld in the Navier-Stokes equation ( . ) r distance vector R i j velocity autocorrelation function, Equation ( . ) Ri j Fourier transform of the velocity autocorrelation tensor R i j rate of collision encountered by one particle Re Reynolds number S(k, t) transfer of energy into mode k, Section . . γ(St, Re λ ) factor relating V rms and u K , Section . .

xviii Γ collision kernel, Section . . ρ mass density σ, σ i j particle velocity gradient tensor, Section . .

  life of a raindrop, from the moment it is formed by condensation to the time it has grown large enough to fall on Earth can be simpli ed as follows. First, it forms by condensation of vapor on tiny aerosol particles.ose are called cloud condensation nuclei (CCNs) and have sizes < µm(Devenish et al.

		. .	he droplet size distribution
		e
	A non-exhaustive list includes the papers of Beard & Ochs	; Blyth	; Devenish et al.
		; Grabowski & Wang	; Jonas	; Pinsky & Khain	; Shaw	; Vaillancourt &
	Yau	.

  on the type of cloud, its age, the region within the cloud, and other factors. In general, the energy dissipation rate ε has been documented to vary in between -m s and at least -m s (see Grabowski & Wang ; Pinsky & Khain , and references therein). Making use of the relation λ u rms

See for example the works of Sundaram & Collins ( ), Wang et al. ( ), and Williams & Crane ( ). But we also note the studies of Rosa et al. ( ) and Woittiez et al. ( ), which suggest, that the combined e ects of gravity and turbulence cannot be understood independently. In fact, Siebert et al. ( b) use a di erent, but conceptually similar approach. depending

  Wilkinson et al. ( ) remind that according to the approach of Maxey ( ) the clustering should be most pronounced for St ∼ . Numerical results, like those shown in Figure . , however show that the collision rate starts to grow explosively for St ≪ . Furthermore, the collision rate does not abruptly fall for St > , but slowly approaches Abrahamson's ( ) prediction for St → ∞. Wilkinson et al. ( ) argue further, that the usual clustering process could only bring together particles on a distance similar to the smallest scales of the ow. In a typical cloud however-see Section . . -there is on average only one droplet per cube of edge size η.

erefore Wilkinson et al. (

  Also the preferential concentration e ect is a ected by the GCA, as Reade & Collins ( ) show. ey compare the RDF obtained from simulations with ghost particles and from simulations, where particles bounce o each other like hard spheres. e two values for g(r) di er with a dependence on the Stokes number St and the particle size a. e inaccuracy of the GCA is by no means restricted to a particular numerical algorithm. In fact, Andersson et al. ( ) and Gustavsson et al. ( ) point out, that all theories based on similar reasoning as introduced by Sa man & Turner (

When the particles separate and collide again a erward, either with each other or with further particles, these contacts are counted as new collisions. Evidently, this approximation leads to an overestimation of the collision rate-this has been noted by Zhou et al. (

).

  Wang et al. ( b) show that actually Sa man & Turner's () second derivation di ers from the rst, in that it makes use of the "cylindrical formulation". e correct approach is however, as Wang et al. ( b) go on to demonstrate, the "spherical formulation".Another remark with practical importance is made by Wang et al. (

  Equation ( . ) provides a complete representation of the velocity eld at any point in space. is is a big advantage over DNS which rely on interpolation schemes to calculate the velocity in between grid points. To obtain the particle trajectories, we integrate Equations ( . ) and ( . ) with a low storage Runge-Kutta scheme of fourth order. is scheme is presented in detail by Gottlieb et al. ( ) based on work by Carpenter & Kennedy (

	. .	ntegration of particle trajectories

Our simulations are based on those described by Ducasse ( ) and Ducasse & Pumir (

). Further details may be found there.

  It may be useful to know upper and lower bounds of f (t) in the interval [ , ∆t], so one can decide whether it is necessary to proceed with the collision checking. To this end one can make use of the following facts,¹ which hold for arbitrary functions g k (t)

		) = α -+ β t -t ∆t = α + (γ -α)	+	t ∆t ∆t t + t ∆t -∆t + γ t t ∆t	t ∆t + δ -+ β t -∆t t -t ∆t ∆t t + t ∆t	+ (β + δ) -	t ∆t	+	t ∆t (A. )
	and	min	k	g k (t) ≥	k	min (g k (t))
		max	k	g k (t) ≤

k max (g k (t)) .

(

Γτ K ( a) St Γ GCA (Φ ) Γ (Φ ) Γ m (Φ ) Γ R (Φ )

ACKNOWLEDGMENTS I thank my supervisor Alain Pumir. I thank my co-supervisor Emmanuel Lévêque and our collaborator Michael Wilkinson. I thank the committee members Jérémie Bec, Michel Lance, Bernhard Mehlig, and Mike Reeks. I thank the Centre Blaise Pascal, the Pôle Scienti que de Modélisation Numérique, and the physics lab at the ENS Lyon. I thank all developers of FLOSS.

I thank all colleagues from the ENS.

I thank all former colleagues and mentors.

I thank all friends I have made in Lyon.

I thank my parents, family, friends and everyone who supported me.

respectively, were introduced. Here I designates the × identity matrix. Simply integrating Equation (B. ) may lead to numerical instabilities. But one has the freedom to scale the matrix P by any non-singular × matrix M in order to prevent these problems. Because P P -= (P M)(P M) - such a scaling does not in uence the results. In fact one may even rescale P at every time-step. Garrett & Li (

) explain that M = P -is an especially convenient choice.

ese observations lead to the following recipe: Integrate Equation (B. ) with an ordinary Runge-Kutta scheme and regularly rescale P such that P = I . Whenever necessary extract σ by calculating σ = P P -.

Garrett & Li ( ) show that with this recipe one can simply step over all singularities without any further knowledge about the details. As can be seen in Figure B. it works correctly and especially it does not display the problem of shi ed singularities, as discussed in Appendix B. . We have implemented this technique in our code and the results shown in Chapter were obtained with it. We also compared the two approaches-the GIP integrator and the sign-ipping-in our numerical simulations of particles entrained in a model ow. But we could not nd any obvious di erences.

Probably the problem of the shi ed singularities regularizes itself in a more complex environment. Nevertheless, the GIP integrator rests a convenient technique for integrating Equation ( . ).

LITERATURE CITED

is is the complete list of cited articles, books, and other resources. For the comfort of the reader I tried to include a digital object identi er ( ) where possible. is should allow a stable (in time) access to articles published online. For other online resources, which do not have a , like books that are freely available on the internet, I included a simple

. ese are however, unfortunately not guaranteed to persist for long times. Actually I noticed some already disappearing over the process of writing this manuscript. A good and usually successful practice in these situations, is to head to a web archiving service like the Wayback Machine (http://archive.org).