利用三維顆粒顯像儀分析水與土運動之研究(3/3)

I 利用三維顆粒顯像儀分析水與土運動之研究(3/3) 摘要 本研究旨在發展一套實驗室之試驗設備，以探討土石流及因地震引發之崩塌 滑動。這些現象存在於台灣本島地區，因台灣地勢處於陡坡又強烈雨型與地震頻 繁之地理位置。這些現象有一些極重要的力學機制，即在飽和顆粒力學，如何由 固態過渡成為液態。由於因為實驗儀器及設備的普遍不足，使得在物理觀念的瞭 解及數值模擬的進行，受到極大的阻礙。 本研究即針對上述的研究動機提出下列的研究方法：（1）設計及建造一套土 石流循環水槽；（2）建立反射因子映象，以探討土石流內部水與沙顆粒之間的變 形及混和；（3）建立三維數值顆粒顯像流速儀，以測量水與土沙顆粒之運動狀態。 透過這些實驗設備及數據的量取，將有助於土石流及土石滑動之微力學之物理機 制及土石流基本流場之瞭解。 我們將應用這套系統（此系統將安置於台大水工所新設置之現代流體力學實 驗室內），來研究下列三大重要主題：（1）因地震引發之邊坡破壞及土石流；（2） 土石流之基本物理機制：（3）土壤液化與土石流之不穩定之現象。吾人相信這套 最新的設備，將可以達到上述三項目標。同時，透過高解析度的量測系統，將可 量到土石流顆粒之間的內部變形，並能將這些量得之實驗數據，用作理論及數值 計算之回饋，也可以將這些數據提供給其他同好研究工作者做其他研究之用。相 信透過本研究計畫，吾人對台灣地區土石流發生之原因及防治對策，將可更上一 層樓。 關鍵詞 土石流；因地震引發之崩塌、滑動；土壤液化；固體－流體之微力學；反射因子 映象；三維數值顆粒顯像流速儀。

II

3D Particle Imaging Study of Liquid-Granular Flows and Flowslides

English summary

The present report documents a three-year research project aimed at developing novel experimental methods for the laboratory study of debris flows and flowslides. Transparent access to the micromechanics of liquid-granular flows is sought through a combination of three components: 1) liquid and granular media with special optical properties allowing unhindered 3D imaging; 2) a tilting flume equipped with two recirculation circuits for the liquid and granular phases; 3) digital imaging algorithms allowing acquisition and analysis of the 3D particle flow field inside the bulk. Key results from the project include:1) the successful preparation of transparent liquid-granular mixtures with small visible cores at the centre of the granular particles, achieved using Refractive Index Matching (RIM) and laser marking; 2) the demonstration that stereo imaging techniques can reliably track the 3D positions and motions of large sets of such marked particles, using high-resolution video and Particle Tracking Velocimetry (PTV) algorithms; 3) the design, implementation and testing of a novel double recirculation flume in which rapid downslope surges of granular-liquid mixtures can be made stationary and characterised in detail.

III

目錄：

1. Background and objectives 1

2. New experimental materials and apparatus 2

2.1. Refractive index matching and laser marking 2

2.2. Fluidisation cell apparatus 6

2.3. Double recirculation flume for the study of downslope granular-liquid flows 8

3. New imaging measurement techniques 13

3.1. Statistical analysis of Particle Tracking Velocimetry (PTV) data 13

3.2. Three-dimensional Particle Tracking Velocimetry (PTV) measurements 14

Goals attained 19

1

Project contents 1. Background and objectives

Subject to heavy rainfall and frequent earthquakes, the upland areas of Taiwan are exposed to high risks of debris surges and flowslides. A variety of applied research efforts are currently directed at mitigating the hazards due to these flows. Many fundamental aspects of their behaviour, however, remain poorly understood. One major reason lies in the difficulty of probing experimentally the flow behaviour of high-concentration granular-liquid mixtures. Intrusive instruments perturb or are perturbed by the flow, while grain occlusion precludes optical access to the flow interior.

The present research seeks to address this challenge by developing novel experimental methods. The methods sought feature three components: 1) special liquid-granular materials allowing optical access and unhindered 3D imaging; 2) new laboratory apparatus allowing rapid phenomena to be made stationary and probed at leisure; 3) post-processing methods allowing insights to be gained from the large sets of data typically acquired by Particle Tracking Velocimetry (PTV) methods; 4) novel imaging algorithms allowing the capture of 3D positions and motions of dense sets of particles in liquid-granular dispersions.

Documented in the present report, the techniques developed to meet these requirements rely on the following principles: 1) the use of Refractive Index Matching to achieve transparency of liquid-granular mixtures, and the laser marking of particle cores to allow their tracking; 2) the assembly of a fluidization cell for pilot tests, and the development of a novel tilting flume equipped with two recirculation circuits for the liquid and granular phases; 3) statistical flow field analysis methods, highlighting correlations of granular motions in space and time; 4) robust stereo imaging methods based on the simultaneous matching of particle positions between views and tracking of particle positions in time.

In the next sections, we first describe the “harware” part of the research project: the development of new experimental materials and laboratory apparatus. We then proceed to describe the “software” part: the development of novel flow field analysis and stereo imaging algorithms.

2

2. New experimental materials and apparatus

The laboratory methods developed involve three components: 1) the preparation of a novel particle-liquid system using Refractive Index Matching (RIM) and laser marking methods; 2) the assembly of a fluidization cell apparatus for the conduct of pilot tests; 3) the design and development of a novel double recirculation flume to study rapid downslope flows of granular-liquid mixtures.

2.1. Refractive index matching and laser marking

A major objective of the overall research programme was to identify a combination of liquid and solid materials that would allow observation of rapid motions inside a dense 3D liquid-granular mixture. This goal was attained through a combination of refractive index matching (RIM) and laser marking techniques. The combination of these two techniques yields a transparent dispersion of particles marked with a small but highly visible core suitable for identification and tracking.

After consideration of various alternatives and testing on small samples, a solid-liquid combination proposed recently by Haam et al. (2000) was selected as RIM system. Solid particles consist of small 7 mm spheres of poly-methyl metacrylate (PMMA), a transparent material more commonly known as acrylic or perspex. The corresponding index-matched liquid is 1-Methyl-4 - (1-methylethyl) - benzene, an organic liquid traded under the name

para-cymene. It has a low viscosity suitable for high Reynolds number experiments. Acrylic

particles of excellent optical qualities are available commercially in Taiwan, while para-cymene can be imported at reasonable cost. Neither material presents excessive flammability or health risks. The characteristics of both media are summarised in Table 1.

property / material para-cymene PMMA

specific gravity 0.855 1.19

viscosity [ Pa s ] 1.023 × 10−3 _n.a.

refractive index 1.489 1.49

Table 1. Summary of the properties of the selected liquid and solid media.

Thanks to the close match between the refractive indices of these two media, a high degree of transparency is achieved for their mixture, eliminating to a large extent the light ray distorsion occurring across the boundaries of the spheres. This allows unhindered optical penetration even for high concentrations of particles. This is illustrated on Fig. 1. The density ratio between the two phases also falls in a range which is well-suited for laboratory investigations of gravity-driven liquid-granular flows.

3

Fig. 1. Refractive index matching. Left: transparent particles visible in air become invisible when bathed in the index-matched liquid; right: a grid is seen undistorted through the matched liquid-

solid mixture. Faint outlines of the solid spheres can be detected upon closer look at the photo. The refractive index of the liquid varies with temperature. As illustrated in Fig. 2, the matching with the acrylic solid spheres is optimal only within a rather narrow range of temperatures close to 15oC. Since this is lower than typical ambient temperatures, refrigeration equipment must be used to cool the liquid before and/or during the experiments. With adequate temperature control, excellent transparency can be achieved.

