Visualization of submerged turbulent jets using particle tracking velocimetry

Over the past few decades, advances have been made in using particle image velocimetry (PIV) and particle tracking velocimetry (PTV) for mapping of Lagrangian velocity and acceleration flow fields. With PIV, Lagrangian trajectories are not measured directly; rather, hypothetical trajectories must be constructed from sequences of Eulerian velocity snapshots. Because PTV directly measures actual trajectories, it provides distinct advantages over PIV, especially for trajectories with abrupt changes in direction. In this work, a novel particle tracking algorithm is described, then applied to track trajectories of tracer particles in submerged turbulent jets. The Reynolds numbers ranged from 1000 to 25,000, thereby covering laminar, transitioning-to-turbulence, and fully turbulent flow regimes. The novel particle tracking algorithm is designed to handle flows with very high particle concentrations, thereby resolving small-scale flow structures. Trajectories are tracked with high velocity gradients, sharp curvatures, cycloids, abrupt changes in direction, and strong recirculation—all of which are inaccessible via construction from PIV sequences. Most trajectories measured in this work are at least 500 camera frames (time steps) long, with many being more than 3000 frames long.

particles in submerged turbulent jets with Reynolds numbers from 1000 to 25,000, thereby covering laminar, transition-to-turbulence, fully developed turbulent flow regimes.
Various types of algorithms have been developed for Lagrangian particle tracking in turbulent flows: fuzzy logic by Hering et al. (1995); neural networks by Ge and Cha (2000), linear array optimization by Sethi and Jain (1987), and extrapolation searches by Adamczyk and Rimai (1988), Papantoniou and Dracos (1990), Shaffer and Ramer (1989), Ramer and Shaffer (1992), Singh et al. (1993), Willneff and Gruen (2002), Willneff (2003), Lüthi et al. (2005), Hoyer et al. (2005), Ouellette (2006), Ouellette et al. (2006), Bodenschatz et al. (2007), Shaffer (2012, 2013), and Gopalan et al. (2016). Dabiri and Pecora (2020) present an overview of PTV algorithms. The basis of the PTV algorithm in this paper was started in the late 1980s with the work of Shaffer et al. (1988) and Ramer and Shaffer (1992), and has been further developed by Shaffer (2012, 2013), and Gopalan et al. (2016) for flows of maximum particle concentrations. To produce Lagrangian PTV data, Shaffer et al. (1988) used a copper-vapor laser to multiply expose a high-resolution cathode-ray tube camera. Their algorithm used linear extrapolation of particle positions in two camera frames at t i and t iþ1 to search for the next image of a particle at t iþ2 in a circular search area (Fig. 1). Ramer and Shaffer (1992) found that the use of a higher-order kinematic models for extrapolation did not improve recognition of trajectories. Willneff and Gruen (2002), Willneff (2003) and Ouellette (2006), Ouellette et al. (2006), Bodenschatz et al. (2007), later arrived at similar conclusions. The following constraints were used to improve the accuracy of trajectory recognition: a minimum number of particle centroids was required for a valid trajectory; the maximum number of particle centroids along a valid trajectory is the number of exposures during a camera frame; the maximum distance between centroids along a trajectory is determined by the maximum velocity in the flow field and the exposure rate. These are similar to the constraints used by Sethi and Jain (1987), Papantoniou and Dracos (1990), Maas et al. (1993), Ouellette (2006), Ouellette et al. (2006), Bodenschatz et al. (2007), Willneff and Gruen (2002), Willneff (2003) and others.

