Use of Particle Tracking Velocimetry for Measurements of Granular Flows: Review and Application

Granular flows are commonly encountered in the materials, metallurgical, cement, fertilizer, chemical, coal, food, and pharmaceuticals industries. Examples of granular flows in geophysical contexts include formation and dispersion of sand dunes in deserts, sediment transport in oceans and rivers, debris flows, and snow avalanches. In general, granular flows involve one or more of three broad regimes: (1) static, (2) slowly deforming or quasi-static flows, and (3) rapid flows. Discussions of these and additional regimes of granular flows have been presented by Bagnold (1954), Savage (1984, 1993), Campbell (1990, 2006), Iverson and Vallance (2001), Goldhirsch (2003), and Hutter (2005) among others. Over the last half-century or so, there have been numerous analytical, experimental, and numerical studies of granular flows. Typical examples of such studies include the works of Bagnold (1954), Savage and Sayed (1984), Lun and Savage (1987), Johnson and Jackson (1987), Jenkins (1987), Sinclair and Jackson (1989), Kobayashi et al. (1989), Campbell (1990), Jaeger and Nagel (1992), Ohyama et al. (1993), Okamoto et al. (1995), Savage (1998), Sela and Goldhirsch (1998), Oda and Iwashita (1999), Ottino and Khakhar (2000), Chou (2000), Komatsu et al. (2001), Bonamy et al. (2002), Jop et al. (2005), and Campbell (2006). However, mathematical models of these flows are not yet fully established. Furthermore, a single (unified) model of all of the aforementioned regimes of granular flow is not available, except for some ad hoc efforts aimed at patching together available models of the individual regimes. Efforts to formulate such models are continuing, and there is a need Use of Particle Tracking Velocimetry for Measurements of Granular Flows: Review and Application ─Particle Tracking Velocimetry for Granular Flow Measurements─