Fig. 2. RIM dependence on temperature. Left: the transparency is better at temperature of 15oC; right: the outline of PMMA spheres is stronger at temperature of 25oC.

Transparency of course is not an end in itself, and our ultimate purpose is to be able to track individual particles located inside the flowing bulk. To achieve this, particles must be made visible again. Previous RIM researchers (e.g. Cui and Adrian, 1997; Haam et al. 2000) have all solved this problem by coating the outside of a small subsample of particles with a

4

chemical dye or a metal plating. As a result, only this small subset (typically 1-5% of all the particles) can be seen and tracked on digital images, leading to sparse velocity measurements. This limits the resolution of the measured flow field and precludes capture of the interacting motions of neighboring grains.

After exploring some coating and inclusion options, we were able to solve the problem in a novel way thanks to a recently developed technology in use in Taiwan's crystal ornament industry. The technique called internal laser etching is used to carve 3D patterns of visible dots inside transparent blocks of lead crystal. Thanks to the cooperation of a local crystal company, we were able to apply this technique to acrylic spheres, and mark a small visible core at the centre of the transparent particles (Fig. 3). The key advantage of this technique is that a particle can be seen without occluding other particles along a given line of sight. As a result, rather than a small subset, all particles can be marked and used for velocity tracking.

Fig. 3. Left photo: visible core laser-etched at the centre of an acrylic sphere. Right photo: flasks of unmarked (left) and laser-marked (right) acrylic spheres bathed in para-cymene.

To process a large number of particles, the marking procedure had to be made efficient. High-precision machining was used to prepare a pair of identical perforated plates each of which can hold ~ 400 particles. While one plate is being loaded with particles, the other is mounted on a computer-controlled x-y-z carriage which sequentially moves particles below the laser beam. In ~ 8 hours of laser bench time, it was possible to etch 20,000 particles in this fashion. The process is illustrated in Figure 4. The laser-marking approach does not seem to have been used before in RIM experiments, and should allow us to obtain high resolution measurements of the motions of interacting grains.

5

Fig. 4. Laser-marking process. Top: sequential laser-marking of the particles placed on the carriage-mounted plate; bottom: close up of the laser lens and carriage during the laser-marking

sequence.

After successful testing of the techniques on small samples liquid and solid materials, a batch of 20,000 acrylic particles and two drums of 200 liters of para-cymene were obtained. Laser marking of the large number of particles was performed, and the materials have been successfully prepared for use in the fluidisation and recirculation tests.

6

2.2. Fluidisation cell apparatus

In a second stage of the project, a fluidisation cell apparatus was assembled to conduct pilot tests of the novel particle-liquid system described above. The aim of these pilot tests was threefold: to demonstrate the feasibility of refractive index matching RIM measurements; to troubleshoot materials and methods to be used later in the large scale flume; to obtain a first set of reference measurements in homogeneous flow conditions.

Figure 5. Fluidisation cell apparatus. From left to right: a) transparent-walled device; b) schematic; c) close-up of the inner cylindrical working section. Temperature control is achieved here with the

cooling coil seen on the right.

Through these pilot tests, the combination para-cymene / acrylic was confirmed as an excellent choice of transparent liquid-granular mixture. Temperature control was achieved by plunging a cooling circuit in the liquid bath (see Fig. 5c) , or by pre-cooling the liquid in a refrigeration unit. The marking of the particle cores using laser etching, described above, was also found to yield excellent results. The best viewing conditions are found to be obtained using oblique backlighting, with stereo views angled at around 10 degrees. Video images obtained in such conditions yield sharply contrasted black cores on a light background. The cores are thus easily identifiable on the digital images, while at the same time being small enough to minimise occlusion effects. An example of a digital frame is shown on Fig. 6. The corresponding tracking measurements are discussed infra in the present report.

7

Figure 6. Sample video image of a 3D fluidisation of particles: the particle cores show up as sharp black speckles on a light background, and are identifiable in both left and right

windows of the stereo image.

The tests have also highlighted various technological issues which require special care. First, wall, sealant, and piping materials must be compatible with the para-cymene liquid. The set of compatible materials was found to include acrylic, nylon, PVC, glass, stainless steel, and resin, but to exclude common sealants and plastics used for water flow experiments. Secondly, conditions of uniform liquid inflow are best approximated using a combination of expansion chamber and filtration layer of heavy particles. These results have helped guide the development of the larger scale flume apparatus discussed in the next section.

8

2.3. Double recirculation flume for the study of downslope granular-liquid flows

The third major component of the project was the design and development of a novel double recirculation flume to study rapid downslope flows of granular-liquid mixtures. The originality of the new apparatus is its ability to independently recirculate the liquid and granular phases involved in the flow. A further constraint is that all materials used for the flume have to be compatible with the special RIM liquid. This aims to insure that the flume will be suitable for the acquisition of 3D measurements using the transparent materials.

Figure 7 shows the conceptual design proposed initially to achieve this purpose, while Figure 8 shows the revised design adopted after a more comprehensive study. Instead of performing the double recirculation using two conveyer belts, the new design uses one conveyer belt to recirculate the solid phase, and an external pumping circuit to recirculate the liquid.

2 3 1 1 3 4 5

Figure 7. Initial design for the proposed double recirculation flume (principle sketch): 1) glass-walled open-channel; 2) internal conveyer belt; 3) external conveyer belt; 4) imaging sensor

9

Figure 8. Revised design incorporating technological and economical constraints, in which the double recirculation is achieved by combining an internal conveyer belt

with an external pumping system.

A first version of the flume apparatus was assembled according to this design. A view of the flume installation is shown on Figure 9. Key components include the adjustable slope, the high-speed conveyer belt for granular recirculation, and the independent pumping circuit for liquid recirculation. All three components of the device are controlled using an electronic panel. Operating parameters can be set precisely and varied over a wide range, in accordance with the design.

10

Figure 9. Installation of the first version of the double recirculation flume.

While the overall design was successfully implemented, difficult mechanical problems were encountered with the conveyer belt. The purpose of this belt is to entrain granular material upstream along the floor of the flume. Driven by a variable speed motor, it forms a continuous loop around the apparatus (see Fig. 9). As this belt is in direct contact with the liquid-granular mixture, its material must be compatible with the para-cymene liquid. This led to the choice of a flexible, stainless steel perforated belt entrained by solid rollers with matched pins. The original design called for a series of rollers of 10 cm diameter placed at the four corners of the flume. To form a continuous loop, the two ends of the belt must be carefully welded together.

The scheme was implemented and found to work fine, but only for a limited number of cycles: repeated bending of the joint as the belt curves around the rollers rapidly leads to fatigue failure (see Fig. 10). Improved endurance is obtained by subjecting the welded joint to an annealing process prior to belt installation, however joint failure nevertheless occurs after a number of cycles of the order of 1,000 to 10,000. This has forced us to reconsider the design of the belt circuit: instead of 4 small rollers at the four corners of the flume, the revised design features two large wheels of 70 cm diameter at both ends of the flume. The enlarged radius is crucial to minimize bending of the belt and was calculated to extend the belt lifetime by at least 2 orders of magnitude.

11

Figure 10. Fatigue failure of the welded joint of the stainless steel belt. The revised design currently being implemented replaces the small rollers at the flume corners by two wheels of large diameter.

Figure 11. Replacement of the small rollers by large wheels to reduce belt fatigue problems. Testing at the workshop.

Figure 12. New version of the double recirculation flume installed in our laboratory. The flume has now undergone a number of tests and is fully operational.

12