Particle tracking
The tracking algorithm used in this study is a three-frame search method with iteratively increasing search areas. Each particle in a frame, i, at, t i , is tested as the first particle, P 1 , of a trajectory. An initial Nearest Neighbor Search Area (NNSA) in frame i þ 1 is defined around P 1 , which is illustrated in Fig. 2. The radius of the NNSA is determined by a search velocity, V search , and the time between camera frames, Dt. The radius of the NNSA is iteratively increased from a minimum value close to zero. Starting V search with a small value will find slowest particles first. Trajectories of the slowest particles represent the highest local particle concentrations. Recognizing and removing the trajectories with the highest local particle concentrations reduces computation time and improves recognition accuracy. The maximum radius of the NNSA is maximum velocity expected in a flow field multiplied by the time between frames.
In order to continue the trajectory recognition outside the NNSA, a Descendant Search Area (DSA) is formed with its center along a line projected from particle P 1 through each particle found in the NNSA. The length of the line is twice the distance between particles P 1 and P 2 (Fig. 3). Particles found inside the DSA are called nth generation descendants. The radii of the DSA, R acc min and R acc max , represent a minimum and a maximum translational acceleration, while the angle, / acc , represents an angular acceleration. The search parameters, R acc min , R acc max , and / acc , are also iteratively varied from a small value to a maximum value. The value of R acc min can be varied in the range of R 1À2 \R acc min \2R 1À2 . The value of R acc max can be varied in the range of 2R 1À2 \R acc min \3R 1À2 . However, it is rarely necessary to expand the DSA to its maximum, because particles will be found at much smaller values of the DSA. In Fig. 3, a candidate particle, P 3 , was found in the upper DSA, but not in the lower DSA shown. The next-generation DSA is then formed using particles P 2 and P 3 . This procedure repeats through future frames until a particle is not found in the DSA. It is not unusual to have missing images of a particle along its trajectory. This can be caused by particles blocking the view of other particles, or by casting shadows onto another particle. In this case, to handle missing particle images, an ''imaginary centroid'' is formed at the center of the DSA. The search then continues by forming the next DSA using the last found centroid and the imaginary centroid. If a particle is not found in the next DSA, the search is terminated. If the length of the trajectory is longer than a minimum value, the trajectory is defined as a ''candidate trajectory.'' When a candidate trajectory is recognized, its centroids are removed from consideration in future searches, which improves the speed and accuracy of trajectory recognition.
Many candidate trajectories can be grown from each particle. This forms a tree of candidate trajectories, as illustrated Fig. 4. This approach has been applied to a wide range of applications, from the flow of 30lm fluorescent particles in 1-cm-diameter micro-turbine blood pump spinning at 10,000 RPM (Kerrigan et al. 1993;Shaffer et al. 1992), to the flow of 60 lm FCC (fluid cracking catalyst) particles in a 18m tall, 30-cmdiameter riser of a circulating fluidized bed Singh et al. 1993;Shaffer 2012, 2013;Shaffer et al. 2013;Gopalan et al. 2016). These tracking algorithms were designed to handle missing particle images, crossing trajectories and conflicts in search areas. Conflicts occur when more than one particle image was found in the search area. Figure 4 shows particles found from a search starting with a first particle in the first frame, P 11 , where the first subscript represents the frame and the second subscript the particles in the frame. The longest trajectory is assumed to be the correct one, i.e., the branch P 11 -P 22 -P 33 -P 42 -P 51 in Fig. 4. Shaffer (2012, 2013) and Gopalan et al. (2016) further developed the tracking algorithm of Ramer and Shaffer (1992) to handle flows of very high particle concentrations, up to maximum particle packing. The accuracy of searches was improved by curve-fitting trajectories using the Cholesky curve-fit method. This negates steering of searches caused by jitter in the location of the centroids of particle images. Lüthi et al. (2005) showed that curve-fitting trajectories in the search process acted as a low-pass filter to remove noise and improve the accuracy of calculating velocity derivatives along trajectories.

