Unifying ideas on mixing and atomization

The mechanisms building the overall concentration distribution in a scalar mixture, and the drops in a spray, are examined successively. In both cases, the distributions belong to a unique family of distributions stable by self-convolution, the signature of the aggregation process from which they originate.


Scope
A dye diffusing in a diluting medium while it is stirred and a set of dispersed droplets in the spray entrained by the wind blowing over a liquid surface are at first sight very dissimilar events. However, discovering the process by which they are built reveals unexpected analogies. As explained below, both the concentration levels in a stirred mixture and the liquid drops in a spray result from a process of random addition which has its counterpart on the shape of the concentration distribution in the mixture, and the drop size distribution in the spray. These distributions are of the form where X denotes the dye concentration normalized by its mean C/ C or the drops' size normalized by the average size d/ d . The parameter n depends on the operating conditions.

Mixing
A mixture is a transient state between the initial segregation of the constituents, and their ultimate homogeneity. The overall mixing process of a drop of dyed fluid in a stirred medium involves two phenomena: a process of dispersion of the drop in the diluting medium by which the phases interpenetrate, and a process of interaction between the dispersed elements from which homogeneity arises. We report here on findings suggesting that the nature of the interaction is of an aggregation type and that this phenomenon is the key to understanding the concentration distribution of a stirred mixture. Let a jet of water plus a diluted fluorescent dye (disodium fluorecein) be discharged in a square, transparent, long duct. For a given duct cross section, the injection diameter d and the velocity of the co-flow at the entrance of the duct can be varied so that the average concentration of the dye in the channel C can be set at will. Since the cross-section of the duct and the average velocity of the mixture in the downstream direction are constant, the average concentration is conserved. The jet velocity u is such that the Reynolds number Re = ud/ν 10 4 .
The flow is made visual by means of a plane argon laser sheet slicing the duct along its axis and, as can be seen on figure 1, the dye rapidly invades the whole duct cross section, erasing its concentration differences as it travels downstream to relax towards a more or less uniform mixture.
After the dye has filled the channel cross section and evolves around a constant average concentration, the distribution P(C) presents a skewed bell shape, which gets narrower around C in time (axial distances are converted to time through the average axial velocity with confidence as the radial velocity profile in a turbulent duct is flat [1]). The shape of P(C) is very well described by a family of one-parameter distributions, namely gamma distributions The parameter n is adjusted at each downstream location for the gamma distribution of (2) to fit the experimental one. It is seen in figure 2 that the fairness of the fit holds for the whole concentration range, down to quite low probability levels, and accounts for the downstream deformation of P(C) through the single parameter n, whose dependence on the downstream location is quite strong: figure 2 suggests a power-law dependence with an exponent close to 5/2. The dependence of n on the jet Reynolds number is, although noticeable, very weak.

Stretching enhanced diffusion
The stirring motions progressively convert a compact blob in a set of sheets of increasing surface and decreasing thickness [2]- [5]. The intersects of these sheets with the visualization plane are visible in figure 1 in the form of ligaments. Let s(t) be the distance between two material particles in the direction z perpendicular to a sheet (figure 3), and σ(t) = ∂ ln s(t)/∂t its rate of compression. If c(z, t) is the scalar concentration profile across the sheet, the convection-diffusion transport equation reduces to a one-dimensional problem provided the radius of curvature of the sheet is large compared to its thickness [6]. In that case, the direction z aligns with the direction of  maximal compression and for a species with diffusion coefficient D By the change of variables τ = D t 0 dt /s(t ) 2 and ξ = z/s(t), (3) reduces [7]- [10] to a simple diffusion equation ∂c(ξ, τ)/∂τ = ∂ 2 c(ξ, τ)/∂ξ 2 . The topology of the stirring motions select particular forms for s(t). Starting with a sheet of initial uniform concentration and thickness s 0 , the maximal concentration for z = 0 can be written as We describe several generic examples: in incompressible flows in two dimensions where the length of material lines grow like γt [11], the mean transverse thickness of the scalar filaments decreases as s(t) = s 0 / 1 + (γt) 2 and thus τ = (Dt/s 2 0 )(1 + ((γt) 2 /3)), providing c(0, t) ∼ (t/t s ) −3/2 for t > t s , with t s ∼ (1/γ)Pe 1/3 , where Pe = γs 2 0 /D is a Péclet number. If material surfaces in three dimensions grow like (γt) 2 , then [12], s(t) = s 0 /(1 + (γt) 2 ) and τ = (Dt/s 2 0 )(1 + 2 3 (γt) 2 + 1 5 (γt) 4 ), providing c(0, t) ∼ (t/t s ) −5/2 for t > t s , with t s ∼ 1 γ Pe 1/5 . For flows in which the length of material lines increases exponentially with time like e γt as realized by a succession of stretching and folding motions in random flows [13], s(t) = s 0 e −γt and τ = (Dt/s 2 0 )(e 2γt − 1) providing c(0, t) ∼ e −γt for t > t s with t s ∼ 1/2γ ln Pe. These timescales are the relevant mixing times as soon as the inverse of the elongation rate γ −1 is smaller than the diffusive time of the sheet constructed on its initial size s 2 0 /D, that is for Pe 1. In this limit, t s is essentially given by the time needed to deform the sheet γ −1 and pure diffusion (for which c(0, t) ∼ (Dt/s 2 0 ) −1/2 ) is enhanced by the substrate motion.