The revised double recirculation flume is shown on Figures 11 and 12. Thanks to the wheel replacement, fatigue problems have been solved. The flume was found to perform well for a number of tests, and is now fully operational. A first set of experiments has been carried out, aimed at identifying the flow regimes that can be obtained under different conditions. One flow pattern of particular interest for debris flow studies is shown on Fig. 13. Both surges and uniform flow regimes are expected to be targets of great interest for further studies.

Figure 13. Granular-liquid surges observed in the double recirculation flume. The flume is inclined and the belt operated at large speeds to make the downslope surge stationary with respect to the laboratory frame of reference. Top: overall view of the surge. Bottom: close-up of the front, where

13

3. New imaging measurement techniques

To complement the novel materials and apparatus described in the previous section, a second major part of the work undertaken during the three-year project consists in developing new imaging measurement techniques. These techniques fall in two different categories: 1) statistical methods used to extract information from the very large datasets typically generated by PTV measurements. 2) 3D imaging algorithms aimed at reconstructing the 3D positions and motions of large sets of particles from stereo images.

3.1. Statistical analysis of Particle Tracking Velocimetry (PTV) data

Particle Tracking Velocimetry (PTV) measurements typically yield very large sets of data. To handle such data sets, statistical flow analysis techniques have been developed and refined in the course of the project. While the new experimental devices described above were still under development, the analysis techniques were first applied to footage obtained using existing equipment. This was done in cooperation with two other universities: 1) a detailed analysis of steady uniform dry granular flows was performed in cooperation with the National Central University (Professor H.-T. Chou); 2) a comprehensive analysis of steady uniform granular-liquid flows was conducted in cooperation with the Università degli Studi di Trento (Professor A. Armanini). Two manuscripts based on this work have been submitted for publication and are included in appendix.

Fig. 14. Liquid-granular flow imaged through a flume sidewall. Left: side view image (mean flow is from left to right); right: convected trajectories (i.e. measured grain trajectories from which the

mean flow component has been subtracted).

We have developed statistical tools for two different purposes: first, to isolate physical signals from random experimental errors; and secondly to identify correlated motions in fluctuating flow fields. As an illustration of what statistical techniques can help achieve, Fig. 14 shows a debris flow run imaged in 2D through a flume sidewall, from which detailed particle trajectories were extracted (this run is taken from a previous study done in cooperation with the Università degli Studi di Trento, Italy). A spatially correlated flow pattern could be

14

extracted from these 2D measurements, and is shown on Fig. 15. Coherent motion components can thus be successfully extracted from a seemingly random granular flow field. A variety of other results obtained by statistical post-processing of PTV data are documented in the two manuscripts included in appendix.

Fig. 15. Correlated flow pattern averaged around a single particle, as sampled inside the fluctuating granular ensemble of Fig. 14. Each unit of the spatial grid correspond to one grain diameter. 3.2. Three-dimensional Particle Tracking Velocimetry (PTV) measurements

Three-dimensional particle tracking algorithms were developed and refined in the course of the project, and tested on the fluidisation cell experiments. Image acquisition is performed using a single camera and a system of mirrors. The acquired images correspond to two distinct viewpoints, which have to be calibrated independently. The calibration target used for this purpose is also shown on Fig. 16. The illumination is provided by backlighting, using spots aimed at translucent acrylic board placed behind the fluidisation cell.

15

Once images have been acquired in this fashion, 3D particle positioning and tracking can be performed. A first set of methods for the 3D positioning and tracking of particles viewed on stereoscopic image pairs was developed and tested in cooperation with the Université catholique de Louvain (Belgium) and Università degli Studi di Trento (Italy). The algorithms are described in detail in the Appendix and we will not repeat the description here. Illustrating the kind of results obtained, Fig. 17 shows the fluctuating 3D trajectories of an ensemble of opaque particles viewed through a transparent wall. Due to particle occlusion, here the depth of field is limited to the vicinity of the sidewall.

Fig. 17. Near-wall trajectories of a fluidised dispersion of opaque grain bathed in water: (a) top view; (b) oblique view; (c) frontal view through the sidewall. Measurements obtained using the stereo

matching and tracking methods described in Appendix.

The same methods were applied to the fluidization cell experiments. Thanks to the RIM materials and laser marking technique, of course, particles inside the 3D dispersion can now be tracked, rather than only those close to the sidewall. When only 10 marked particles were introduced in a fluidisation of 1,000 transparent particles, the methods generated good results (see Fig. 18). However, when the full set of 1,000 particles were laser-marked, two problems were encountered. First, the temporal and spatial resolution of the video camera was no longer sufficient to resolve the motions of all the particles. Secondly, the imaging algorithms themselves were not sufficiently robust, yielding poor tracking statistics.

16

Figure 18. 3D particle tracking of 10 marked particles in a fluidization of 1,000 transparent spheres bathed in a liquid of identical refractive index.

To solve these problems, new footage was acquired using a high-resolution video camera on loan from the Taiwan representative of the camera maker. Images acquired using this camera are shown on Fig. 19. To take advantage of this higher resolution and higher frame rate, new algorithms were developed. Instead of positioning particles in 3D first, and then tracking them, the new methods combine the 3D positioning and tracking in a single step. This is aided by a preliminary step, shown on Figure 20, in which particles are first tracked in the 2D image plane.

Figure 19. Sample high-resolution video image: the particle cores show up as sharp black speckles on a light background, and are identifiable in both left and right windows of the stereo image.

17

Figure 20. Fluidisation of 1,000 marked particles. Preliminary 2D particle tracking in the image plane.

Results using this higher resolution and improved algorithms were found to be excellent. Instead of being able to track only a small subset of particles (10), the new methods are able to reliably position and track a majority of the 1,000 particles of the fluidization. Position results and the mean circulation pattern averaged from the tracking results are shown on Fig. 21. This pattern would be very difficult to measure by any other means.

Figure 21. Fluidisation of 1,000 marked particles. Left: matched particle positions; right: mean particle circulation in a meridian plane. The particle diameters in the left plot are the true diameters of the solid

18

The results obtained demonstrate the unique potential of the methods developed in the project. By combining Refractive Index Matching, laser marking, and stereo imaging techniques, it is demonstrated for the first time that dense 3D systems of particles bathed in a liquid can be tracked, unveiling the details of the rapid motions inside the dispersion. We are currently working on getting these results ready for a paper submission.

19

Goals attained

From the point of view of method development, the goals of the research project have been fully attained:

1) Matched liquid and solid materials allowing optical access to the interior of high-concentration granular dispersions have been successfully identified and prepared. 2) A novel double recirculation flume was successfully designed, assembled and tested. 3) Post-processing techniques for the detailed statistical analysis of large flow field data sets

have been derived and implemented.

4) Robust algorithms for the 3D tracking of particle positions and motions in systems of many particles have been developed and tested.

A highlight of the present work is the successful tracking of a majority of particles from a dense fluidisation of 1,000 marked transparent particles in a liquid of matched refractive index. To the best of our knowledge, tracking results of this kind had never been obtained before.

From the point of view of method application, however, the original goals of the present research project have only been partially attained. It was initially hoped that the various techniques described above could be jointly applied, in order to obtain detailed 3D measurements of liquid-granular surges in the double recirculation flume. This unfortunately turned out to be unrealistic within the three-year time frame. Each component was successfully developed, but the integration of the various components was not achieved in time. A major factor was the late date at which the double recirculation flume became operational (in the third year of the project). This was due to necessary modifications in the design: first, to insure that all flume materials could withstand the RIM liquid; secondly, to resolve fatigue problems by installing large wheels instead of small rollers.