Flows
The subject matter of the current study is the simultaneous schlieren and PTV imaging recorded in the near field of the submerged water jets. Fluorescent dye visualization was also done simultaneously and is described in Ibarra et al. (2020). Submerged jets were created in a glass water tank with dimensions ðL Â W Â HÞ of 240 cm Â 120 cm Â 120 cm at a water level of 105 cm. The exit diameter of the submerged jets was D ¼ 1:38 cm. The flows were seeded with 45lm silver-coated hollow ceramic spheres (Potter Industries Inc., AG-SF-20, 0.8 g/cm 3 ). For these particles, while experiencing the highest velocity gradients in these jet studies, the Stokes number was 0.01. This ensured that the seed particles accurately followed the fluid flow. High-speed videos of the submerged jets were taken at 2000 frames per second.
Four flow conditions were studied, with Reynolds numbers of 1000; 3000; 12,700; and 25,200. The Reynolds number, Re ¼ 4Q=pmD, is based on the flow discharge rate Q and the kinematic viscosity m ¼ 0:01 cm 2 /s of water at ambient temperature. This range covers laminar, transitional and turbulent flow discharge regimes. The corresponding complete PTV results as videos are available as supplemental material online at Videos (2020).
The schlieren system used two concave mirrors in a Z-configuration. Two additional flat mirrors were used to extend the total beam path 10 m from light source to cameras. Each concave mirror had a diameter of 45 cm and a focal length of 400 cm. The beam path was 12 off the normal to the PTV laser sheet in order to allow clear 90 access for the cameras. A light-emitting diode (LED) light source (Leica KL 1500LED) was used for illumination. The light beam was collimated using a matched achromatic doublet lens pair Search Area Fig. 1 Extrapolation and search areas used by Ramer and Shaffer (1992) (Thorlabs MAO:103030-A), a pinhole and a microscope objective. Both a single knife edge and two orthogonal knife edges were applied after the light beam is split by a cubical beam splitter. Further details of the schlieren system and its application may be found in Ibarra et al. (2020). Figure 5 shows a schlieren image of a submerged jet at Re ¼ 1000. The jet transitions to turbulence about six diameters downstream of the exit. High-speed video of the particle flow was recorded synchronously with a second camera for PTV analysis. The particle laden flow was illuminated by a 4-mm-  Fig. 4 Candidate tree of particles found in search areas initiated from particle P 11 thick laser sheet through the center of the jet from a 10W Argon-ion laser in TEM-00 mode. Particle trajectories were recognized with the algorithm of Shaffer (2012, 2013) and Gopalan et al. (2016) described above. Figure 6 shows recognized trajectories for the flow shown in Fig. 5 for 3000 frames. Trajectories are pseudo-colored according to the average velocity along the trajectory. A total of 56,400 trajectories were recognized in this flow. Most of the trajectories in the jet core are more than 1000 frames long, with many trajectories away from the jet core up to 3000 frames long. Figure 7 shows a closeup view of the trajectories in the region undergoing a transition-to-turbulence at Re ¼ 1000. In the shear   Figure 8a shows an example trajectory, recirculating in the shear layer at Re ¼ 1000. An 8 Â 8 pixels window is shown in yellow in Fig. 8a. PIV with interrogation window sizes of, for example, 16 Â 16 pixels, would not be able to resolve the fine scale recirculation structures in this trajectory. The velocity along a trajectory can be written in terms of a magnitude (speed), u, and the direction of the tangential unit vector e t : u ¼ ue t . The material acceleration along the trajectory can be also be written as a magnitude and direction of a tangential unit vector (Novaro and Scarano 2013) where Du/Dt is the magnitude of the material acceleration along a trajectory. Figure 8b, c shows the speed and acceleration, respectively, along the trajectory shown in Fig. 8a. To remove jitter in the measurement of particle centroids along trajectories, the speed shown is a moving average of 5-10 frames. The cusps of the cycloid in Fig. 8a correspond to the minima in the speed history in Fig. 8b. Clearly, as the particle reverses direction, its speed goes through the minimum value of zero at the cusps, becoming momentarily stationary. It is intriguing not to see comparable signatures in the acceleration history in Fig. 8c. However, a moving box-averaging of about 50 samples wide suggests negative excursions at the cusps of the cycloid in Fig. 8a. Figure 9 is a schlieren snapshot at Re ¼ 3000. Figure 10 shows 6704 trajectories found in 500 frames. More than 30,000 trajectories were found over 3000 frames, but when overlaid, are too dense to distinguish. Increased turbulence level with increasing Re hinders tracking particles for longer periods. Figure 11a shows an example trajectory at the edge of the jet, a cycloid, with abrupt changes in direction. Figure 11b, c shows the speed and acceleration along that path, respectively. The vanishing speed on the trajectory in Fig. 11b corresponds to the cusp of the cycloid while the second minimum point corresponds to the second sharp turn in the trajectory in Fig. 11a. The acceleration trace in Fig. 11c shows negative values at the cusps and substantially positive values at the middle. Figure 12 shows a schlieren snapshot at Re ¼ 12;700 and Fig. 13 shows an overlay of 9884 trajectories recognized in 500 frames. Figure 14a shows an example, Trajectory 35,905, which is 352 frames long. With the measured coordinates of the trajectory, the speed and acceleration along the trajectory can be readily   Fig. 14b, c. Here again, velocity and acceleration are moving-averaged over 5-10 frames, to remove jitter caused by uncertainties in the locations of particles. The particle on the sample trajectory in Fig. 14a experiences a rather high deceleration as it goes through the loop in the trajectory. Figure 15 shows a schlieren snapshot at Re ¼ 25;200 and Fig. 16 shows an overlay of trajectories recognized in 1000 frames. The trajectories outside the turbulent jet show a steady entrainment flow pattern. The pattern loses its coherence as the particles reach the edge of the jet. Figure 17a shows example Trajectory 477 at the edge of the jet, which is 999 frames long. Figure 17b, c also shows the speed and acceleration along Trajectory 477, respectively. The trajectory starts outside the jet in the entrainment region. As the particle meets the turbulent interface, it goes through rapid acceleration/deceleration, and nearly comes to rest before reversing its direction. Once in the turbulent flow, it undergoes large changes in both in speed and acceleration. Figure 18 shows a map of time-averaged magnitude of acceleration along trajectories at Re ¼ 12;700. This is direct measurement of the acceleration along actual particle trajectories from PTV-not an estimation based on construction of hypothetical trajectories from series of PIV maps. For display of the acceleration in Fig. 18, the field-of-view is divided into grid of 64 cells by 26 cells. The magnitude of acceleration of particle trajectories as they pass through each cell is averaged for 3000 frames.
Material acceleration is important because it can be used to calculate instantaneous pressure fields by integration of the pressure gradient term in the Navier-Stokes equations for incompressible flow (Lui and Katz 2006;Novaro and Scarano 2013;van Gent et al. 2017) For incompressible flows of high Reynolds numbers, away from boundaries, the material acceleration is dominant and the viscous term can be neglected, so In this work, a novel particle tracking algorithm was described and applied to track Lagrangian trajectories of tracer particles in submerged turbulent jets. Lagrangian trajectories with high velocity gradients, sharp curvatures, cycloids, abrupt changes in direction, and strong recirculations were tracked-all of which are all inaccessible via reconstruction from PIV sequences. The algorithm is capable of tracking flows with very high particle concentrations, thereby resolving small scale flow structures. Spatial resolution is on the order of a pixel. To demonstrate the particle tracking algorithm, Lagrangian trajectories of tracer particles were tracked in submerged turbulent jets with Reynolds numbers in the range of 1000 to 25,000, thereby covering laminar, transitioning-to-turbulence, and fully turbulent flow regimes. Most trajectories were at least 500 camera frames (time steps) long, with many being more than 3000 frames long. This allows direct calculation of the Lagrangian acceleration along real trajectories. Lagrangian maps of real trajectories can be used to calculate instantaneous pressure fields, either solving the Poisson equation or direct integration of the pressure gradient term in the Navier-Stokes equations. Instantaneous pressure maps are valuable in a wide range of aerodynamics and fluid dynamics applications. Funding USDOI Bureau of Environmental Enforcement and USDOE Energy Policy Act Program.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.