Composition rule and kinetic equation
However, the sheets interact as they move in the flow so that their diffusive boundaries interpenetrate to give rise to new sheets whose concentration profile is the addition of the original ones. This elementary interaction rule is a consequence of the linearity of the Fourier diffusion (2) and dictates the kinetic evolution equation for P(C, t).

Linearity.
The concentration profile c(ξ, τ) of a substance evolving by diffusion is obtained at later times from its initial spatial profile c(ξ, 0) by the summation of the impulse response of the diffusion equation (c(ξ, 0)dξ/2 √ πτ)e −ξ 2 /4τ over space The concentration profile resulting from the merging of two lamellae 1 and 2, each having a profile c 1 (ξ, τ) and c 2 (ξ, τ) is thus the addition of the original profiles If Q(c 1 ) is the probability to find a sheet with concentration c 1 in the medium and Q(c 2 ) the probability to find c 2 , then the probability to find c after the sheets have merged is The above composition operation defines a convolution of the original distribution.

Kinetic equation.
Let Q(c, t) be concentration distribution of the elementary sheets and Q(s, t) = ∞ 0 Q(c, t)e −sc dc its Laplace transform. The change of Q(c, t) between t and t + t is composed of two ingredients: firstly, the decay of c due to the stretching of the sheets results in a global shift of Q(c, t) towards the low concentration levels without altering its shape as (−∂/∂c)( dc/dt Q) t = −γ(∂/∂c)(cQ) t. Secondly, the sheets interact with nearby neighbours and because of the irregular stirring motions, the addition of the concentration levels is made at random among those available in the population Q(c, t) at time t which therefore evolves, through this process, by self-convolution as ⊗2 . The time it takes to complete a convolution is the time needed to coalesce two sheets t ≈ γ −1 = −c(0, t)/(dc(0, t)/dt); thus the continuous version of this convolution step can be written, in the Laplace space, asQ(s, Taking the limit t → 0, the global rate of change ofQ(s, t) accounting for the two ingredients is thus given by The above interaction rule leads, irrespective of the initial condition for Q(c, 0), to decaying exponential distributions [14] of mean c = c(0, t) as defined in (4). The concentration C results from clusters of several of these interacting sheets. If these clusters are made of n independent coalesced sheets, then the distribution of the concentration P(C, t) is the nth convolution of Q(c, t), i.e.P(s, t) =Q(s, t) n whose overall kinetic equation derives from (8) as With γn = dn/dt, the asymptotic solution of (9) isP = (1 + C (s/n)) −n , providing a gamma distribution of order n for P(C, t), which is a convolution of n exponentials. The average concentration C of the mixture is conserved when the damping factor γ balances exactly the coalescence rate, when n = 1/c(0, t). This condition is realized provided n ∼ (t/t s ) 5/2 in this three-dimensional flow, as expected from figure 2. 2

Change of dimension
By changing the dimensionality of the flow from three to two, it is expected that as long as the stirring motions are sufficiently irregular to ensure an addition of the concentration levels at random, the distribution is, as in the three-dimensional case, generated by self-convolution. It is, by contrast, also expected that the dependence of n on time will be different. We illustrate this difference by a simple experiment of stirring a blob of dye with a rod in a thin layer of a very viscous fluid, by a two-dimensional quasi-periodic protocol which mimics the motion of a straw in a milk shake. The diluting fluid is pure glycerol and the drop is made of the same fluid coloured with indian ink. A number of parallel cuts is made in one direction, and then the same number at right angle, this operation defining one cycle (figure 4).
In this low Reynolds number flow (the typical Reynolds number of the motion of the rod is Re = us 0 /ν = 10 −1 ), the fluid is deformed by the passage of the rod on a scale which is given by its own size s 0 . The length of material lines is equal to the distance travelled by the rod in the medium, and it is actually observed that the net contour length of the deformed scalar drop increases in proportion to the number of cycles (figure 5). Since L/L 0 = 1 + γt ≈ 1 + 17.5 × p, the mixing time, or, alternatively, the number of cycles p needed to start mixing the drop in the surrounding medium is thus given by γt s ≈ 17.5p s ∼ (γs 2 0 /D) 1/3 = (Re Sc) 1/3 with γ ≈ u/s 0 the typical elongation rate, u ≈ 4.6 cm s −1 the rod velocity, giving, with s 0 ≈ 2 mm and Sc ≈ 10 6 , a mixing cycle p s of the order of 2.5 (see figure 5).
The maximal rate of stretch is obtained for fluid particles close to the rod trajectory, while the protocol leaves nearly unstretched fluid parcels, which therefore keep a concentration close to the initial concentration as the number of cycles is increased. Concomitantly, fluid particles are brought together in the wake of the rod, and they coalesce. The amplitude of the slicing movements is constant through the cycles, so that the average concentration of the dye must be conserved; these are the ingredients required for the aggregation scenario we have described to occur. The distributions displayed in figure 6 are actually reasonably well described by the gamma functions family of (2) with an average concentration C constant. The maximal concentration in each single scalar filament decreases as c(0, t) ∼ (t/t s ) −3/2 and the orders of the gamma distributions consistently follow n ∼ (number of stirring cycles) 3/2 , in this two-dimensional flow. As soon as a mixture is stirred in a more or less irregular manner, the addition is made at random among the levels available in the current distribution. This operation implies that the concentration distribution evolves by self-convolution and, when it actually does, gives an a posteriori precise definition of what 'random' means. This process, which basically amounts to a convolution of exponentials [15], [16]- [19], might certainly be generic of the construction mechanism of the skewed distributions with exponential tails widespread in various instances including turbulent convection [20], grid turbulence [21,22], shear layers [12], randomly stirred two-dimensional flows [23]- [25] or jets [26]. This scenario not only accounts for the distributions tail [27]- [30], but also for their whole shape and evolution.

Atomization
The disintegration and dispersion of a liquid volume by a gas stream is a phenomenon which embraces many natural and industrial operations. The entrainment of spume droplets by the wind  [31,32]. In order to compute the rate of exchanges of solutes between the ocean and the atmosphere, or to estimate the size of a combustion chamber, it is frequently desirable to have a precise knowledge of the liquid dispersion structure, in particular its distribution of droplet sizes as a function of the external parameters.

Exponential tails
The broad size statistics is a salient and fundamental feature of natural sprays formed in an uncontrolled way. Spume droplets [33,34], atmospheric aerosols [35], rain drops [36], volcanic Tephra, fuel sprays [37,38] all display a wide, highly skewed distribution of sizes, the most probable droplet sizes being close to the smallest ones and the probability of finding a drop size twice or three times larger than the mean being not vanishingly small. A generic character of these distributions is an exponential tail at large sizes. For instance, Simmons [37,38] notes that, for a large collection of industrial sprays, the distribution of sizes p(d ) is universal in shape and that its tail is well fitted by an exponential fall-off. The existing models invoked for this fragmentation process essentially rely on cascade ideas [39], following the early suggestion of Kolmogorov [40], leading to log-normal statistics of the fragment sizes. Notable exceptions in this context are the work of Longuet-Higgins [41], which shows how a simple geometrical model of ligament random break-up produces broad, skewed size distributions without resorting to sequential cascade arguments, and that of Cohen [42], which shows how pure combinatory and thermodynamic arguments [43] lead to a Poisson distribution for the fragment volumes.

Experiments
The flow configuration designed to address this problem consists of a round water jet surrounded by a coaxial air flow, an axisymmetric geometry convenient for visualization purposes. Each stream is potential at the exit of the injector, with a very weak residual turbulence. The velocity of the water stream is of the order of 1 m s −1 and that of the air stream is varied up to u = 50 m s −1 .
The spray was analysed on frozen images using a short time illumination system, and high speed video. As suggested by figure 7, at the root of the disintegration process is a shear between the light, fast stream and the slow, dense liquid. Provided the Weber number (We = ρu 2 δ/σ) and the Reynolds number (Re = uδ/ν where ρ, ν denote the density and viscosity of the gas and σ the surface tension of the liquid) are large enough [44,45], this shear induces an instability of a Kelvin-Helmoltz type forming axisymmetric waves on the liquid jet interface. Their wavelength and growth rate are controlled by the thickness δ of the velocity profile in the gas stream [44]- [47]. When the amplitude of these primary undulations is large enough, they undergo a transverse destabilization of a Rayleigh-Taylor type [48,49] caused by the accelerations ξ d Figure 8. Double flash exposure of a ligament just before and after breakup. The duration of a flash is 5 µs, and the interval between the two flashes is 1.6 ms. The resulting droplets are highlighted by a white rim.
perpendicular to the liquid-gas interface imposed by the passage of the waves. The resulting modulation of the wave crests is further amplified by the air stream forming ligaments, which ultimately break by capillary instability [50]- [52]. Provided a ligament is stretched in the wind at a rate which overcomes the rate of the capillary instability σ/ρd 3 0 based on its initial volume V = d 3 0 , where d 0 depends on δ and We [45], this volume V remains constant as it deforms and finally detaches from the liquid bulk with thickness ξ, a function of the operating conditions (air velocity, liquid surface tension, gas/liquid density ratio). This overall process is usually referred to as 'stripping' [53]. Since the ligaments are elongated at breakup (see figures 8 and 9), their typical transverse thickness ξ, which is all the more thin that the air velocity is fast, is smaller than d 0 .
A key observation is that, although the ligaments are very thin when they detach from the liquid bulk, they give rise to drops whose size is substantially larger (figure 9), and which scales like d 0 . In other words, no matter how thin a ligament is at breakup, it will form drops of the order of the size which sets the ligament volume.
The situation considered here is fundamentally different from the classical Rayleigh breakup of a liquid thread of uniform thickness in a quiescent environment: the reorganization of the liquid volume in the ligament while it stretches is a superposition of remnant motions from the liquid bulk, motions due to the transient growth and damping of capillary waves [54], motions induced by the deformation of the ligament due to perturbations in the gas stream, etc. These are so complex that they are out of reach of a microscopic analysis. However, capillary forces ultimately fragment the ligament in several blobs. When two liquid blobs of different sizes d 1 and d 2 (with, say, d 1 < d 2 ) are connected to each other, they aggregate due to the Laplace pressures difference ∝ σ(1/d 1 − 1/d 2 ). The time it takes for the coalescence to be completed is of order ρd 3 1 /σ which is also the time it takes for the neck connecting the two blobs to destabilize and break [51,52]; the confusion of these two timescales induces a nice 'coalescence cascade' [55]. For the very same reason, the blobs constitutive of the ligament tend, as they detach, to coalesce, thereby forming bigger and bigger blobs along the ligament. This is why the final mean drop size is larger than the average thickness of the ligament ξ just after it has been released from the liquid bulk. As long as the ligament is attached to the liquid bulk and is stretched in the wind, the capillary instability of its core is strongly damped. This is true for all modes whose instability rate is smaller than the stretching rate [56]. The time it takes for the ligament to detach from the liquid bulk is given by the capillary time based on its initial size T = ρd 3 0 /σ. As soon as it has detached and is no more (or much less) stretched, the capillary breakup and coalescence period develop on a shorter timescale of the order of the capillary time based on its thickness ξ.

Size distribution
Our interest is to understand the statistics of the drop sizes in the spray and our observation is that this distribution is determined by the large excursion wing of the distribution associated to d(z,t) Figure 10. An isolated ligament just before breakup covered with blobs of various sizes d matching its local thickness. The evolution of the size distribution n(d, t) is governed by (13). a single ligament breakup ( figure 14). The ligaments' dynamics is thus crucial to understand the overall spray size distribution.

Linearity.
The dynamics of the rearrangement of the fluid particles along the ligament is well described by the small slope approximation conservation laws [57] linearized around the initial state {u, d} = {0, ξ} where u denotes the axial velocity along the ligament, d its diameter (figure 10), and ξ its initial (taken as uniform) value. Obviously, this system leads to a linear evolution equation for d which describes the dynamics of the (possibly unstable) capillary waves along the ligament. Note that this solution describes most of the breakup period, the time lapse of the-fundamentally d d' Figure 11. Sketch of the blob interaction scheme in the ligament.
nonlinear-pinching period close to the droplet separation being of the order of 10 −3 the overall breakup period (∼ ρξ 3 /σ, see e.g. [52,58]). Most of the time in this process is spent moving the fluid particles apart around the initial ligament shape. The phenomena in the vicinity of the separation of the fluid particles are comparatively much more rapid.
For the same reason discussed in section 2.2, the amplitudes of the waves building the corrugations of the ligament are thus likely to obey a convolution rule.

Kinetic equation.
We model the dynamics of the set of travelling waves which overlap at random along the ligament by dividing it into a set of interacting nearby blobs. Let n(d, t) dd be the number of blobs constitutive of the ligament whose size is within d and d + dd at time t during the interaction period. The total number of blobs constitutive of the ligament at time t is N(t) = n(d, t) dd, its length L(t) = dn(d, t) dd and its volume V = d 3 n(d, t) dd (see e.g. figures 10 and 11). Let the random motions in the ligament result in ν independent layers and let q(d , t) be the distribution of the sub-blobs' sizes d in each layer with q(d , t) dd = 1. The layers are adjacent to each other across the ligament section so that ν d = d with d = dn(d, t) dd/N(t). In the course of the coalescence-like process between adjacent interacting fluid particles described above, the sub-blobs overlap and the distribution of sizes in each layer an instant of time later q(d , t + t) will result from the interaction of blobs of various sizes in the current distribution q(d , t) at time t. The average thickness of the layers d /ν is the typical mean free path of the agitation motions in the ligament. If we conjecture that the interaction is made at random with no correlation between the sizes interacting, the evolution of q(d , t) is then directed by a convolution process [14] q Since the layers are assumed independent, the distribution of the blobs' size d itself is thus n(d, t) = N(t)q(d , t) ⊗ν . The corresponding evolution equation for n(d, t) is with γ = 1 + 1/ν. Time t is counted from the moment when the ligament detaches from the liquid bulk (t = 0), and is made nondimensional by T ξ = ρξ 3 /σ, the capillary time based on the initial average blob size ξ = d 0 where d 0 = dn(d, 0) dd/N(0) (see (11)). The structure of (13) and prefactors are such that the net volume of the ligament V is conserved. The interaction parameter γ is determined by the compatibility of (13) with the initial distribution of the blobs along the ligament by γ = d 2 0 / d 2 0 with d 2 0 = d 2 n(d, 0) dd/N(0). A uniform thread of constant thickness (made of many, very thin independent layers) has γ = 1 and a corrugated ligament is such that γ > 1.
The asymptotic solution of (13) for p B = n(d, t)/N(t) is a gamma distribution of order ν = 1/(γ − 1). The very mechanism goes back to the discovery by von Smoluchowski that systems like (12) evolving by self-convolution generate exponential distributions [16], as already noted in section 1. The distribution solution of (13) is indeed a convolution of ν exponentials, providing [15] where d = dn(d, t) dd/N(t) is the current average blob diameter. These gamma shapes closely fit the experimental distributions of the blob sizes before breakup ( figure 12) and the drop sizes after ligament breakup with an average diameter d 0.4d 0 independent of the air velocity (figures 14(a) and (b)). The order ν increases slightly with the air velocity ( figure 14(c)). These facts indicate that the rearrangements and coalescence between the blobs tend to restore the average diameter d from ξ to d 0 , or a fraction of d 0 . This is made at the expense of a reduction of the number of blobs N(t) which, according to (13)    For large ν, that is for smooth and uniform ligaments giving rize to a narrow size distribution centred around ξ, one has ln( d /ξ) ∼ 1/ν. For small ν, that is for corrugated ligaments inducing a broad size distribution, one has The above anticipated trend is not incompatible with figure 14(c); this process interestingly suggests that thinner, but still corrugated ligaments formed by faster winds, or when the capillary breakup is slowed down by an increased liquid viscosity [52], produce drops with a narrower distribution (the standard deviation of the gamma distribution is ∼1/ √ ν), within logarithmic corrections.
For the given operating conditions, the diameter d 0 is itself distributed among the population of ligaments, although this distribution p L (d 0 ) is narrower than p B (d/d 0 ). The size distribution in the spray p(d ) is thus a mixture of the distribution of ligament size p L (d 0 ) and of the universal distribution of sizes after the ligament break- This composition operation stretches the large excursion wing of p B (d/d 0 ) over nearly the whole range of sizes d, and the size distributions in the spray shown in figure 14(d) thus coincide with the exponential fall-off The prefactors of the exponential slopes are about 3.5, like the order ν of the ligament gamma distributions (figure 14(c)), and increase slowly with the air velocity, as does the ratio d /ξ ( figure 14(c) and equation 15). The exponential shape of the global distribution, and the value of their argument have thus to be understood as the large-sizes behaviour, and the order ν of the gamma distributions coming from the ligament-breakup process, respectively. That step thus appears as the crucial step, building-up the broad statistics in the spray. The slight increase of the exponential slopes with air velocity (insert) reflects the variation of the gamma orders ν on d /ξ.

Conclusion and final remarks
Phenomena exhibiting the same statistics are likely to be ruled by the same microscopic mechanisms. We have purposely put in relation the mixing of a dye in a stirred medium, and the atomization of a liquid in a dispersed spray. The evolution of the concentration levels in a scalar mixture follow the Fourier equations, which are linear. Drops in a spray are formed via the destabilization of elongated objects, ligaments, whose dynamics results from the overlapping of capillary waves, mostly described by their linear version around the initial shape of the ligaments. From this fundamental property of linearity, common to our two-apparently very different-physical examples, we have derived a general evolution rule for the distributions of the quantities subjected to an addition rule, i.e. concentration levels in a mixture, and drop sizes in a spray. These distributions are stable with respect to the process giving birth to them, that is the self-convolution construction mechanism; hence the gamma distributions. Let us end with some specific remarks: it is instructive to note that a scalar mixing experiment performed at a very low Reynolds number produces composition fields very similar to those obtained at much larger Reynolds number (figures 2 and 6). The reason is that the evolution mechanism is the same. The motion of the rod in a two-dimensional viscous fluid plays the role of the random stirring motions present in high-Reynolds-number flows. Their role is to ensure the independence-in the statistical sense of (9)-of the addition of the concentration levels; a particularly simple paradigm for the impact of turbulence on mixing.
The fragmentation mechanism we have described which, somewhat surprisingly, consists of a coalescence process, is representative of situations where drops 'go with the wind', as in spume or airblast sprays. It is, in fact, generic of all situations where drops come from the capillary destabilization of a strongly corrugated ligament, which was our basic ingredient. In a remote context, the disintegration of big nuclei has been suggested to obey a similar scenario in the celebrated 'drop model' for nuclear fission [59]. Since the formation of smaller drops from a liquid volume imply that it shapes in ligaments for capillary forces to produce breakup [50]- [52], this mechanism might therefore be relevant to sprays formed by splashes or mutual droplet collisions [60]. Those occur for instance among drops with different terminal velocities in rain for which the gamma distribution has been identified long ago as a convenient empirical fitting distribution of the drops' sizes [35,36].