As a result, the imaging techniques could not be applied to surges observed in the recirculation flume (it was not yet online). Instead, other applications were used for testing of the imaging algorithms: 1) downslope flows observed in existing devices at other universities (through cooperation with the National Central University, Taiwan, and the Università degli studi di Trento, Italy) were used to test and exploit the statistical flow field analysis methods; 2) fluidisation cell tests performed at our own laboratory were used to test the RIM techniques and 3D stereo imaging algorithms.

Addressing this limitation, an integrated exploitation of the methods developed in the present project is currently in progress, thanks to continued support by the NSC.

20

Appendixes

Four papers (two published and two tentatively accepted for publication) document the most mature aspects of the techniques underlying the present research effort. They are included as appendixes to the present report.

Capart, H., Young, D.L., and Zech, Y. (2002) Voronoï imaging methods for the measurement of granular flows. Experiments in Fluids 32(1), 121-135

Spinewine, B., Capart, H., Larcher, M., and Zech, Y. (2003) Three-dimensional Voronoï imaging methods for the measurement of near-wall particulate flows. Experiments in

Fluids 34(2), 227-241.

Perng, T.H., Capart, H., and Chou, H.T. Granular motions and patterns in slow uniform flows driven by an inclined conveyor belt. Tentatively accepted, Granular Matter.

Armanini, A., Capart, H., Fraccarollo, L., and Larcher, M. Rheological stratification of liquid-granular debris flows down loose slopes. Tentatively accepted, Journal of Fluid

VoronoõÈ imaging methods for the measurement of granular flows

Abstract A set of digital imaging methods derived from the VoronoõÈ diagram is proposed and tested on various liquid±granular ¯ow applications. The methods include a novel pattern-based particle-tracking algorithm,as well as estimators of the three-dimensional granular concentra-tion from two-dimensional images. The proposed algo-rithms are able to resolve individual grain motions even for rapid shear ¯ows involving dense,¯uctuating granular ensembles. Full automation is achieved,allowing the der-ivation of accurate statistics from large sets of individual measurements,as well as the construction of complete sets of grain trajectories. Results are presented for different applications: homogeneous ¯uidization,steady uniform debris ¯ow,and unsteady debris surges.

1

Introduction

Flows of granular materials are notable for the diversity of their behaviours and their involvement in a wide range of geophysical and industrial processes. While much about

them can be learned from computational simulations (Campbell 1990; Kalthoff et al. 1997),there is a corres-ponding need for detailed experimental measurements. A variety of techniques are available for multiphase ¯ow measurements (Bachalo 1994; Chaouki et al. 1997),and those applied more speci®cally to granular ¯ows include force and impact measurements (Savage and McKeown 1983; Zenit et al. 1997),acoustic probes (Bennett and Best 1995),tracked transmitters (Dave et al. 1999),and mag-netic resonance imaging (Nakagawa et al. 1993). The present work focuses on digital imaging techniques applied to the analysis of monocular image sequences. These are typically acquired by ®lming a ¯ow free-surface from above or by imaging the ¯ow through a transparent wall. As granular ¯ows are characterized by high solid fractions,opacity of the material generally prevents optical penetration beyond a few grain diameters. In many instances,this does not prevent one from making mean-ingful measurements,and it is not necessary to resort to sophisticated techniques such as magnetic resonance imaging or refractive-index matching (Cui and Adrian 1997) to penetrate the ¯ow interior. One would therefore wish to resort to the powerful imaging techniques of experimental ¯uid mechanics (Adrian 1991) to provide nonintrusive,whole-®eld coverage of the ¯ow kinematics of the visible grains.

Granular ¯ows are,however,characterized by three features that pose problems to imaging velocimetry algo-rithms: dense dispersions of grains,¯uctuating motions produced by interparticle contacts,and sharp ¯ow gradi-ents on the scale of only a few grain diameters. The sim-plest methods of particle-tracking velocimetry (PTV) rely on minimum displacement matching (Guler et al. 1999) and fail when the interframe displacement becomes sig-ni®cant with respect to the mean particle interdistance. This limit is quickly reached for dense,rapid granular ¯ows. More sophisticated particle-tracking methods in-volving trajectory-based matching (Sethi and Jain 1987) are easily offset by the uncorrelated granular velocity ¯uctuations. The techniques of particle imaging veloci-metry (PIV) are robust with respect to large displacements of dense dispersions (Willert and Gharib 1991) but have dif®culties in dealing with intense shear (Huang et al. 1993). Finally,the direct-correlation (Fujita et al. 1998) or spatial ®ltering methods (Uddin et al. 1998) do not resolve individual grain motions,an information which is most important if one seeks to test microstructural theories.

Similar dif®culties are encountered for concentration estimation. When grains become closely packed,occlusion

Experiments in Fluids 32 (2002) 121±135 Ó Springer-Verlag 2002 DOI 10.1007/s003480100351

Received: 12 December 1999 / Accepted: 21 June 2001 Published online: 29 November 2001

H. Capart (&)

Fonds National de la Recherche Scienti®que and Department of Civil Engineering,

Universite catholique de Louvain,

1 place du Levant,1348 Louvain-la-Neuve,Belgium D. L. Young

Department of Civil Engineering and Hydrotech Research Institute,National Taiwan University,

158 Chow Shan Rd,Taipei 10617,Taiwan Y. Zech

Department of Civil Engineering, Universite catholique de Louvain, 1 place du Levant,

1348 Louvain-la-Neuve,Belgium

The ¯uidization cell tests were performed at the Laboratoire du GeÂnie Civil,Universite catholique de Louvain,Belgium. The uniform debris ¯ow experiments were conducted by L. Guarino at the Hydraulics Laboratory of the Civil and Environmental Engineering Department of the UniversitaÁ degli Studi di Trento, Italy. H. C. bene®ted there from the hospitality and support of professors A. Armanini and L. Fraccarollo. The dam-break measurements were obtained with the help of Y.H. Liu at the Hydrotech Research Institute,National Taiwan University, Taiwan,and with support of a travel grant from the Fonds National de la Recherche Scienti®que,Belgium.

effects hamper the correspondence between two-dimen-sional observations of the grain number density and their three-dimensional volumetric concentration. Suggestions have been made for such a correspondence (Capart et al. 1997),but have not been tested in detail. Proposals based on luminance measurements (Louge and Jenkins 1997; Kanda et al. 1999),on the other hand,require case-by-case calibration and a careful control of illumination,which is dif®cult to obtain even in laboratory conditions.

Various strategies have been adopted by researchers to cope with these dif®culties. A ®rst possibility is to di-minish tracking problems by restricting the investigations to dilute dispersions or slow deformation rates,and to avoid concentration estimation problems by dealing with two-dimensional dispersions of spheres or disks moving between two closely spaced plates (Elliott et al. 1998; Wildman et al. 1999). Another widely used approach is to track only a small proportion of coloured ``tracer'' parti-cles included in the dispersion (Natarajan et al. 1995; Liu et al. 1997; Hsiau and Jang 1998),leading to a corres-sponding loss of resolution. Finally,a third strategy con-sists in performing the tracking manually or with manual supervision (Drake 1991; Capart and Young 1998),limit-ing the number of measurements to the bounds of human patience.

In the present work,we develop and test imaging methods that are aimed at addressing the above short-comings. Both for velocimetry and concentration,the methods derive from a pattern-based principle. The idea is that the local pattern formed by neighbouring grains will remain stable over a certain time,even for a rapidly deforming,¯uctuating granular phase,and that it can therefore serve as a match template. Furthermore,it can also be used to characterize the local degree of packing of the dispersion. The speci®c tool chosen to describe local patterns is the VoronoõÈ diagram. Various properties of the VoronoõÈ diagram endow the approach with advantages over other pattern-based methods (Haynes and Turner 1992; Song et al. 1996,1999; Ruan et al. 1999),and this will be discussed in detail in the presentation of the algorithms that forms the ®rst part of the paper. In the second part,the methods are tested on selected liquid± granular ¯ows,including homogeneous ¯uidization, steady uniform channel ¯ow,and dam-break induced debris surges.