INTRODUCTION
Granular flows are commonly encountered in the materials, metallurgical, cement, fertilizer, chemical, coal, food, and pharmaceuticals industries. Examples of granular flows in geophysical contexts include formation and dispersion of sand dunes in deserts, sediment transport in oceans and rivers, debris flows, and snow avalanches. In general, granular flows involve one or more of three broad regimes: (1) static, (2) slowly deforming or quasi-static flows, and (3) rapid flows. Discussions of these and additional regimes of granular flows have been presented by Bagnold (1954), Savage (1984Savage ( , 1993, Campbell (1990Campbell ( , 2006, Iverson and Vallance (2001), Goldhirsch (2003), and Hutter (2005) among others.
Over the last half-century or so, there have been numerous analytical, experimental, and numerical studies of granular flows. Typical examples of such studies include the works of Bagnold (1954), Savage and Sayed (1984), Lun and Savage (1987), Johnson and Jackson (1987), Jenkins (1987), Sinclair and Jackson (1989), Kobayashi et al. (1989), Campbell (1990), Jaeger and Nagel (1992), Ohyama et al. (1993), Okamoto et al. (1995), Savage (1998), Sela and Goldhirsch (1998), Oda and Iwashita (1999), Ottino and Khakhar (2000), Chou (2000), Komatsu et al. (2001), Bonamy et al. (2002), Jop et al. (2005), and Campbell (2006). However, mathematical models of these flows are not yet fully established. Furthermore, a single (unified) model of all of the aforementioned regimes of granular flow is not available, except for some ad hoc efforts aimed at patching together available models of the individual regimes. Efforts to formulate such models are continuing, and there is a need for experimental data that can be used for the verification of such models. The work presented in this paper is a part of a recently completed investigation which was undertaken to obtain experimental data for a flow that spanned the three above-mentioned flow regimes [Jesuthasan (2005)].
The use of particle tracking velocimetry for measurements of granular flows is discussed in this paper. A brief review of publications that deal with particle image velocimetry (PIV) and particle tracking velocimetr y (PTV) is presented first. It is followed by a brief review of papers that discuss the applications of PIV and PTV to granular flows. Next, the application of PTV to granular flows is demonstrated in the context of a recently completed experimental investigation of free-surface flows of almost spherical, slightly polydisperse, ceramic particles immersed in air [Jesuthasan (2005)]. Flows of this granular material down the upper inclined surface of a wedge-shaped static pile of the same material, formed naturally and contained in a narrow channel between two parallel vertical glass plates, as shown photographically in Fig. 1, are considered here. These free-surface granular flows are relatively simple and affordable to set up in a university research context, but they encompass all three of the aforementioned regimes. PTV allows non-invasive measurements of such flows, and the results could be used to aid fundamental works aimed at the formulation, testing, and refinement of a unified theor y (one that covers the whole range of static, slowly deforming, and rapid flow regimes), at least for dry cohesionless materials. This was the primary motivation for the work on which this paper is based [Jesuthasan (2005)]. Some sample results obtained from PTV measurements of these flows in the statistically steady and fully developed region are presented in this paper.

REVIEW OF PIV AND PTV TECHNIQUES
PIV and PTV techniques have been used for measurements of single-phase fluid flows for the last 15-20 years. Reviews of PIV in this context have been presented by Adrian (1991), Willert and Gharib (1991), and Raffel et al. (1998). PIV is a technique where the whole flow field is visualized (measured) at a multitude of points by seeding it with suitable tracer particles, illuminating these particles with a sheet of laser light that is pulsed, obtaining images of these particles on photographic film or a video array detector, transferring these images to a computer, and using suitable computational algorithms and procedures to deduce the velocity field and related information. If the particle images are digitally acquired and stored, using a high-speed digital camera, rather than recorded on conventional photographic film, then this technique is sometimes referred to as digital particle image velocimetr y (DPIV) [Willert and Gharib (1991)]. Adrian (1991) classifies low-image density PIV as PTV. In PTV, images of the flow obtained with multiple exposures are used to track particle displacements that are short compared to the mean spacing between individual par ticles tracks; fur thermore, since the number of particle displacements per unit area of the image is relatively small, it is feasible to track individual particles. As obser ved by Adrian (1991), PTV has its roots in flow-visualization techniques such as particle-streak photography and stroboscopic photography [Van Dyke (1982)]. In modern PTV, the quantitative results are obtained using computerized analyses of the seed-particle images. PTV is effective for investigations of flows with large velocity gradients, a condition that is known to pose problems for conventional PIV methods [Huang et al. (1993a and1993b)]. Two well-established PTV algorithms are the four-frame in-line tracking method [Kobayashi et al. (1989); Hassan and Canaan (1991)], that uses four consecutive images of the tracked par ticles, and the binary image cross-correlation method [Uemura et al. (1989)], which makes use of only two such consecutive images. In the four-frame tracking algorithm, the movement of the tracer particles is detected frame-by-frame, while evaluating the geometrical consistency of every particle path. Typically, an iterative scheme is used to select the best-match par ticle trajector y, involving extrapolation of the particle displacement and searching for the nearest neighbor. On the other hand, the binary image crosscorrelation technique is viewed as a variation of the standard cross-correlation PIV, in which the correlation functions are computed for each interrogation window centered on the first frame of particles using an adaptive shifting scheme.
Additional two-frame algorithms that employ a particle-cluster matching concept have been proposed [Okamoto et al. (1995), Ishikawa et al. (1997)]. In this method, particles from the first and second frame are conceived to form a cluster with their respective neighboring particles, and the selection of the best match particles is conducted on the basis of a deformation index defined for the relationship between the clusters in the two frames. Another two-frame technique is based on the use of a suitable cost function in particle tracking algorithms: An example is the algorithm proposed by Ohyama et al. (1993), in which the best-match particle pairs are determined by using a fitness function that minimizes the total sum of the squares of particle displacements. The neuralnetwork PTV proposed by Labonté (1999) can also be classified in this category of algorithms. Although time consuming, this algorithm seems more efficient (reliable) in identifying and tackling unpaired particles between the frames.

REVIEW OF APPLICATIONS OF PIV AND PTV TO GRANULAR FLOWS
Over the last 15 years, continual improvements of computer hardware and software have greatly facilitated the use of PIV and PTV for measurements of granular flows. Today, digital PIV and PTV systems and procedures, and modern high-speed digital cameras, enable the speedy acquisition and processing of massive amounts of data with excellent spatial resolution.
As was mentioned earlier, PTV is better suited than PIV for handling flows with steep velocity gradients, such as those commonly encountered in granular flows. Nevertheless, PIV has been used successfully for measurements of granular flows. Examples of such efforts include the works of Warr et al. (1994), Lueptow et al (2000), Hanes and Walton (2000), and Tischer et al. (2001). Jain et al. (2002) have discussed the difficulties associated with the application of PIV to thin granular shear layers. In their rotating-drum experiments, the granular shear layers adjacent to the free surface were only 5-10 particle diameters deep, leading to relatively steep velocity gradients. Thus, they resorted to the use of a hybrid PTV/PIV technique proposed by Cowen and Monsmith (1997). Larcher (2002) has performed experiments on suspensions of solid particle in liquids flowing down inclined chutes. He measured mean streamwise velocities of the particles, granular temperature, and solids fraction profiles over the depth of the flow. Larcher also described his attempts to use PIV techniques to determine the mean velocity profiles. However, at high shear rates, problems of the kind discussed by Jain et al. (2002) were encountered. Larcher also used a particle-tracking technique based on Voronoï diagrams around the par ticle centers in paired frames: The particles in one frame and the next were paired based on a goodness of match between the shapes of the Voronoï cells. The particle velocities were found by dividing the corresponding particle displacements by the time inter val between the frames. Larcher (2002) found this Voronoï approach to work effectively at shear rates higher than those possible with the standard PIV method. A detailed description of particle tracking based on Voronoï imaging methods is presented in Capart et al. (2002). An extension of these ideas to a three-dimensional Voronoï imaging method is given in Spinewine et al. (2003).
Over the last ten years or so, there have been several experimental investigations of free-surface granular flows in which optical techniques have been used to obtain measurements of particle velocities, and, in some cases, also their fluctuations. Two different experimental set-ups have been primarily used in these studies. In one, the granular material flows down a static pile (or heap) of the same material, usually contained between two vertical walls. The granular material is released from a hopper on top of the pile, and its flow rate out of the hopper is controlled: examples include the works of Komatsu et al. (2001) and Khakhar et al. (2001), and also the recent work of Jop et al. (2005). In the second set-up, a cylindrical tumbler, or drum, whose axis is horizontal (with respect to the gravitational acceleration vector), is halffilled with the granular material and rotated about its axis at a controlled angular speed to create the freesurface granular flow: examples include the works of Jain et al. (2002) and Bonamy et al. (2002). The rotating-tumbler arrangement is very convenient experimentally, since with a finite mass of particles, the steady-state free-surface flow can be continued indefinitely for all practical purposes. However, in the work of Jesuthasan (2005), on which this paper is based, the granular flow down a static pile was the chosen configuration for the experimental set-up (Fig. 1). The motivation for this choice is that under statistically steady conditions, the free-sur face granular flow down a static pile can achieve a thickness of the flowing layer that is independent of the distance from the hopper (fully developed condition), whereas in the rotating-tumbler set-up, this thickness varies continually along the direction of the flow. Furthermore, PTV was selected for measurements of the granular flows of interest. Full details of this technique and the related procedures are available in Jesuthasan (2005).

APPLICATION OF PTV TO A FREE-SUR-FACE GRANULAR FLOW
PTV measurements of a free-surface granular flow down the upper inclined surface of a wedge-shaped static pile of the same material, contained in a narrow channel between two parallel vertical glass plates, were undertaken in this work. A photographic illustration of this flow is given in Fig. 1. In the design of the experimental set-up, some ideas contained in the works of Komatsu et al. (2001) and Khakhar et al. (2001) were borrowed and extended. An overview of the experimental apparatus and techniques, a synopsis of the PTV procedures, and some sample results from this investigation are presented in this section.
For full details of this research, the readers are referred to the work of Jesuthasan (2005).

Over view of the Experimental Apparatus
and PTV set-up A photograph of the experimental apparatus is provided in Fig. 2. It consists of the following main elements located in sequence along the flow path of the granular material: (i) a reservoir hopper that was designed to safely hold a sufficiently large amount (up to 20 L) of the granular material; (ii) a heightadjustable spout with a double-slider gate mechanism (schematically illustrated in Fig. 3) which was used to control the flow rate of the granular material; (iii) a test section in which the free-surface flow was set up and measured (see 38.1 mm, and 50.8 mm. For each of these separation distances, the following three values (nominal) of the mass flow rate per unit width were considered: 0.811 kg/s-m, 1.847 kg/s-m, and 3.332 kg/s-m. For these ranges of dimensional parameters, statistically-steady fully developed flows were established inside the test section in a reliable and repeatable manner [Jesuthasan (2005)].
A high-speed camera system (Redlake Motion-Pro 10000) was used to acquire digital images of the granular flows of interest in the fully developed region. A photograph of this high-speed camera and the system that was used to mount and traverse it is shown in Fig. 4. Illumination of the front surface of the granular flows of interest was done using two 250 W halogen photoflood variable-intensity lights (Lowel Pro-Light). Image processing and a commercial PTV software package (DiaTrackPro 2.3) were used to obtain the particle trajectories, and special programs were written in Matlab to obtain the corresponding instantaneous and ensemble-averaged velocity distributions and related results.

Synopsis of the PTV Procedures
The high-speed digital camera, after it was properly positioned and aligned to capture images of the granular flows, was first used to acquire images of a plumb ball and a ruler placed on the inside surface of the front glass plate of the test section. The images of the plumb ball were used to establish the orientation of the camera with respect to the horizontal, and then the static angle of repose of the granular particles was determined. The images of the ruler were employed to determine a conversion factor for relating data in pixels counts to physical length dimensions. These analyses were done using an image-processing computer program called ImageTool, which is available on the Internet at http://ddsdx.uthscsa.edu/ dig/itdesc.html. This program allows angular measurements with a precision of about 0.05 o .
As was stated above, a commercial PTV software package (DiaTrackPro 2.3, developed and distributed by Semasopht, Swizterland) was purchased and used in this work. The initial step in the PTV analysis using this software package is to load the first image frame (of a sequence of such frames) in the appli-  cation window. The image is then convoluted by a Gaussian filter of variable strength such that centers and boundaries of the particles are highlighted [Willert and Gharib (1991)]. The locations of the center of mass of the particles in the image are then identified based on local brightness maxima. Incorrectly identified particles are selected and eliminated by imposing suitable image thresholds [Jain et al. (2002); Larcher (2002)]. Once this is done for the first image frame, the aforementioned processing steps are automatically applied to all successive image frames in the sequence. At this point, the tracking algorithm is initiated and particles in each pair of successive frames are sequentially matched [Adrian (1991); Cowen and Monismith (1997); Larcher (2002)] to produce a series of particle displacement vectors. This sequence of steps is illustrated by the images provided in Fig. 5.
The PTV data processing procedures are summarized in the flow diagram given in Fig. 6. First, the high-speed digital camera is used to acquire sequential images of the granular flow. Then the acquired raw digital images are rotated and cropped, using programs written in Matlab, to aid the PTV analyses using DiaTrackPro 2.3. Finally, specific data processing for the granular flows of interest is done, again using programs written in Matlab.

Data Processing
Upon successful completion of PTV analysis via DiaTrackPro 2.3, the resulting particle displacement data were used to compute the corresponding instantaneous velocities. The following simple approximation was used in this step: where ri (m 1) (xi (m) , yi (m) ) is the position of the i th particle in the measurement plane at time t (m) , Vi (m) (ui (m) , vi (m) ) is the corresponding in-plane velocity, and x and y are oriented in the directions parallel and perpendicular to the mean flow direction, respectively. The particle displacements yielded by the PTV analysis are in pixels counts: The conversion factor mentioned in the previous subsection was used to relate the image pixel counts to physical lengths.
Once the instantaneous velocities were computed and stored for each image pair in the sequence, a program written in Matlab was used to compute the following granular flow field properties: mean velocities, mean-square velocity fluctuations, granular temperature, and the Savage-Jeffrey parameter.
First, the flow region of interest was subdivided (d) final image after particle selection; (e) superimposed particle displacement vectors after PTV analysis. Fig. 6 Flow diagram of the PTV and data processing procedures.
into non-overlapping horizontal slices (or bins) of thickness, δ, equal to one particle diameter, as shown schematically in Fig. 7. The coordinates of the centers of the measured particle displacement vectors were used to distribute them into these bins. They were then processed to obtained ensemble-averaged results for each bin. The ensemble-averaged streamwise velocity, k, and transverse velocity, v k, within the k th bin, at its centroid, Yk, were calculated as follows: In this equation, Nk is the total number of displacement vectors in the k th bin. Second-order statistics such as the square of the velocity fluctuations, the granular temperature, and the Savage-Jeffrey parameter were then evaluated. Within each bin, the streamwise and transverse fluctuating velocities were defined as follows: The means of the squares of the streamwise and transverse velocity fluctuations were then calculated as follows for each bin: In this work, the z-direction fluctuating velocity component was assumed to be comparable in magnitude to the y-direction fluctuating velocity component, and the granular temperature was calculated as follows: The Savage-Jeffrey parameter, R, was calculated as shown below: In this equation, d is the mean particle diameter. The streamwise velocity gradient in Eq. (6) was determined using cubic-spline interpolation of the binned ensemble-averaged velocity data.

Sources of Error and Benchmarking
The PTV operations are prone to errors. Particle identification from one frame to the next could be erroneous. Limited accuracy in particle positions induces physically unrealistic disturbances to the velocity data. Particle mismatch between successive frames due to incorrect tracking leads to so-called "outliers". As a result, erroneous and/or spurious velocity vectors could contaminate the field data.
Another set of errors is associated with the experimental set-up. The high-speed digital camera and related optics, as well as the lighting system, operate within their own levels of uncertainties. The camera lens and zoom level can contribute to image distortions, compromising the accuracy of the acquired images. An appropriate choice of frame speeds in image recording is also important. Long time steps between successive images damps high-frequency fluctuations. On the other hand, short time steps lead to small displacements and large pixel errors.
For detailed and authoritative discussions of errors involved in PTV measurements, the reader is referred to the works of Adrian (1991), Willert and Gharib (1991), and Cowen and Monismith (1997).
Benchmarking of the PTV software, experimental set-up, and procedures is needed in order to evaluate the effects of the aforementioned errors on the measurements. In this work, the PTV software was validated via standard digital images obtained from the Japanese Standard Image Project, following the recommendations of Okamoto et al. (2000). These standard image sets are offered free of charge on the website of the Visualization Society of Japan (http:// www.vsj.or.jp/piv). Additional details of these and related benchmarking procedures are available in the work of Jesuthasan (2005).

Sample Results
In this subsection, the intention is to demonstrate the type of results that can be obtained using PTV measurements of granular flows. Thus, with reference to the test section shown in Fig. 1, sample results are presented only for one normal separation distance, W 25.4 mm, between the ver tical glass plates, and the following three different mass flow rates per unit width: m /W 0.775 kg/s . m, 1.799 kg/ s . m, and 3.288 kg/s . m. For additional results and related discussions, the readers are referred to the work of Jesuthasan (2005). The static angle of repose for the separation distance of W 25.4mm was determined to be θrep 20.55 o (standard deviation of 0.44 o ).
The results presented in this subsection are based on measurements acquired and ensemble-averaged over a sampling period of 10 s. It was established by Jesuthasan (2005) that this sampling period is sufficient to provide statistically repeatable results for four different frame speeds: 250 fps, 333 fps, 500 fps, and 1000 fps. Furthermore, Jesuthasan (2005) determined that the optimal frame speeds for the cases considered here are the following: 333 fps for m /W 0.775 kg/s . m; 333 fps for m /W 1.799 kg/s . m; and 500 fps for m /W 3.288 kg/s . m.
In the presentation of the sample results, X denotes the streamwise coordinate; and Y is the coordinate normal to X, directed into the flowing layer towards the underlying static pile, and zero at the free surface. The ensemble-averaged free surface velocity is denoted by Usurf, and the angle of inclination of the free surface of the flowing granular material is denoted by θfs.
The X-component of the velocity of the flowing layer has its maximum values at the free surface (Usurf) and it decays to zero with increasing Y. The exact Y location where this velocity component reaches zero is difficult to determine from the PTV measurements. It was assumed here (somewhat arbitrarily) that the bottom of the flowing layer of granular material (Y H) corresponds to the position where U/Usurf 0.001.
Distributions of the mean (ensemble-averaged) streamwise and transverse velocity components (U and V) with Y are given in Figs. 8a and 8b, respectively. As was discussed earlier, U reaches its maximum value at the surface (Y 0) and decays down to zero with increasing Y. The transverse mean velocity component, V, is zero throughout, for all practical purposes, indicating that the ensemble-averaged flow of the granular material is parallel to the X direction. Values of the tangent of the angle of inclination of the free surface, θfs, the depth of the flowing layer of granular material, H, as defined in Equation (7), and the maximum velocity, Umax Usurf, are given in Table 1. Distributions of the mean-square velocity fluctuations in the X and Y directions ( (u ) 2 and (v ) 2 ) and the granular temperature, T, as defined in Equation (5), are given in Figs. 9a, 9b, and 9c, respectively: Each of these ensemble-average quantities achieves its maximum value at the free surface and decays to zero with increasing Y; at corresponding locations, the values of (v ) 2 are lower than (u ) 2 by a factor of two to three; and all of these quantities increase with increasing values of m /W.
Distributions of the Savage-Jeffrey parameter, R, are presented in Fig. 10   mensionless parameter requires inputs of the gradient of the streamwise velocity, UX, Y, and the squareroot of the granular temperature, T 1/2 , as shown in Equation (6). The values of UX, Y, were calculated using cubic-spline fits to the nodal values of streamwise velocity, Uk, available at the centroids of the bins, Yk .
In the vicinity of the free-surface and also as the static pile is approached, the errors in these calculated velocity gradients could be substantial. Furthermore, as the static pile is approached, both the velocity gradients and the granular temperature go to zero, thus the calculation of R using Equation (6)   sented in Fig. 10: The variation of R with Y is linear at the lower mass flow rates per unit width; and at the highest mass flow rate per unit width, R assumes an essentially constant value. These results seem to indicate that for the free-surface granular flows considered in this work, when the wall friction effects are significant (with respect to the shear and normal stresses in the interior of the granular flow) R vs. Y is almost linear; and R is essentially constant for the other cases.

CONCLUDING REMARKS
A brief review of publications related to PIV and PTV techniques was presented first in this paper. That was followed by a brief review of papers on the applications of these techniques to granular flows. The application of PTV to a demonstration problem involving a free surface granular flow was presented next, along with some sample results. The aforementioned reviews and results show that PTV is a convenient and reliable technique for non-invasive measurements of granular flows, and these measurements can be processed to obtain results that could be useful for checking and refining mathematical models of such flows. Mean-square fluctuating transverse veloc-ity component ( (v ) 2 ) 1/2 Root-mean-square of the transverse fluctuating velocity V Mean velocity component in the Y direction (transverse) W Separation distance between the vertical glass plates (gap width) X Streamwise coordinate oriented along mean flow direction (same as x) Y Coordinate normal to X, directed into the flowing layer (same as y) Z Transverse coordinate to the flow, directed from the front glass wall to the back glass wall (same as z) θfs Free-surface inclination angle with the horizontal θrep Static angle of repose of the granular material

ACKNOWLEDGMENTS
The authors gratefully acknowledge sponsorship of this work by Natural Resources Canada (NRCan), under the auspices of the Climate Change Technology Innovation Initiative, through a project on granular multiphase flow directed by the Canadian Hydraulics Centre (CHC) of the National Research Council (NRC). Thanks are also due to the Natural Sciences and Engineering Research Council of Canada (NSERC) for a scholarship to the first author and research grants to the second and third authors.