1 Introduction

Fluid turbulence has been studied mostly from an Eulerian framework. One of the reasons is the difficulty in measuring the Lagrangian behavior of turbulence. Describing turbulence from the Lagrangian reference frame, albeit challenging, promises to provide much more insight to the physics of the phenomenon (Toschi and Bodenschatz 2009). Lagrangian tracking can be used to measure maps of material acceleration, from which instantaneous pressure fields can be calculated (van Gent et al. 2017). Instantaneous pressure maps are valuable in a wide range of aerodynamics and fluid dynamics applications.

Regarding measurement of material acceleration, van Oudheusden (2013) states: “In this case the accuracy is limited by the reliability with which the fluid path can be reconstructed, based on the sequential velocity fields, i.e. the extent to which it is representative of the path of a true material particle. Especially the curvature of the trajectory may be a limiting factor.” This paper describes a novel particle tracking approach that is capable of directly measuring trajectories with sharp curvature, abrupt changes in direction, and strong recirculation. Reconstruction of hypothetical trajectories from sequential PIV velocity fields is not required. Trajectories are measured with spatial resolution on the order of a pixel. To illustrate this measurement technique, and the results that can be derived from it, experiments were conducted to track 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), Gopalan and 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 Gopalan and 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.

Currently, most tracking algorithms are based on the three-frame or four-frame extrapolation-search methods, as discussed by Gopalan and Shaffer (2012, 2013), Gopalan et al. (2016), Dracos (1996), Willneff and Gruen (2002), Willneff (2003), Ouellette (2006), Ouellette et al. (2006), Bodenschatz et al. (2007), Cierpka et al. (2013), and Dabiri and Pecora (2020). A limited search area is created around the predicted particle position in a future frame. Conflicts (more than one particle in a search area) are handled by applying additional constraints. For example, Dracos (1996), Willneff and Gruen (2002), Willneff (2003), and Lüthi et al. (2005) applied the additional constraints that the particle in a search area with the smallest acceleration is the most probable match.