2

Principle and algorithms

Consider a sequence of images depicting the ¯ow of an ensemble of grains. Rather than directly correlating win-dows of pixel values,the present work deals with grain images by ®rst abstracting them into point-like particles. The granular ensemble is thus reduced to a set of feature-points corresponding to grain centroids that are

dispersed in the image plane. Particle identi®cation can be achieved using a variety of segmentation or ®ltering methods,and the particular algorithm used for the applications is described in Appendix 1. After this oper-ation,the analysis of the motions and patterns of sets of point-like particles constitutes the general object of the VoronoõÈ methods.

2.1

VoronoõÈ construction and properties

Let a dispersion of n feature-points Pioccupy positions riˆ (xi, yi) (i ˆ 1,¼, n) in the two-dimensional plane (Fig. 1). The VoronoõÈ construction designates the tiling of the plane into n polygonal regions (or ``cells'') such that each polygon Vi encompasses the region of the plane that is nearest to Pithan to any other feature-point. Feature-points characterized by VoronoõÈ cells sharing an edge are termed ``natural neighbours'' of each other. The graph that connects natural neighbours further de®nes a second tesselation,dual to the VoronoõÈ diagram: the Delaunay triangulation,composed of triangles Dj.

The two constructions present many remarkable prop-erties that have led to recent applications in a variety of ®elds including cell biology (Marcelpoil and Usson 1992), computational mechanics (Braun and Sambridge 1995), astronomy (Bernardeau and van de Weygaert 1996),and molecular hydrodynamics (Espanol 1998). A general over-view of the VoronoõÈ diagram,its properties,and applica-tions is given in Okabe et al. (1992),and a mathematical introduction in Preparata and Shamos (1985). In the present context,the VoronoõÈ and Delaunay diagrams are of interest primarily because they provide useful local structures,or ``tokens'',which can be exploited for pattern characteriza-tion and matching (Ahuja 1982). The most obvious among these structures are the VoronoõÈ cells and Delaunay trian-gles themselves,with properties such as area and perimeter.

Structures of a second type,most useful for pattern-matching,de®ne local ``neighbourhoods''. The triplets of feature-points that are vertices to a common Delaunay triangle provide one such local neighbourhood (Song et al. 1999). Another structure of this type,key to the present work,is the VoronoõÈ 1-star (Fig. 2). The ®rst vertex star (or simply 1-star) Siof feature-point Piis de®ned as the set of its natural neighbours,including itself (Senechal 1995). The 1-star can be visualized as a ``star of spokes'' origi-nating at feature-point Pi. (Likewise,the 2-star of feature-point Pican be de®ned as the set of feature-points that are natural neighbours to the 1-star of feature-point Pi,and so on). As shown in detail in the present work,the VoronoõÈ 1-star constitutes a very useful token for pattern-matching in a particle-tracking context. A comparison with other neighbourhood de®nitions is provided in Sect. 2.2.

Fig. 1. VoronoõÈ diagram (Ð) and its dual,the Delaunay trian-gulation (- - -),constructed on a random dispersion of feature-points Pi(o)

The VoronoõÈ and Delaunay tesselations are character-ized by the following desirable properties:

1) Geometric properties. First,except for unlikely degenerate cases,the construction is unique for a given set of feature-points. Secondly,the construction is local, i.e.,the position of remote points does not affect the local structure of the diagram (Sibson 1981). The tokens de-scribed above can thus be construed as reliable descrip-tors of the local point pattern. Thirdly,the construction is adaptive with respect to variations in the local density of points. In particular,natural neighbours of a feature-point tend to ``surround'' it evenly,regardless of a pos-sible density gradient in one direction or another (Ahuja 1982).

2) Kinematic properties. The VoronoõÈ diagram is stable to continuous deformation. This means that if the feature-points move along continuous trajectories,then the shape of the VoronoõÈ cell will also evolve gradually,and neigh-bourhood relations will change one by one (except for unlikely degenerate cases). If the trajectories of the fea-ture-points are known,this can be used to continuously update the VoronoõÈ diagram rather than reconstruct it from scratch at discrete times (Albers et al. 1998). This property makes the VoronoõÈ 1-star a particularly attractive token for ¯owing dispersions of particles,as branches of the star can be expected to deform gradually,except for isolated topological events that will affect one branch of the star at a time. Delaunay triplets are not so attractive in this respect,since they are affected much more severely by these ``swap'' events (Fig. 3). The VoronoõÈ diagram also remains stable in the case of addition or suppression of a feature-point,as in the case of particle occlusion.

3) Computational properties. The construction presents a number of desirable computational properties. First, algorithms are available that can construct the VoronoõÈ diagram in O(nlogn) operations (where n is,as before,the number of feature-points). The particular algorithm used in the present work is based on Fortune's plane sweep method (Fortune 1987). Secondly,once the construction is obtained,a nearest-neighbour search (useful for matching operations) can be performed in O(logn) operations (Preparata and Shamos 1985). Finally,local reconstruction of the VoronoõÈ diagram (after removal of a spurious par-ticle image,for example) can be performed in O(1) oper-ations. These features make it possible to devise ef®cient implementations of the algorithms detailed in Sect. 2.2.

2.2

The matching problem

Consider now a set of moving particles,with positions sampled on two image frames acquired at instants t1and t2. Supposing that particle images are indistinguishable from each other (particles are of the same size,colour, etc.),the information available is limited to the two sets of feature-point positions ri,1(i ˆ 1,¼, n1) and rj,2(j ˆ 1,¼, n2). Based on this information,the objective is to match particles between the two frames. Formally,one seeks a pairing {i(k), j(k)} associating particle image Pion frame 1 and particle image Pjon frame 2 to one and the same physical particle k. Particle velocity vectors can then be obtained from:

vk ˆrj…k†;2_t2 r_t1i…k†;1 …1† If the time interval Dt ˆ t2)t1 is small enough,then a reasonable algorithm is to associate to each feature-point Pi,1of frame 1 the nearest feature-point of frame 2. This is the ``minimum displacement'' algorithm (JaÈhne 1995; Guler et al. 1999),which can be written formally match…Pi;1† ˆ min_P

j;2 …distP…Pi;1; Pj;2†† : …2†

Statement (2) is to be interpreted as follows: for a given point Pi,1on frame 1,the best match among points Pj,2on frame 2 is chosen as the one that minimizes the ``point-distance'' distP(Pi,Pj) ˆ [(xi)xj)2+(yi)yj)2 ]1/2,i.e.,the Fig. 3a, b. Local VoronoõÈ diagram (Ð) and Delaunay triangles (- - -) constructed on a deforming set of points (o) at times a t1

and b t2. Because of the deformation,``swap'' events may change

the local con®guration. Delaunay triangles are much more strongly affected by these topological changes than the VoronoõÈ cells (and their associated 1-stars)

Fig. 2. The VoronoõÈ 1-star Si(Ð) associated with one of the

vertices Pi(o) of the VoronoõÈ diagram (- - -)

standard Euclidean distance. Although such a general notation is not necessary at this stage,it will be useful in the remainder of the presentation as other distances are to be introduced.

To make the ``minimum displacement'' algorithm slightly more robust and avoid multiple matches to one and the same feature-point,a re®nement is to set up the following optimization problem: ®nd the global pairing {i(k), j(k)} that minimizes the objective function

X min…n1;n2†

kˆ1

distP…Pi…k†;1; Pj…k†;2† …3†

