Lagrangian tracking of colliding droplets
aa r X i v : . [ phy s i c s . f l u - dyn ] J un Experiments in Fluids manuscript No. (will be inserted by the editor)
Lagrangian tracking of colliding droplets
RV Kearney · GP Bewley
Received: date / Accepted: date
Abstract
We introduce a new Lagrangian particle tracking algorithm thattracks particles in three dimensions to separations between trajectories ap-proaching contact. The algorithm also detects low Weber number binary col-lisions that result in coalescence as well as droplet break-up. Particles areidentified in two-dimensional high-resolution digital images by finding sets ofcircles to describe the edge of each body. This allows identification of particlesthat overlap in projection by over 80% even for noisy images and without in-voking additional temporal data. The algorithm builds trajectories from three-dimensional particle coordinates by minimizing a penalty function that is aweighted sum of deviations from the expected particle coordinates using in-formation from four moments in time. This new hybrid algorithm is validatedagainst synthetic data and found to perfectly reproduce more trajectories thanother commonly used methods. Collisions are detected with 95% accuracy forparticles that move on average less than one tenth the distance to their nearestneighbor.
Keywords
Particle tracking · Image analysis · Particle collisions
The organization and coalescence rate of inertial droplets in turbulent envi-ronments is of great importance in a number of natural and industrial sys-
This material is based upon work supported by the National Science Foundation underGrant No. 1605195.Reece V. KearneyCornell University 124 Hoy Road, Ithaca NY 14850Tel.: (845)-802-3001E-mail: [email protected] P. BewleyCornell University RV Kearney, GP Bewley tems. In combustors, the size distribution of droplets affects the rate of evap-oration of fuel and thus combustion efficiency, as in Betelin et al. (2012). Inclouds, droplets must grow large enough in order to fall as rain. Inabilityto predict droplet growth across the so-called “size gap” (Xue et al., 2008;Grabowski and Wang, 2013), the regime of droplet sizes from about 10 mi-crons to 100 microns in which droplets grow faster than condensation anddifferential gravitational settling can explain, is one of the greatest uncertain-ties in climate and weather modeling (Devenish et al., 2012).Inertial droplets and their near-collision and coalescence have been stud-ied extensively. Analytical and numerical work predicted the growth rate ofdroplet populations. Smoluchowski (1916) introduces a stochastic collectionequation that includes the notion of a collision kernel, a measure of the col-lision rate. Saffman and Turner (1956) shows how the collision kernel canbe expressed for droplets in turbulence. Maxey and Corrsin (1986) describeshow particles settling in turbulence cluster into the straining regions of theflow. Direct numerical simulations (DNS) of inertial particles have been in-strumental in showing that turbulence induces clustering as a function ofthe particle inertia and Reynolds number of the flow and can enhance col-lision rates (Reade and Collins, 2000; Chen et al., 2006; Salazar et al., 2008;Wang et al., 2000; Ayala et al., 2008; Ireland et al., 2016a,b; Olivieri et al.,2014), and many of these findings have been successfully validated in experi-ments (Sumbekova et al., 2017; Obligado et al., 2014). Combined experimen-tation, simulation, and modeling have also illuminated the behavior of in-ertial particles sedimenting under the influence of gravity (Lapp et al., 2012;Aliseda et al., 2002; Good et al., 2014; Salazar et al., 2008; Kawanisi and Shiozaki,2008; Ireland and Collins, 2012; Siewert et al., 2014)Experiments measuring collisions to complement these results are lacking.Three-dimensional Lagrangian particle tracking (3D LPT) is used routinelyin research to measure fluid motions (Xue et al., 2008; G¨ulan et al., 2012;Holzner et al., 2008; Hoyer et al., 2005; Virant and Dracos, 1997; Ott and Mann,2000), and the motions of non-tracer particles in fluid (Lapp et al., 2012;Aliseda et al., 2002; Bewley and Saw, 2013; Salazar et al., 2008), but theseworks do not give direct measurements of collisions from 3D LPT.This deficit in experimental data is partly due to the difficultly of measuringthe motions of particles when they are close together (Ouellette et al., 2006).To avoid this difficulty, researchers have successfully conducted experimentsto characterize the preferential concentration and collision rate of populationsof droplets by other means (Bateson and Aliseda, 2012; Duru et al., 2007).Researchers have also recorded the coalescence outcome of collisions betweentwo droplets given the Weber number and impact geometry (Qian and Law,1997), but these studies do not explain how the dynamics of droplets areaffected by background flow or provide a statistical description of collisionrates given ambient flow conditions.Schanz et al. (2016) developed a family of new tracking algorithms knownas Shake the Box (STB). Whereas older methods search for particles in cameraimages independent of the tracking phase, STB methods use the time history agrangian tracking of colliding droplets 3 of each measured trajectory to inform the search for particles in subsequentimages. The predicted location of each particle is projected onto a cameraimage and “shaken” to optimize the difference between the measured pixelintensity and the anticipated effect of the particle that has moved to thatlocation. This allows trajectories to come much closer together. The emphasisof this algorithm has been on increasing the maximum particle seeding densitythat can be tracked, and to our knowledge, the methods have not been appliedexplicitly to study the motions of colliding particles.We believe that the only experiment that directly measured collisions ofdroplets in a turbulent environment is reported in Bord´as et al. (2013). Thesoftware used in this experiment, however, required user intervention to discernwhether or not collisions had occurred, and the uncertainty in the result is notclear. Ultimately, to understand the near-contact dynamics of particles, a newparticle tracking tool is necessary.
Particle tracking can be broken into three steps: the identification of parti-cles in two dimensional images on each camera, stereoscopic reconstruction ofthree-dimensional coordinates from several simultaneous camera images, andthe building up of 3D trajectories in time to determine what particle is thesame from frame to frame. In this paper we describe advancements in thefirst and third steps of this process. We do not discuss direct improvements tothe second step of particle tracking, though the process of stereoscopic recon-struction of particle coordinates can benefit from more precise two-dimensionalparticle identification and particle size measurement. Additionally, erroneouslyconstructed three-dimensional particles or “ghost particles” that result fromerrors in the reconstruction process can be more easily identified with betterparticle tracking.2.1 Particle IdentificationIn typical particle image velocimetry (PIV) and LPT applications, particles aresmall; their diameters are on the order of a few pixels or less. These methodsare often used in order to measure the motions of the flow, so minimizing theparticle size in order to make them better tracers of the flow is important.Because we are interested in inertial particles, we instead take advantage of ahigher resolution regime of particle images wherein particles have tens of pixelson their circumference. This allows us to use the particle size as a considerationto assess whether or not a collision has occurred during the collision detectionphase.The primary difficulty of differentiating between particles at small separa-tions in a two-dimensional image is identifying the distinct projections withina body as well as identifying their positions and sizes in pixel space. Consider a
RV Kearney, GP Bewley spherical particle at some position in three-space x w . This particle is recordedby a camera and a two-dimensional projection is created on the sensor. Theprojection of the particle leaves a circular spot. When two or more particlesare close together, individual spots can result from several overlapping circularprojections. We will call this agglomerated bright spot a single body.The new method discussed here uses only the edge data from each body.Many methods exist for finding edges in digital images (Canny, 1983; Deriche,1987; Sobel, 1968; Trujillo-Pino et al., 2013). We used the following process:1. Apply a Gaussian blur to the image to reduce the effect of noise (Kobayashi et al.,1991).2. Calculate the gradient of pixel intensity of the image.3. Find the pixel-accuracy boundary of the body using the coordinates ofpixels above some brightness threshold, T .4. For each pixel-accuracy boundary point, interpolate on the grid in thedirection of the gradient to find the subpixel location where intensity fallsbelow threshold T .Given the edge coordinates of a body, the task is to find the set of circlesthat best describes the body. This is a combinatorially complex process, andso we must formulate a reasonable way to continue. The algorithm we pro-pose, which we name “Pratt-Walking”, is described pictorially in Figure 1 andproceeds as follows.1. Consider a set of edge data that consists of N two-dimensional coordinates.2. A contiguous set of Q < N edge data points are selected and a circle is fitto the data using the method of Pratt (1987).3. This process repeats for all N sets of Q contiguous data points, yielding N circles (see 1, Top right ).4. Erroneous circles are discarded (circles with centers outside the bound-ing rectangle of the body and circles whose edges extend far outside thebounding rectangle of the body; see 1,
Bottom left ).5. The remaining circles are binned according to their center positions witha bin-width of 1 pixel.6. Contiguous sets of bins that represent a number of circles greater than somethreshold H are associated with populations of circles that have similarcenter coordinates. The properties of these circles are averaged to generatea reduced family of circles (see 1, Bottom right ).7. A residual is calculated that gives the average distance from the edge dataand the nearest edge of a circle in the reduced family. If this residual issmaller than some threshold G , then the solution is accepted.2.2 Particle TrackingGiven the coordinates of a set of particles in space and time, particle track-ing refers to the process of identifying which particle is the same from frame agrangian tracking of colliding droplets 5 Fig. 1
Top left:
An example of two synthetic particle images that are overlapping by 60%
Top right:
Circles fit by the Pratt method (Pratt, 1987) to contiguous data points along theboundary of the body ( Q = 11) Bottom left:
The circles that remain after those with centersoutside the body were discarded
Bottom right:
The results of averaging the properties offamilies of circles that have center coordinates within the same peak of the accumulatormatrix to frame. Here, we refer to a linked set of coordinates in time as a trajec-tory. When there are many particles during each of many time steps, con-structing trajectories is an NP-hard multidimensional assignment problem(Veenman et al., 2003). To avoid this computational difficulty, tracking algo-rithms use only a small subset of particle coordinates at a time (Ouellette et al.,2006).When particles are far apart from one another relative to how far theytravel, tracking their positions is easy since the number of physically realizabletrajectories is small. As the particle density increases or the particles movefaster, holding all else constant, the number of physically reasonable choices forthe next point in each trajectory increases. Regardless of the tracking methodchosen, tracking performance degrades as particles move faster relative to thedistance between them (Ouellette et al., 2006). If measurements of collisionsare to be taken, it is impossible to avoid the problem of particle coordinatescoming arbitrarily close together, and so it is necessary to form a robust,quantitative metric to determine the best match from the identified particlecoordinates in the next time step for each trajectory. Collision detection willtake place in a distinct step prior to the tracking process, as described inSection 2.3.Consider a trajectory t n that has been successfully tracked until position x n , i at time τ i . At time τ i , t n also has some estimate for its velocity v n , i .The algorithm seeks to extend this trajectory into the next time step τ i +1 bycalculating a penalty function associated with each set of particle coordinatesthat are known at time τ i +1 . This penalty function consists of the weightedsum of three components: RV Kearney, GP Bewley P = w d + w d + w d (1)Where the w i are weightings. The d m are defined for t n matching withparticle coordinates at position y j , i+1 : d = | x n , i − y j , i+1 | (2) d = | x n , i + v n , i dτ − y j , i+1 | (3) d = | y j , i+1 + v n , i+1 dτ − y k , i+2 | (4)where dτ = τ i +1 − τ i . y k , i+2 is chosen from all possible coordinates ofparticles at time τ i +2 such that d is minimized.These three d m correspond to three tracking methods described in Ouellette et al.(2006), so we call it a “hybrid method.” d is the distance between the cur-rent trajectory’s coordinates and the coordinates of a candidate particle at thenext time step. It is identical to the tracking cost of the “Nearest Neighbor”heuristic. d is the difference between the predicted and observed position ofparticle coordinates at the next time step. It is identical to the tracking costof the “Minimum Acceleration” method. d is the difference between the pre-dicted and observed position of particle coordinates in two time steps. It isidentical to the tracking cost of the “Four Frame: Best Estimate” with an es-timated acceleration of zero at all times. Determining the optimal weights totrack a given population of particles is not trivial, but we will show (in Section3.3) that even for unoptimized, reasonably chosen weights the hybrid methodoutperforms other methods.We also apply a threshold on the maximum distance a particle can travelin one frame: | x n , i − y j , i+1 | < M (5)Excluding matches based on the maximum distance a particle can travelyields a reduced set of potential matches. We use velocity estimates based ontwo-point, forward finite differences. At the beginning of trajectories, it is notpossible to formulate an estimate for the velocity using a finite difference, so weset d = M/
2. For this reason, trajectories achieve substantially higher valuesof P , the penalty function, near their beginning compared with trajectoriesconsisting of two or more points The optimal pairs of trajectories and candi-dates are those that minimize the sum of P across every trajectory, and theycan be found with the Munkres Algorithm (Bourgeois and Lassalle, 1971).2.3 Collision DetectionConsider two spherical droplets of radius r and r that collide and coalesce.Neglecting fragmentation, the radius of the daughter droplet is r = ( r + r ) .In addition, neglecting losses to heat and external forces, the momentum of agrangian tracking of colliding droplets 7 the center of mass of the system must remain the same. These are the criteriawe used to determine whether or not a collision occurred.For each trajectory at time τ i , before the track is extended into frame τ i +1 as described in the previous section, the anticipated path of the trajec-tory is constructed. The algorithm then searches for intersections and near-intersections of these predicted paths. If one is found, the algorithm searchesfor a daughter at time τ i +1 at the location predicted by conservation of mo-mentum. The algorithm also searches for the daughter at time τ i +2 , again atthe location predicted conserving momentum. If all these conditions are sat-isfied, then a collision is deemed to have taken place, and a link is formedbetween the two parent trajectories and the daughter trajectory. It is onlynecessary to check these criteria for trajectories that are closer than some dis-tance d < M + r + r , where r and r are the radii of the particles. Whereasprevious tracking algorithms have produced particle tracks that have only onebeginning and one end, this allows trajectories to converge so that a trajectorymay have many beginnings, each associated with parents, and one end, asso-ciated with the ultimate daughter. Droplet bursting can also be detected byfeeding particle coordinate data in reverse-chronological order. This method issuited for detecting collisions between liquid droplets at low Weber number. r located some distance dr away from the particle center is given by I ( dr ) = I max − ( I max − I min ) dr r (6) I max was chosen to be 180 and I min was chosen to be 10, consistent with an8-bit image with good contrast between particle and background. Zero-meanGaussian noise was added to the image to test the robustness of the algorithm.The standard deviation of the noise was set as a percentage of the maximumbrightness. Sample images of the synthetic particle at two different noise levelscan be seen in Figure 2.Figure 3 shows the results of the algorithm on an image containing a singleparticle. The results of the algorithm described here are compared to a circularHough Transform as implemented in the MATLAB function imfindcircles . Wealso compare against a simple method that searches for local maxima abovethreshold T in the pixel brightness distribution after a Gaussian blur. Thismethod does not achieve sub-pixel precision and is intended as a benchmarknaive and fast process The results shown are averaged over 10,000 trials. In RV Kearney, GP Bewley
80 100 120 x (pixels) y ( p i x e l s )
80 100 120 x (pixels) y ( p i x e l s ) Fig. 2
Left:
Synthetic particle image in 1% noise conditions.
Right:
Synthetic particleimage in 20% noise conditions. each trial, the synthetic particle image is generated in the center of a 200 by200-pixel frame and offset by a random, sub-pixel amount in both the x and y-directions. Pixel intensities are bounded such that they can only achieve valuesbetween 0 and 255 (consistent with 8-bit images). We use the following valuesfor the algorithm parameters (see Section 2.1 for definitions): T = 55; Q =11; H = 11; G = 1. The standard deviation of the Gaussian filter is √ P ID ( n ) as the probability of the algorithm detecting that thereare n particles in the image. The accuracy of finding the center of the particleis given by dxr p , where dx is the distance between the measured particle centerand its actual center and r p is the actual particle radius. The accuracy ofsize estimation is given by dr p r p where dr p is the difference between the actualparticle radius and the measured particle radius.Pratt-Walking identifies a wider range of particles of varying size withgreater reliability and superior center-finding precision than either of the othermethods considered at all noise levels. The Hough method is able to give betterestimates for particle sizes than Pratt-Walking for a small regime of particlesize (between 4 and 5 pixels in radius) at low noise conditions only. The localmaximum method is also able to more reliably identify particles with radiismaller than 3 pixels than Pratt-Walking at the lowest noise conditions, butits error in measuring the particle center increases steeply in this regime.Note that the radius of the smallest particle identifiable by Pratt-Walkingdecreases with increasing noise, which may be counterintuitive. This is be-cause the method relies on data at the edge of the body, where the true pixelintensity should be close to zero. While the noise added is Gaussian with amean value of zero, a pixel’s minimum value saturates at zero, so the averagepixel intensity imparted by the noise near the boundary is positive. This alsoleads to systematic overestimates in the size of the particle, as shown in thebottom of Figure 3.Pratt-Walking performs worst in terms of particle identification, with poorerperformance than the Hough Transform method for particles larger than 8 pix-els in radius for moderate noise conditions (see 3: Top ). This is a result of the agrangian tracking of colliding droplets 9 r p (pixels) P I D ( ) PrattHoughMaximum
Method r p (pixels) -3 -2 -1 x / r p Noise Level r p (pixels) -0.4-0.200.2 r p / r p Fig. 3
Top:
The probability of correctly identifying a particle in a simulated image ofan isolated particle for different noise levels and particle sizes.
Middle:
The average errorin measuring the center location of the particle.
Bottom:
The average error in measuringthe particle radius. A negative value indicates an overestimate of particle size. The localmaximum method does not produce an estimate of particle size.
Legend:
Symbols corre-spond to the identification method used. Color corresponds to the standard deviation ofthe Gaussian noise added to the images. In low-noise conditions, the center-finding errorfor Pratt-Walking is about 20 times smaller than the local maximum method and about 10times smaller than the Hough Transform method. choice of Q and H . Increasing Q and H would give better estimates for all threeof the metrics considered in Figure 3 for large particles, but would degradeperformance for small particles. In general, care must be taken in selecting thevalues for Q , H , and T based on the size and brightness of the particles beingtracked.Both the Hough Transform and Pratt-Walking exhibit radius error thatchanges as a function of the particle size and noise conditions. The near-constant error of Pratt-Walking at low-noise conditions is due to the choice of T , which was selected to be high enough to ignore ambient noise throughoutthe image, but this causes an underestimate of particle size. As discussedabove, Pratt-Walking overestimates the size of small particles because noiseenergizes the intensity of bodies at their border. All of these errors can beeasily corrected post-measurement with calibration on particles of known size.Note that the local maximum method does not produce an estimate for theparticle size. When there is one local maximum in pixel intensity in a body, a natural estimate for the particle size is the square root of the area of the bodyabove threshold T , but when multiple maxima exist, additional segmentationmethods are necessary to produce a size estimate, and we do not consider suchmethods here.3.2 Two-Particle Identification TestNext, the algorithm was tested on synthetic data that contains the projectionof two particles that overlap, a frequent occurrence in densely seeded flows.The overlap is defined as Ov = r + r − d r , r ) , (7)where d is the center-to-center distance between the two particles, and r and r are the radii. The overlap is defined such that it is equal to 0 when theparticles are just touching edges with their centers entirely outside one anotherand 1 when one particle completely obscures the other. r is held fixed at 10pixels while r is varied so that different radius ratios can be achieved.The results for this test are shown in Figure 4 again compared to a HoughTransform and the local maximum methods. The results for each set of condi-tions were averaged over 10,000 trials. Noise was fixed at 1% of the maximumpixel intensity.Here, the center error was the mean error averaged between both particles.The radius error was averaged in the same way. We used the following valuesfor the algorithm parameters: T = 55; Q = 8; H = 8; G = 2. The standarddeviation of the Gaussian filter was √ r r =0 .
4. Even in this case, while the reliability of Pratt-Walking is lower, it stillmaintains superior center-finding accuracy. The failure of Pratt-Walking forsmall radius ratios is due to insufficient information on the border of the bodythat is associated with the smaller particle. In general, Pratt-Walking requires Q + H > N p where N p is the number of edge data points associated with asingle particle. For images of small particles with very little noise, the localmaximum method works best, and indeed this method still has a place asa useful tool for characterization of flows seeded with tracer particles at lownumber densities. It is less suited for studies that require precise estimation ofparticle location and size.The radius error is positive for Pratt-Walking and the Hough Transformmethod at all values of overlap except at high overlap of very unlike sizedparticles, indicating a consistent underestimate of size. The magnitude of theerror is smaller for Pratt-Walking than the Hough Transform for each radiusratio. agrangian tracking of colliding droplets 11 P I D ( ) PrattHoughMaximum
Method -3 -2 -1 x / (r + r ) r /r =1.0r /r =0.8r /r =0.4 Radius Ratio
Overlap r / (r + r ) Fig. 4
Low noise (1%) cases.
Top:
The probability of correctly identifying two particlesin a simulated image of two overlapping particles.
Middle:
The average error in finding thecenter locations of the particles.
Bottom:
The average error in finding the particle sizes. Thelocal maximum method does not produce an estimate of particle size.
Legend
Symbols andcolor correspond to the identification method used. Line type corresponds to the ratio ofthe radii of the two particle images. Pratt-Walking is able to able to identify particle pairimages at overlap at least 40% higher than the Hough Transform method. The magnitude ofthe center-finding error for Pratt-Walking is at worst the same size as the error for the othertwo methods. At best, the center-finding error for Pratt-Walking is 15 times smaller thanthe error of the Hough Transform method and 20 times smaller than the local maximummethod.
Pratt-Walking, the Hough Transform, and the local maximum methods arealso compared at moderate noise levels (20% max pixel intensity). The resultsare shown in Figure 5.Noise does not substantially degrade the probability of identifying the cor-rect number of particles for Pratt-Walking for
Ov < .
6, though it does dra-matically harm the performance of the local maximum method, pushing theprobability of correct identification below 0 . r r = 0 . Ov < .
5, thoughPratt-Walking maintains functionality through higher overlap values. NeitherPratt-Walking nor the Hough method exhibit errors over 0.1 for regimes in P I D ( ) PrattHoughMaximum
Method -1 x / (r + r ) r /r =1.0r /r =0.8r /r =0.4 Radius Ratio
Overlap r / (r + r ) Fig. 5
Moderate noise (20%) cases.
Top, Middle and Bottom: as in Figure 4. The accuracyin center-finding and size-estimation for all three methods is decreased by the addition ofnoise. Pratt-Walking is able to able to identify particle pair images at overlap about 30%higher than the Hough Transform method. The local maximum method does not reliablyidentify particle pair images at any overlap at this noise level. The magnitude of the center-finding error for Pratt-Walking is about the same as the Hough Transform method andbetween 1.5 and 5 times smaller than the local maximum method. which the algorithm correctly identifies particles more than half the time. Thelocal maximum method exhibits errors above 0.1 for all values of the overlap.Noise also damages estimates for the particle sizes for the Hough methodand Pratt-Walking, pushing radius estimates farther from 0, especially at highoverlap. Errors associated with Pratt-Walking remain below 0.2 in absolutevalue for all overlap values and all radius ratios.The local maximum method took the least time to process, and the Pratt-Walking and the Hough methods took about 2.5 and 3.8 times as much timeas the local maximum method, respectively.3.3 Particle TrackingIn order to quantify the performance of our new hybrid particle-tracking al-gorithm, we tested it on synthetic data as in Ouellette et al. (2006). The datawere generated by a particle-laden turbulent flow solver called HiPPSTR anddeveloped by Ireland et al. (2013). For an isotropic turbulence at an average agrangian tracking of colliding droplets 13 R λ = 52 in a periodic cubic domain of length 2 π , statistical stationarity ismaintained through deterministic forcing. 6400 tracer particles were presentwithin the volume at all times. Tracer particles have no inertia, and so have aStokes number St = τ p /τ η = 0 where τ p is the relaxation time of the particleand τ η is the Kolmogorov time scale.A measure for the difficulty of tracking a population of particles is givenby ξ which is defined as ξ = δxδd , (8)where δx is the average displacement of a particle from one frame to the nextand δd is the average separation between a particle and its nearest neighbor.The quality of tracking can be measured by E track , which is defined as E track = T imperfect T total , (9)where T imperfect is the number of tracks the algorithm failed to produce per-fectly, and T total is the total number of tracks in the data set. A perfectlytracked trajectory must begin at the same time as an actual trajectory butneed not end at the same time. This way, an initial break in a trajectory isnot penalized more than subsequent breaks. Additionally, a perfectly trackedtrajectory must contain no spurious points. This measure is identical to theone considered in Ouellette et al. (2006). ξ is altered by under-sampling the particle positions so that the averagedisplacement of particles seen by the tracking algorithm is increased, resultingin trajectories that span the same spatial distance but have less informationalong that space. When a particle passes the periodic boundary and returns onthe opposite side of the simulation volume, we consider it to be a new trajectoryfor purposes of calculating ξ and E track . This results in about 86,000 uniquetrajectories.We use a simple finite difference to estimate all particle velocities. Ratherthan using the Munkres algorithm to resolve conflicts between trajectories,we use a simple first-come first-serve heuristic where the candidate with thebetter (smaller) score for the penalty function obtains the next point in thetrajectory. We use the following weights in calculating the penalty function(see Equation 1): w = 1; w = 5; w = 4. We set M = 0 . f s where f s is thesampling frequency.Figure 6 shows the mean values of E track for the hybrid algorithm comparedto the best-performing and worst-performing methods from Ouellette et al.(2006); these are respectively the Nearest Neighbor (NN) and 4 Frame: BestEstimate (4BE) heuristics. Note that our implementation of 4BE extendedall trajectories that already contained two or more sets of coordinates beforeextending new trajectories, which enhances the performance of the method.Our results for NN and 4BE differ from those in Ouellette et al. (2006) becauseour methodology for altering and computing ξ differs.The hybrid method achieves lower values of E track than the NN or 4BEmethods for all ξ . The NN heuristic constructs imperfect trajectories whenever -1 -5 -4 -3 -2 -1 E t r a ck NN4BENew Hybrid
Fig. 6
We measured the performance of particle tracking methods with the tracking error, E track , which is the fraction of imperfect tracks in the dataset, and which increases as afunction of ξ , which measures how far particles travel relative to their distances from eachother. We compared our hybrid method to select algorithms presented in Ouellette et al.(2006), of which the Nearest Neighbor (NN) algorithm was the worst-performing and the 4Frame: Best Estimate (4BE) method was the best-performing. the separation distance between two trajectories is smaller than the distancetraveled between time steps. This kind of event increases in frequency withincreasing ξ .The hybrid method outperforms the 4BE heuristic in two ways: by betterinitiating trajectories and by better continuing longer trajectories. For newtrajectories, 4BE defaults to the NN heuristic, the worst-performing algorithm.In contrast, the hybrid method uses three frames of data to initiate trajectories,as discussed in Section 2.2.We ran our tracking tests while enforcing the same perfect initiation foreach of the tracking algorithms, so that every trajectory began with two correctsets of coordinates. As expected, every algorithm performed better. The differ-ence in E track between the perfect-initiation test and the standard-initiationtest gives the initiation error, and the final value of E track gives the error dueimperfect continuation of longer trajectories. At tracking difficulty ξ = 0 . E track for 4BE was reduced by 0.15 to a final value of 0.29. E track for thehybrid method was reduced by 0.011 to a final value of 0.014. In other words,we found that even with perfect initiation, the hybrid method outperformsthe 4BE method by an order of magnitude. The proposed hybrid method per- agrangian tracking of colliding droplets 15 forms better in part because it incorporates the best features of each method,while excluding the large uncertainty inherent in acceleration estimates usedby 4BE.3.4 Collision DetectionSimilar to E track , we define metrics for the quality of collision detection C g = N g N total , (10) C b = N b N total , (11)where N total is the total number of collisions that occur in the subdomain, N g is the number of collisions correctly identified, and N b is the number of falsecollisions identified (that is, the algorithm declares a collision occurred whenin fact none did).The performance of the algorithm on DNS data with particle collisions isshown in Figure 7. HiPPSTR was modified to incorporate a collision algo-rithm (Li Sing How and Collins, 2020) and simulate coalescence assuming aunit collision efficiency, meaning that every particle crossing trajectory resultsin coalescence, whereby mass and momentum (but not energy) are conserved.6400 inertial particles of identical size were allowed to reach statistical equi-librium over five large-eddy turnover times before the collision algorithm wasactivated, whereupon coalescence events were allowed to occur. The particlesalso had an initial Reynolds number Re p = r u η ν = 0 .
16 where r is the initialparticle radius before any collisions took place, u η is the Kolmogorov veloc-ity scale, and ν is the kinematic viscosity of the carrier fluid. At the end ofthe simulation, 6308 particles remained. The particles all had an initial ra-dius of 0.0042 and St = 0 .
1. To correctly identify a collision, we require thatthe collision occur at the correct time and location. Because coalescence isinstantaneous in the simulation, collisions always occur between time steps,so we only require that the collision time be identified to within the intervalbetween the preceding and following time step. We require a spatial accuracyof no more than 0 . M . We break the domain into 27 subdomains of equalsize without overlap in order to decrease the computational expense, whichgoes as the square of the number of contemporaneous trajectories. Addition-ally, we consider each daughter particle that results from a collision to be thestart of a new trajectory for purposes of calculating E track . This results inseveral hundred trajectories per subdomain over the course of the simulation.We perform a bootstrap process with 10,000 iterations on the 27 subdomainsat each sampling rate to find the mean and 95% confidence intervals of themean at different values of tracking difficulty. We use the following weights incalculating the penalty function: w = 1; w = 5; w = 4. We set M = 0 . f s .For values of ξ smaller than 0.2, collisions are found nearly perfectly. Nofalse collisions are detected in any case, and C g remains nearly constant at -1 -4 -3 -2 -1 E track C g C b Fig. 7
We applied our algorithms to synthetic data to quantify the increase in trackingquality compared to other methods and to measure the accuracy in collision identification.The bars show 95% confidence intervals for the mean of each quantity. When ξ ≈ . C g , is equal to the numberof false collision detections, measured by C b . E track and ξ are reduced here compared toFigure 6 because the decomposition of the domain into 27 subdomains artificially lowers ξ without changing the motion of the synthetic particles. ξ increases, C g decreases smoothly to about 0.01 at ξ = 1. When ξ issmall, the few errors in collision detection lead to a corresponding, small butfinite tracking error.Note that the values of E track here are somewhat lower at the same value of ξ than in the previous simulation with tracer particles. This is a consequenceof subdividing the domain, which has two effects: less information is availablenear the edge of each subdomain, resulting in a degradation of E track , and thespace outside a subdomain is regarded as empty, causing a decrease in ξ .These data suggest that detecting collisions imposes a more stringent re-quirement on the tracking difficulty than does simply tracking particle motion.At ξ = 0 .
4, only 2% of the measured trajectories are imperfect, but more thanhalf of the collisions are missed. agrangian tracking of colliding droplets 17
We introduce a new particle identification method, “Pratt-Walking”, and hy-brid tracking algorithm for three-dimensional LPT. The former is able to dif-ferentiate between two overlapping particles to over 80% overlap, is robust tonoise, and outperforms the circular Hough Transform as well as local maximummethods. The hybrid algorithm is capable of tracking particles to separationsapproaching collision, tracks an order of magnitude more reliably than the 4Frame: Best Estimate, and is able to measure collision rates with 95% accuracyfor tracking difficulty ξ < .
2. Detecting collisions with high certainty imposesa greater restriction on the maximum tracking difficulty than does high-fidelitytracking of non-colliding trajectories; where maintaining E track < . ξ < .
7, maintaining C g > . ξ < .
2. This work represents thefirst successful effort to develop a systematic tool for measuring collision ratesof particles suspended in fluids.The hybrid tracking algorithm presented here is generalizable and modular.By relying on minimizing a penalty function that is the weighted sum of manyelements, the method maintains flexibility for tracking distinct populationsof particles. While here we use only three quantities in the penalty function(based on the positions of identified particles in four frames of data), thiscan be extended to any number of characteristics of the particles (size, shape,etc.) and an expectation for how they should evolve with time with marginalincrease in computation time for each addition. It is necessary to constructthis objective function as a robust metric for differentiating between particlesas they approach contact in order to find collisions between trajectories.
Acknowledgements
We would like to thank Melanie Li Sing How and Lance Collins ofCornell University for providing DNS data of inertial coalescing particles in turbulence fortesting the tracking algorithm presented here.
References
Aliseda A, Cartellier A, Hainaux F, Lasheras JC (2002) Effect of preferentialconcentration on the settling velocity of heavy particles in homogeneousisotropic turbulence. Journal of Fluid Mechanics 468:77–105, DOI 10.1017/S0022112002001593Ayala O, Rosa B, Wang LP, Grabowski WW (2008) Effects of turbu-lence on the geometric collision rate of sedimenting droplets. Part 1. Re-sults from direct numerical simulation. New Journal of Physics 10, DOI10.1088/1367-2630/10/7/075015Bateson CP, Aliseda A (2012) Wind tunnel measurements of the preferen-tial concentration of inertial droplets in homogeneous isotropic turbulence.Experiments in Fluids 52(6):1373–1387, DOI 10.1007/s00348-011-1252-6Betelin VB, Smirnov NN, Nikitin VF, Dushin VR, KushnirenkoAG, Nerchenko VA (2012) Evaporation and ignition of droplets in combustion chambers modeling and simulation. Acta As-tronautica 70:23–35, DOI 10.1016/j.actaastro.2011.06.021, URL http://dx.doi.org/10.1016/j.actaastro.2011.06.021
Bewley, Gregory P and Saw, Ewe-wei (2013) Observation of the sling effectNew Journal of Physics 15, DOI 10.1088/1367-2630/15/8/083051Bord´as R, Roloff C, Th´evenin D, Shaw RA (2013) Experimental determinationof droplet collision rates in turbulence. New Journal of Physics 15, DOI10.1088/1367-2630/15/4/045010Bourgeois F, Lassalle JC (1971) An extension of the Munkres algorithm for theassignment problem to rectangular matrices. Communications of the ACM14(12):802–804, DOI 10.1145/362919.362945Canny J (1983) Finding Edges and Lines in Images. PhD thesis, MassachusettsInstitute of TechnologyChen L, Goto S, Vassilicos JC (2006) Turbulent clustering of stagnation pointsand inertial particles. Journal of Fluid Mechanics 553:143–154, DOI 10.1017/S0022112006009177Deriche R (1987) Using Canny’s criteria to derive a recursively implementedoptimal edge detector. International Journal of Computer Vision 1(2):167–187, DOI 10.1007/BF00123164Devenish BJ, Bartello P, Brenguier JL, Collins LR, Grabowski WW, IjzermansRH, Malinowski SP, Reeks MW, Vassilicos JC, Wang LP, Warhaft Z (2012)Droplet growth in warm turbulent clouds. Quarterly Journal of the RoyalMeteorological Society 138(667):1401–1429, DOI 10.1002/qj.1897Duru P, Koch DL, Cohen C (2007) Experimental study of turbulence-inducedcoalescence in aerosols. International Journal of Multiphase Flow 33(9):987–1005, DOI 10.1016/j.ijmultiphaseflow.2007.03.006Good GH, Ireland PJ, Bewley GP, Bodenschatz E, Collins LR, Warhaft Z(2014) Settling regimes of inertial particles in isotropic turbulence. Journalof Fluid Mechanics 759:R3, DOI 10.1017/jfm.2014.602Grabowski WW, Wang Lp (2013) Growth of Cloud Droplets in a Tur-bulent Environment. Annu Rev Fluid Mech 45:293–326, DOI 10.1146/annurev-fluid-011212-140750G¨ulan U, L¨uthi B, Holzner M, Liberzon A, Tsinober A, Kinzelbach W(2012) Experimental study of aortic flow in the ascending aorta via Par-ticle Tracking Velocimetry. Experiments in Fluids 53(5):1469–1485, DOI10.1007/s00348-012-1371-8Holzner M, Liberzon A, Nikitin N, L¨uthi B, Kinzelbach W, Tsinober A (2008)A Lagrangian investigation of the small-scale features of turbulent entrain-ment through particle tracking and direct numerical simulation. Journal ofFluid Mechanics 598:465–475, DOI 10.1017/S0022112008000141Hoyer K, Holzner M, L¨uthi B, Guala M, Liberzon A, Kinzelbach W (2005) 3Dscanning particle tracking velocimetry. Experiments in Fluids 39(5):923–934, DOI 10.1007/s00348-005-0031-7Ireland PJ, Collins LR (2012) Direct numerical simulation of inertial particleentrainment in a shearless mixing layer. Journal of Fluid Mechanics 704:301–332, DOI 10.1017/jfm.2012.241 agrangian tracking of colliding droplets 19
Ireland PJ, Vaithianathan T, Sukheswalla PS, Ray B, Collins LR (2013) Com-puters & Fluids Highly parallel particle-laden flow solver for turbulence re-search. Computers and Fluids 76:170–177, DOI 10.1016/j.compfluid.2013.01.020, URL http://dx.doi.org/10.1016/j.compfluid.2013.01.020
Ireland PJ, Bragg AD, Collins LR (2016a) The effect of Reynolds numberon inertial particle dynamics in isotropic turbulence . Part 1 . Simulationswithout gravitational effects. Journal of Fluid Mechanics 796:617–658, DOI10.1017/jfm.2016.238Ireland PJ, Bragg AD, Collins LR (2016b) The effect of Reynolds numberon inertial particle dynamics in isotropic turbulence . Part 2 . Simulationswith gravitational effects. Journal of Fluid Mechanics 796:659–711, DOI10.1017/jfm.2016.227Kawanisi K, Shiozaki R (2008) Turbulent effects on the settling velocity ofsuspended sediment. Journal of Hydraulic Engineering 134(2):261–266, DOI10.1061/(ASCE)0733-9429(2008)134Kobayashi H, White JL, Abidi AA (1991) An Active Resistor Networkfor Gaussian Filtering of Images. IEEE Journal of Solid-State Circuits26(5):738–748, DOI 10.1109/4.78244Lapp T, Rohloff M, Vollmer J, Hof B (2012) Particle tracking for polydispersesedimenting droplets in phase separation. Experiments in Fluids 52(5):1187–1200, DOI 10.1007/s00348-011-1243-7Li Sing How M, Collins LR (2020) Direct Numerical Simulation of near contactmotion and coalescence of inertial droplets in turbulence. Manuscript inpreparationMaxey, Corrsin (1986) Gravitation Settling of Aerosol Particles in Ran-domly Oriented Cellular Flow Fields. Journal of the Atmospheric Sciences43(11):1112–1134Obligado M, Teitelbaum T, Cartellier A (2014) Preferential concentration ofheavy particles in turbulence. Journal of Turbulence 15(5):293–310, DOI10.1080/14685248.2014.897710Olivieri S, Picano F, Sardina G, Iudicone D, Brandt L (2014) The effect ofthe Basset history force on particle clustering in homogeneous and isotropicturbulence. Physics of Fluids 26(4), DOI 10.1063/1.4871480Ott S, Mann J (2000) An experimental investigation of the relative diffusionof particle pairs in three-dimensional turbulent flow. Journal of Fluid Me-chanics 422:207–223, DOI 10.1017/S0022112000001658Ouellette NT, Xu H, Bodenschatz E (2006) A quantitative study of three-dimensional Lagrangian particle tracking algorithms. Experiments in Fluids40(2):301–313, DOI 10.1007/s00348-005-0068-7Pratt V (1987) Direct Least-Squares Fitting of Algebraic Surfaces. ComputerGraphics 21(4):145–152Qian J, Law CK (1997) Regimes of coalescence and separation indroplet collision. Journal of Fluid Mechanics 331:59–80, DOI 10.1017/S0022112096003722Reade WC, Collins LR (2000) Effect of preferential concentration on turbulentcollision rates. Physics of Fluids 12(10):2530–2540, DOI 10.1063/1.1288515
Saffman PG, Turner JS (1956) On the collision of drops in turbulent clouds.Journal of Fluid Mechanics 1(1):16–30, DOI 10.1017/S0022112056000020Salazar JP, De Jong J, Cao L, Woodward SH, Meng H, Collins LR (2008)Experimental and numerical investigation of inertial particle clustering inisotropic turbulence. Journal of Fluid Mechanics 600:245–256, DOI 10.1017/S0022112008000372Schanz D, Gesemann S, Schr¨oder A (2016) Shake-The-Box: Lagrangian parti-cle tracking at high particle image densities. Experiments in Fluids 57(5):1–27, DOI 10.1007/s00348-016-2157-1Siewert C, Kunnen RPJ, Meinke M, Schr¨oder W (2014) Orientation statis-tics and settling velocity of ellipsoids in decaying turbulence. AtmosphericResearch 142:45–56, DOI 10.1016/j.atmosres.2013.08.011Smoluchowski M (1916) Drei vortrage uber diffusion, brownsche bewegungundkoagulation von kolloidteilchen. International Journal of Research in Phys-ical Chemistry and Chemical Physics 17:557–585Sobel I (1968) A 3 x 3 isotropic gradient operator for image processing. In:Stanford Artificial Intelligence ProjectSumbekova S, Cartellier A, Aliseda A, Bourgoin M (2017) Preferen-tial Concentration of Inertial Sub-Kolmogorov Particles. The roles ofmass loading of particles, St and Re. Physical Review Fluids 2(2):1–18, arXiv:1607.01256v1
Trujillo-Pino A, Krissian K, Alem´an-Flores M, Santana-Cedr´es D (2013)Accurate subpixel edge location based on partial area effect. Image andVision Computing 31(1):72–90, DOI 10.1016/j.imavis.2012.10.005, URL http://dx.doi.org/10.1016/j.imavis.2012.10.005http://dx.doi.org/10.1016/j.imavis.2012.10.005