2 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_{\mathrm{search}}\), and the time between camera frames, \(\varDelta t\). The radius of the NNSA is iteratively increased from a minimum value close to zero. Starting \(V_{\mathrm{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_{\mathrm{min}}^{\mathrm{acc}}\) and \(R_{\mathrm{max}}^{\mathrm{acc}}\), represent a minimum and a maximum translational acceleration, while the angle, \(\phi _{\mathrm{acc}}\), represents an angular acceleration. The search parameters, \(R_{\mathrm{min}}^{\mathrm{acc}}\), \(R_{\mathrm{max}}^{\mathrm{acc}}\), and \(\phi _{\mathrm{acc}}\), are also iteratively varied from a small value to a maximum value. The value of \(R_{\mathrm{min}}^{\mathrm{acc}}\) can be varied in the range of \(R_{1-2}< R_{\mathrm{min}}^{\mathrm{acc}} < 2R_{1-2}\). The value of \(R_{\mathrm{max}}^{\mathrm{acc}}\) can be varied in the range of \(2R_{1-2}< R_{\mathrm{min}}^{\mathrm{acc}} < 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.

Fig. 1
figure 1

Extrapolation and search areas used by Ramer and Shaffer (1992)

Full size image
Fig. 2
figure 2

Nearest Neighbor Search Area (NNSA) at \(t_{i+1}\)

Full size image
Fig. 3
figure 3

Descendant Search Areas (DSA)

Full size image
Fig. 4
figure 4

Candidate tree of particles found in search areas initiated from particle \(P_{11}\)

Full size image

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 30\(\mu\)m 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 \(\mu\)m FCC (fluid cracking catalyst) particles in a 18m tall, 30-cm-diameter riser of a circulating fluidized bed (Shaffer et al. 1992; Singh et al. 1993; Gopalan and 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.

Gopalan and 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.

Fig. 5
figure 5

Schlieren snapshot of submerged jet at \(\hbox {Re}=1000\)

Full size image
Fig. 6
figure 6

Trajectories found in 1000 frames for \(\hbox {Re}=1000\). 56,400 trajectories are shown

Full size image
Fig. 7
figure 7

Close-up view of trajectories in the transition region at \(\hbox {Re}=1000\)

Full size image
Fig. 8
figure 8

\(\hbox {Re}=1000\): a Trajectory 32,696, b speed and c acceleration along the trajectory

Full size image

3 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\times W \times H)\) of 240 cm \(\times\) 120 cm \(\times\) 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 45\(\mu\)m 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.

Fig. 9
figure 9

Schlieren snapshot of submerged jet at \(\hbox {Re}=3000\)

Full size image
Fig. 10
figure 10

Overlay of 6704 trajectories for 500 frames at \(\hbox {Re}=3000\)

Full size image
Fig. 11
figure 11

\(\hbox {Re}=3000\): a Trajectory 6672, b speed and c acceleration along the trajectory

Full size image
Fig. 12
figure 12

Schlieren snapshot of submerged jet at \(\hbox {Re}=12{,}700\)

Full size image

Four flow conditions were studied, with Reynolds numbers of 1000; 3000; 12,700; and 25,200. The Reynolds number, \(\hbox {Re}=4Q/\pi \nu D\), is based on the flow discharge rate Q and the kinematic viscosity \(\nu =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\(^\circ\) off the normal to the PTV laser sheet in order to allow clear 90\(^\circ\) 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 (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 \(\hbox {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-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 Gopalan and 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 close-up view of the trajectories in the region undergoing a transition-to-turbulence at \(\hbox {Re}=1000\). In the shear layers, trajectories undergo small-scale recirculation and form cycloids of varying degrees including flow reversal as evidenced by the abrupt changes in direction. Figure 8a shows an example trajectory, recirculating in the shear layer at \(\hbox {Re}=1000\). An \(8 \times 8\) pixels window is shown in yellow in Fig. 8a. PIV with interrogation window sizes of, for example, \(16 \times 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 \(\mathbf{e}_t\): \(\mathbf{u}=u\mathbf{e}_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)

$$\begin{aligned} {D\mathbf{u} \over Dt} \approx {\mathbf{u}(s_{t_o} +\varDelta t) - \mathbf{u}(s_{t_o}) \over \varDelta t } = {Du \over Dt}{} \mathbf{e}_t \end{aligned}$$
(1)

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 \(\hbox {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 \(\hbox {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 \(\hbox {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 calculated and are shown in 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 \(\hbox {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.

Fig. 13
figure 13

Overlay of 9884 trajectories recognized in 500 frames at \(\hbox {Re}=12,700\)

Full size image
Fig. 14
figure 14

\(\hbox {Re}=12{,}700\): a Trajectory 35,905, b speed and c acceleration along the trajectory

Full size image
Fig. 15
figure 15

Schlieren snapshot of submerged jet at \(\hbox {Re}=25{,}200\)

Full size image

Figure 18 shows a map of time-averaged magnitude of acceleration along trajectories at \(\hbox {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)

$$\begin{aligned} \nabla p= -\rho {D\mathbf{u} \over Dt} + \mu \nabla ^2\mathbf{u}. \end{aligned}$$
(2)

For incompressible flows of high Reynolds numbers, away from boundaries, the material acceleration is dominant and the viscous term can be neglected, so

$$\begin{aligned} \nabla p= -\rho {D\mathbf{u} \over Dt} . \end{aligned}$$
(3)
Fig. 16
figure 16

Trajectories overlain for 1000 frames at \(\hbox {Re}=25{,}200\)

Full size image
Fig. 17
figure 17

\(\hbox {Re}=25{,}200\): a Trajectory 477, b speed and c acceleration along the trajectory

Full size image
Fig. 18
figure 18

Map of time-averaged magnitude of material acceleration at \(\hbox {Re}=12{,}700\)

Full size image

4 Closing remarks

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.