that is,a sum over all selected pairs of the distances be-tween matched particles Pi(k),1 and Pj(k),2. This is a stan-dard bijective graph optimization problem,which is dif®cult (and computationally expensive) to solve thor-oughly for large numbers of points. An approximate solution is,however,easily found using the Vogel algo-rithm. The algorithm consists in considering for each particle image the best match and the second best match, then constructing a reasonable global optimum by picking particle pairs in the order of maximum difference between ®rst and second best choices.

The simple match algorithm sketched above presents two shortcomings. First,from a computational point of view,it is unnecessarily expensive by requiring the tabu-lation of the distance between each feature-point of frame 1 and every feature-point of frame 2,however far apart from each other. It would thus be desirable to narrow down the match candidates using some suitable criterion. Secondly,a more serious problem is that the algorithm fails as soon as the interframe displacement of particles Dr is not small with respect to the mean particle interdistance dr. In that case,nearest points from two successive images are likely to correspond to two different physical particles, and a ``goodness-of-match'' criterion more robust than the minimum displacement is required. In Capart (2000),it is shown that for the case of a rigidly-moving random dis-persion of points,one must have Dr/dr < 0.35 for mini-mum displacement matches to be reliable at the 90% con®dence level. When the minimization improvement (3) is used,the condition relaxes to Dr/dr < 0.7,still quite a stringent constraint in the case of rapidly moving,dense granular dispersions.

To address these problems,a ®rst possibility is to use information from the particle's previous trajectory. Pro-vided the past trajectory is known,a predictor step can be used to estimate the likely position of a particle on the next frame,and restrict the match candidates to those en-countered in a limited search window around this posi-tion. Path regularity can then be used as an indicator of goodness-of-match (Sethi and Jain 1987; Jain et al. 1995). A re®nement consists in solving a complex multi-frame optimization problem. Such methods,used for instance for the radar monitoring of air traf®c (Brookner 1998),are successful for a number of PTV applications (Malik et al. 1993; Ushijima and Tanaka 1996) where particles are sparse and their trajectories well-behaved. However,when particle density is high and velocity ¯uctuations are important,these methods are highly unstable. This was

experienced by the authors in their ®rst attempts at imaging granular ¯ows (Capart et al. 1997; Capart and Young 1998),with the consequence that close manual supervision was necessary to obtain reasonable results.

A second possibility is to use spatial rather than tem-poral information,and focus on pattern regularity rather than path regularity. Provided that particle neighbour-hoods can be de®ned,these neighbourneighbour-hoods can be used, on the one hand,to restrict the match candidates,and on the other hand,to extract point patterns that can be compared from one image to the next. These two opera-tions can be respectively referred to as the ``screening'' and ``selection'' processes. Various ways of de®ning neigh-bourhoods have been developed (a review is given in Ahuja 1982) and applied to particle tracking applications. With reference to Fig. 4a±b,a common de®nition is to take a circular window of radius R around feature-point Piof image 1,and consider as match candidates the feature-points of image 2 belonging to this window (or R-neigh-bourhood). A second circular window (possibly with a different radius r) can then be used to de®ne a template of neighbours (or r-star) of particle image Pi,which can be compared with similar templates associated with the var-ious match candidates on frame 2 (Fig. 4a±b). Such an approach has been used by Haynes and Turner (1992),

Fig. 4a±d. Screening (left) and selection (right) operations per-formed with circular windows (above) and VoronoõÈ 1-stars (below): a a window of radius R (Ð) is used to select the match candidates and a window of radius r (- - -) de®nes the template to be matched; b goodness-of-match is evaluated by examining the area overlap between disks of radius q centred on each feature-point of the template; c the 1-star Sj0,2(Ð) of the

minimum-displacement match Pj0,2is used to select the match candidates,

and the point's 1-star Si,1(- - -) de®nes the match template;

d goodness-of-match between stars Si,1and Sj,2is estimated by

looking at interpoint distances according to the Hausdorff-type measure of (5). Window-based algorithms (a) and (b) require the de®nition of three parameters R, r, q,whereas VoronoõÈ methods (c) and (d) are non-parametric and are free to adapt to variations in the local state of the point dispersion

Lloyd et al. (1992),and Ruan et al. (1999),choosing a particle area-overlap criterion for goodness-of-match between the templates (Fig. 4b). A similar approach could be devised with neighbourhoods (and stars) de®ned in terms of the K (k) nearest neighbours (de®ning K-neigh-bourhoods and k-stars). However,both these approaches have the drawback that the radii or numbers of nearest-neighbours considered must be de®ned a priori,and cannot adapt to variations or local gradients in point density (Ahuja 1982). The data structures necessary for constructing and handling the neighbour lists,which take the same form as the Verlet lists of molecular dynamics (Frenkel and Smit 1996),are also rather heavy to manipu-late in comparison with a structure such as the VoronoõÈ diagram.

Recently,Song et al. (1996,1999) proposed to use the triangles of the Delaunay tesselation as match tokens. The present work proposes to resort to the VoronoõÈ 1-star both to de®ne a search neighbourhood and to provide a match template.

2.3

VoronoõÈ match algorithm

With reference to Fig. 4c±d,the proposed match algorithm involves the VoronoõÈ construction for both screening and selection operations. It consists ®rst in ®nding on frame 2 the minimum-displacement match Pj0,2corresponding to feature-point Pi,1 of image 1. The match candidates are selected as the feature-points on image 2 that belong to the VoronoõÈ 1-star Sj0,2of point Pj0,2. The VoronoõÈ 1-star of each of the candidates is then compared with the 1-star of feature-point Pi,1 to evaluate the goodness-of-match. For-mally,the best match of point Pi,1 is thus

match…Pi;1† ˆ min_P

j;22Sj0;2…distS…Si;1; Sj;2†† …4†

where Sidesignates the VoronoõÈ 1-star of point Pi,and distS(Si,Sj) represents a suitably de®ned ``star-distance'' re¯ecting the degree of discrepancy between the patterns formed by two stars Siand Sj. For such a function,it is proposed to choose the median of the distances between the star extremities once the star centres have been made to coincide (Fig. 4d). Formally,this can be written distS…S1; S2†ˆmedian

k1ˆ1::m1 k2minˆ1::m2 …rk1;1 r0;1† …rk2;2 r0;2†

 

…5† where vertex k1ˆ 0 (resp. k2ˆ 0) is the centre of star S1 (resp. S2),vertices k1ˆ 1..m1(resp. k2ˆ 1..m2) are the extremities of star S1(resp. S2),and |r| ˆ (x2+y2)1/2is the usual Euclidean norm. Expression (5) features three suc-cessive stages: 1) a translation making the star centres coincide,allowing their shapes to be compared in a com-mon frame of reference; 2) an inner loop based on the minimum norm,whereby for each extremity rk1,1of star S1, the nearest extremity of star S2is found,and their distance from each other is added to a list; and 3) an outer loop based on the median norm,whereby the median value of the list is adopted as an overall measure of the discrepancy between the two stars S1and S2. The above de®nition allows comparison between the shapes of any two stars,without

requiring that they have the same number of branches. It corresponds to the Hausdorff distance (e.g.,Goodrich et al. 1999),except for the use of the median rather than the maximum in the outer norm. The choice of the median is made in order to increase the robustness of the compari-son: even the correctly matched stars (from a physical point of view) are expected to differ signi®cantly at some of their extremities (when topological swap events result from ¯ow deformation,when particle occlusion occurs,when particles reach the boundaries of the domain,etc.). Note that the star-distance de®ned above is ``directed'',i.e.,in general,distS(S1,S2) ¹ distS(S2,S1). An undirected equiva-lent is easily de®ned as 1/2[distS(S1,S2) + distS(S2,S1)],but this was not found necessary in the present context.

Rather than simply setting j(k) ˆ match(i(k)),an opti-mization problem can again be set up to improve the matching and avoid multiple pairings with one and the same particle. This can be performed once again by seeking the global matching {i(k), j(k)} that minimizes the objective function

X min…n1;n2†

kˆ1

distS…Si…k†;1; Sj…k†;2† …6†

where an arbitrarily large distance value is attributed to pairings involving match candidates that the screening process has not retained (i.e.,such that Pj,2 =2 Sj0,2). Approximate solutions to the optimization problem can again be obtained using the Vogel algorithm. It is shown in Capart (2000) that,for a rigidly moving random dispersion of points,screening the match candidates on the basis of their membership in the 1-star of the minimum displace-ment neighbour leads to a limitation in range of

Dr/dr < 1.5 at the 90% con®dence level. This is more than twice the range of the optimized minimum displacement algorithm,but may still be too restrictive in certain situ-ations. In such cases,the VoronoõÈ algorithm can be gen-eralized to consider the VoronoõÈ 2-star or higher-order stars of the feature-points. This was,however,found unnecessary in the applications examined in Sect. 3. Qualitatively,the VoronoõÈ match algorithm has the effect of pairing the neighbour feature-points that are characterized by similar VoronoõÈ cells on successive im-ages. For the granular ¯ows examined,it has been veri®ed that VoronoõÈ cells tend to remain suf®ciently stable along trajectories to allow the procedure to be successful. A close-up view of VoronoõÈ diagrams corresponding to ac-tual granular ¯ow images illustrates this point in Fig. 5. More applications and results are detailed in Sect. 3. 2.4

Natural-neighbour spatial filtering

The algorithm above evidently does not succeed in matching all particle images correctly. It is thus desirable to have an automatic ®ltering procedure to remove mis-matches from the set. In the next paragraph,a procedure for performing such a ®ltering based on multiple-image trajectories is presented. Relying on information derived from a single pair of image frames,it is also possible to devise spatial ®ltering procedures. Various proposals exist for this purpose. Most of them involve comparing each

velocity vector with a local average obtained from neigh-bouring feature-point velocities. An interesting alternative (JaÈhne 1995; Song et al. 1999) consists in exploiting ¯ow invariants associated with polygonal shapes: an example is the area of Delaunay triangles constructed on passive particle tracers,which should be conserved in a diver-gence-free two-dimensional ¯ow. In the present work,we are interested in the ¯ow of a granular phase that is not locally divergence free (the granular dispersion is com-pressible),and we adopt the ®rst approach. The con-struction of local averages,however,can again be based on the VoronoõÈ construction.

The adopted approach resorts to the so-called ``natural-neighbour weights'',introduced by Sibson (1981) for in-terpolation purposes. The idea consists in subdividing each VoronoõÈ tile Viinto m subtiles Wk,corresponding to the intersections of Viwith the VoronoõÈ polygons obtained if feature-point Piis removed from the set (Fig. 6). Weights can then be de®ned as

kk ˆarea…Wk†_area…Vi† …7†

These weights present various remarkable properties outlined in Sibson (1981). In particular,they can be used to derive C2 _{interpolants (Sibson 1981),which have been} used in a particle-tracking context by Lloyd et al. (1995). By construction,the weights kksum to 1 and re¯ect the proximity of point Piwith each of the extremities of its associated 1-star (Fig. 6). It is thus possible to de®ne the

``star-averaged'' velocity associated with feature-point Pias the weighted sum

v i ˆ

Xm kˆ1

kkvk …8†

where indices k ˆ 1..m refer to the m extremities of the 1-star of Pi,and the vkare obtained from (1) once the matching has been performed. The above-de®ned average does not involve the velocity viat the star-centre,yet it presents the remarkable property of yielding vi*ˆ viif the velocities are governed by ®rst-order variations (i.e., constant gradients) around feature-point Pi(Sibson 1981). One can then de®ne a second-order spatial difference operator as

D2_viˆ vi vi

area…Vi† …9†

where Viis the VoronoõÈ cell associated with point i. Expression (9) can be veri®ed to reduce to the classical ®nite difference Laplacian if the local tesselation happens to be an orthogonal grid. A robust spatial ®ltering proce-dure can then be devised by imposing |D2_v

i| < tol (where tol is a user-de®ned tolerance),and removing one by one the outliers by starting with those that create the maxi-mum local ``curvature''. Implementation of the subtiling algorithm is complicated,but reasonably ef®cient (Sibson 1981). Once the subtiling is obtained,however,it can also be used for other purposes,for example to interpolate velocities onto a regular grid (as in Lloyd et al. 1995). 2.5

Trajectory reconstruction and temporal filtering

If a sequence of many images is recorded,displacement data obtained from each pair of frames according to the above algorithms can,of course,be concatenated and edited by using the more complete temporal information. Once successive match correspondences are known,pure concatenation of trajectory segments is trivial. Parallel to this operation,however,trajectories can be checked and edited using split±merge operations. The principles involved are similar to those of path-coherence based particle tracking (e.g. Ushijima and Tanaka 1996). The important difference is that these operations can now be performed on the basis of a pattern-based ``skeleton'' of trajectories rather than from scratch.

Temporal ®ltering can be performed ef®ciently using the 4-point stencils shown in Fig. 7 a. A subtrajectory Ti 1..Ti‡1 composed of points Pi)1..Pi+2 located at positions ri)1..ri+2for successive times ti)1..ti+2is consi-dered. Focusing on the central link Tiof the subtrajectory, a robust local measure of path coherence is given by the ``trajectory-distance''

distT…Pi; Pi‡1† ˆ min…D ; D‡† ˆ min jr_i rij; jr‡

i‡1 ri‡1j

…10† where ri)ˆ 2ri+1)ri+2and r‡_i‡1 ˆ 2ri)ri)1are respectively backward and forward extrapolations based on the sur-rounding segments. One can again ®lter out likely mis-matches by imposing distT(Pi,Pi+1) < tol. The motive for using the minimum norm is illustrated in Fig. 7b: to avoid Fig. 6. Subtiling of the VoronoõÈ cells used to de®ne the Sibson

weights

Fig. 5a±c. Overview of the VoronoõÈ match algorithm: a image of a granular ¯ow abstracted into point-like particles (+); b VoronoõÈ diagrams constructed on these points (thin lines) and on the points of the next image of the sequence (thick lines); c dis-placement vectors (true scale) obtained by matching the VoronoõÈ 1-stars. The result amounts to matching visually concordant VoronoõÈ cells

discarding correct links,only those presenting excess discrepancies with respect to both the forward and back-ward extrapolations are ®ltered out.

Once trajectories are ``split'' at their weak links,one can attempt ``merge'' operations for trajectory ends that ap-pear to ®t together well according to the same criterion. An elegant way to perform the ``split'' and ``merge'' operations in one single sweep consists in recasting the problem as one of ®nding for each pair of frames the global match {i(k), j(k)} that minimizes the objective function

X min…n1;n2†

kˆ1

distT…Pi…k†;1; Pj…k†;2† …11†

and keeping only the matches that satisfy a given toler-ance. By starting from the ``trajectory skeleton'' obtained by the VoronoõÈ matching process,trajectory-based im-provements can be made in substeps that do not require the handling of more than two frames at a time. In this way,one avoids the algorithmic complexity involved in keeping track of various multiframe trajectories.

In a granular ¯ow context,trajectory information can be quite useful,for instance,to visualize the granular motions or to extract autocorrelation statistics. These are of theoretical importance in granular ¯ows,for instance,to estimate collisional frequencies or compare self-diffusion statistics with those of computational simulations. 2.6

Granular concentration estimation

Because they tile the plane without gaps,the VoronoõÈ cells and Delaunay triangles can be used to provide a local measure of the planar density g of feature-points observed in an image. One only needs to take the reciprocal of the areas,which corresponds to 1 particle in the case of a VoronoõÈ cell,and 1/2 particle in the case of a Delaunay triangle. When dealing with two-dimensional granular dispersions (for instance,disks or spheres constrained to move between two closely spaced walls),this constitutes an immediate estimate of concentration. In contrast,when

three-dimensional dispersions of grains are imaged through a sidewall,the estimation of volumetric concen-tration / from two-dimensional images is a nontrivial task. Various authors (Louge and Jenkins 1997; Kanda et al. 1999) propose to use luminance information for this purpose. Even in laboratory conditions,however,illumi-nation conditions are dif®cult to control precisely. In the present work,it is therefore sought to use exclusively the information associated with particle positions,which is much less dependent on lighting.

Situations are different for the dilute and dense cases. In the dilute limit,for instance,in the case of sparse disper-sions of particles restricted to a thin sheet (e.g.,a laser light sheet),one can relate the 3D concentration to the 2D density of visible particles by assuming that their positions obey a Poisson process (Adrian 1991). A short analysis sketched in Appendix 2 yields:

g ˆ_pd4₂ 1 exp 3₂/Dz_d

 

 

…12† where g is the number of visible particle centroids by unit image surface, / is the volumetric solid concentration,and where it was assumed that spherical particles of diameter d are contained in a viewing volume of thickness Dz. From (12) it is immediately apparent that as soon as the ratio /Dz/d grows,occlusion effects cause the surface density g to lose its sensitivity to changes in /.

Fortunately,when the dispersion becomes dense, excluded volume effects intervene,and the particle positions are no longer governed by a random Poisson process. When the mean particle interdistance becomes of the order of the diameter,neighbouring particles are forced to organize with respect to each other in a type of glassy state (Allen and Thomas 1999). This creates short-range correlations between grain positions and opens a possibility of trying to link solid concentration with local descriptors of particle con®gurations. Various approaches can be used to construct such descriptors (for a review and comparison,see Wallet and Dussert 1998). The VoronoõÈ diagram provides once again a possible tool,assessed by Wallet and Dussert (1998) to be one of the best from the point of view of discrimination power and stability (the best results for their tests were,however,obtained with minimal spanning tree approaches). In the present work, three VoronoõÈ-based indicators are tested.

The ®rst is an estimator for the point density g:

gi ˆ_area…Vi†1 …13†

where Viis the VoronoõÈ cell enclosing point i. The second is the ``roundness factor'',which provides a descriptor of the shape of VoronoõÈ polygons:

ni ˆ_{‰perimeter…Vi†Š}4parea…Vi† 2 …14† Finally,a third indicator re¯ects the local VoronoõÈ ``area disorder'' and is de®ned as:

vi ˆ_{1 ‡ r=l}1 …15†

Fig. 7a±b. Trajectory ®ltering: a four-point stencil used to eval-uate the local path coherence; b by resorting to the minimum norm in (10),link T3can be discarded without discarding link T2

where r2 _ˆXm kˆ1 kk‰area…Vk‡ Wk† lŠ2 …16† and l ˆXm kˆ1 kkarea…Vk‡ Wk† …17†

where indices k ˆ 1..m designate the m natural neighbours of point i,and where the weights kkand polygons Vkand Wkare de®ned as in (7) (see also Fig. 6). The ``area dis-order'' estimator can be interpreted as the normalized variance of the areas of the ``petals'' of the VoronoõÈ ``¯ower'' shown in Fig. 8. The estimator is a local version of the global area disorder de®ned by Marcelpoil and Usson (1992). The weighted averages are chosen in anal-ogy with the treatment of mixture densities in Duda and Hart (1973). Note that in contrast with (13),estimators (14) and (15) are dimensionless and scale-invariant.

For the three indicators above,relationships with con-centration / are sought in the following power-law form:

/ /_rcpˆ g g_rcp !a ˆ g0a …18a† / /_rcpˆ n n0 nrcp n0  b ˆ n0b …18b† / /rcpˆ v v0 vrcp v0 !_c ˆ v0c …18c† where indices ``rcp'' and ``0'' designate the state of random close packing and the dilute state,respectively; /rcp 0.64 for spheres (Allen and Thomas 1999); n0 0.72 and v0 0.80 are obtained from Monte-Carlo simulations (as they concern the dilute state,the values do not vary with particle shape); random close-packing values grcp, nrcp,and vrcp,as well as exponents a, b,and c remain to be determined from calibration tests.

3

Applications and results

The algorithms above were tested on a number of liquid± granular ¯ow experiments,which also motivated and

guided the developments. In the present report,three applications are selected as testing grounds for the algo-rithms: (i) homogeneous ¯uidization tests; (ii) steady uniform debris ¯ow; (iii) transient debris surges. These contrasted applications make it possible to highlight different features of the proposed methods.

3.1

Homogeneous fluidization

In order to test the concentration estimators,it is desirable to obtain states of homogeneous dispersion for various values of the solid fraction. A most convenient way of achieving this is to set up ¯uidization cell experiments. These consist of subjecting a static array of grains to an upward ¯uid ¯ux,causing the granular assembly to expand until the voidage is such that the interphase drag balances the submerged weight of the granular phase. At some stage,when concentrations become lower than the random loose-packing concentration /rlp 0.55,the grains become mobile with respect to each other and un-dergo disordered ¯uctuating motions. Provided one can avoid the regions of instabilities (Batchelor 1988),a sta-tistically homogenous state is obtained in which the grains randomly explore a variety of local con®gurations. Such conditions are the ideal ones to test and calibrate pattern-based concentration estimators.

In the present work,tests were performed with light spherical grains ¯uidized by a water current. The grains (arti®cial pearls) have a diameter d ˆ 6.1 mm and a rela-tive density q/qwˆ 1.048,where q is the density of the granular material and qwis the density of water. The cylindrical ¯uidization cell has a height of 25 cm and a diameter of 10 cm. The lower part of the cell is ®lled with a 5-cm-deep layer of small lead spheres in order to diffuse the incoming water current. Finally,the cylinder is ®tted with a plane rectangular observation window of dimen-sions 5 by 10 cm allowing one to ®lm the dispersion without optical distortion. Varying the upward water velocity in a range of 1 to 3 cm/s enables observations of ¯uidized granular dispersions having concentrations be-tween / ˆ 0.2 and / ˆ 0.55. The relationship bebe-tween concentration and ¯uidizing ¯ux was veri®ed to corre-spond reasonably well with the Richardson±Zaki empirical relation (Richardson and Zaki 1954). As the granular motions are rather slow in this case,tracking of the par-ticles constitutes a simple matter,and was carried out only to verify that the state of the dispersion was close to sta-tistically homogenous. The measured kinetic energy of the velocity ¯uctuations was found to agree well with the relationship proposed by Batchelor (1988) for homoge-neous ¯uidization conditions. These veri®cation data are presented in Capart (2000).

The main purpose of the ¯uidization experiments is to test the pattern-based concentration indicators. Once particle centroids are located with the algorithms of Appendix 1,the VoronoõÈ diagram can be constructed on these feature-points,and estimators (13)±(15) obtained for each of the VoronoõÈ cells. To prevent edge effects from biasing the statistics,the criterion of Kenkel et al. (1989; see Okabe et al. 1992) is used to discard VoronoõÈ cells close to the image boundaries. It consists in excluding Fig. 8. The union of a neighbour cell Vkwith subtile Wkforms a

``petal'' of the VoronoõÈ ``¯ower''. The variance of the petal areas can be used to construct a local ``area disorder'' estimator 128