Three-dimensional Voronoı¨ imaging methods for the measurement
of near-wall particulate flows
B. Spinewine, H. Capart, M. Larcher, Y. Zech
Abstract A set of stereoscopic imaging techniques is proposed for the measurement of rapidly flowing disper-sions of opaque particles observed near a transparent wall. The methods exploit projective geometry and the Voronoı¨ diagram. They rely on purely geometrical principles to reconstruct 3D particle positions, concentrations, and velocities. The methods are able to handle position and motion ambiguities, as well as particle-occlusion effects, difficulties that are common in the case of dense disper-sions of many identical particles. Fluidization cell experi-ments allow validation of the concentration estimates. A mature debris-flow experimental run is then chosen to test the particle-tracking algorithm. The Voronoı¨ stereo methods are found to perform well in both cases, and to present significant advantages over monocular imaging measurements.
1
Introduction
Flows of disperse phases are involved in a wide variety of situations of scientific and engineering interest. These include liquid-entrained gas bubbles, aerosols, dry gran-ular flows, fluidized beds of particles, and liquid–saturated
particulate currents. In many of these situations, the individual elements (bubbles or grains) can be approxi-mated as rigid bodies undergoing distinct motions, and the dynamic system can be abstracted into an evolving configuration of particle positions (Campbell 1990; Zhang and Prosperetti 1994; Kang et al. 1997). The system be-havior in the dense limit is often of particular interest, featuring active particles interacting with each other, and adopting flow-induced preferential arrangements (Fortes et al. 1987; Savage and Dai 1993; Rouyer et al. 2000). The present work addresses some of the experimental chal-lenges associated with such flows through the development of novel stereo imaging techniques.
We consider the following typical measurement set-up (Fig. 1). The flow of interest involves a dense ensemble of identical opaque particles immersed in a transparent fluid and imaged by twin cameras through a transparent plane (e.g., flume side-wall). Because of the high concen-tration of particles, optical peneconcen-tration within the bulk is limited to a depth of a few grain diameters. Measurements sought include the 3D positions of the visible grains, the in-plane (i.e., parallel to the image plane), and out-of-plane (normal to the image out-of-plane) particle velocities, and the volumetric particle concentration. Two main problems must be faced in deriving such measurements: (1) particle-pairing ambiguities, which hamper both ste-reoscopic matching and velocimetric tracking; this is especially true for dense, rapidly sheared dispersions of identical particles undergoing irregular motions; (2) par-ticle occlusion, which biases estimates of the volumetric particle concentration; in the dense case, the apparent concentration of visible particles is not equivalent to the actual concentration of physical particles.
In some instances, it is possible to circumvent these difficulties. Problems caused by matching and tracking ambiguities can be diminished by focusing only on the positions and motions of a subset of marked particles or tracers. Problems caused by occlusion, on the other hand, can be bypassed by using refractive index matching techniques (Cui and Adrian 1997) or full volumetric scanning as in nuclear magnetic resonance imaging (Phillips et al. 1992; Nakagawa et al. 1993; Seymour et al. 2000). A number of other non-intrusive methods have been used for the dynamic characterization of opaque solid–liquid flows, including positron emission particle tracking (Fairhurst et al. 2001; Wildman et al. 2001), and diffusing wave spectroscopy (Menon and Durian 1997). These techniques, however, are not always feasible and present their own drawbacks. The use of subsets of marked Received: 8 August 2002 / Accepted: 12 September 2002
Published online: 19 December 2002 Springer-Verlag 2002
B. Spinewine (&)
Fonds pour la Recherche dans l’Industrie et l’Agriculture and Department of Civil and Environmental Engineering, Universite´ catholique de Louvain, place du Levant 1, 1348 Louvain-la-Neuve, Belgium
E-mail: [email protected] Tel.: +32-10-472123
Fax: +32-10-472179 H. Capart
Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan M. Larcher
CUDAM – Dipartimento di Ingegneria Civile ed Ambientale, Universita` degli Studi di Trento,
Via Mesiano 77, 38050 Trento, Italy
Y. Zech
Department of Civil and Environmental Engineering, Universite´ catholique de Louvain, place du Levant 1, 1348 Louvain-la-Neuve,
Belgium
Experiments in Fluids 34 (2003) 227–241 DOI 10.1007/s00348-002-0550-4
particles, as in PEPT for example, can yield good
Lagrangian statistics for the motion, but misses the spatial correlation and particle arrangement effects that act on the scale of a few particle diameters. Refractive index match-ing, in addition, is only possible for liquid–solid mixtures of very special optical properties. The approach of the present paper is to stick with more widely available test materials and digital cameras, and to deal with the ambi-guity and occlusion problems in the context of stereo imaging methods.
Stereo imaging techniques applied to dilute disper-sions (passive tracers seeding a fluid flow) have been described in the review of Adrian (1991) and in the studies of Malik et al. (1993), Ushijima and Tanaka (1996), and Virant and Dracos (1997). A method to obtain concen-tration estimates from statistical distributions associated with stereo imaging techniques is described in Murai et al. (2001). Monocular imaging methods, on the other hand, have been applied to dense flows by various experiment-ers. Most of these works have been restricted to 2D ana-logues, e.g., a monolayer of disks or spheres sandwiched between two parallel plates (Drake 1991; Azanza et al. 1999; Wildman et al. 1999), or to the imaging of a subset of marked particles (Natarajan et al. 1995; Hryciw et al. 1997; Hsiau and Jang 1998). In Capart et al. (2002), monocular imaging methods were proposed for the measurement of the discrete kinematics of full sets of visible grains in a dense, rapidly sheared dispersion imaged near the side-wall. The issues associated with stereo imaging of dense 3D dispersions do not seem to have been specifically ad-dressed before.
The approach adopted follows the work of Capart (2000) and Capart et al. (2002). The core of the proposed methods consists in exploiting the special properties of Voronoı¨ diagrams (Ahuja 1982; Okabe et al. 1992). These define spatial tessellations of the 2D plane or 3D space into cells centered around individual feature points (see Sect. 3). After abstraction of digital images into discrete sets of point-like particle positions, the Voronoı¨ diagram is used for the three main steps of the analysis: (1)
stereoscopic matching of particles; (2) estimation of volumetric concentration; (3) pattern-based particle tracking. Presentation of these various methods consti-tutes the central section of the present paper. This is preceded by an outline of the 2D particle-identification technique and generic 3D ray tracing concepts, which are needed as basic tools. The final part of the paper is devoted to the application of the methods to two partic-ulate flows in water. The first case, a homogeneous flui-dized bed of light particles, is of particular interest for the validation of concentration estimation methods. The second case is a mature debris flow of PVC particles fea-turing steep concentration and velocity gradients, and constitutes a challenging test for the particle-tracking techniques. Monocular imaging results for similar flows were presented in Capart et al. (2002), while some pre-liminary stereo results were described in Douxchamps et al. (2000) and Spinewine et al. (2000).
2
Particle identification and ray geometry
This first section outlines the basic methods required to pinpoint particles on 2D images and to reconstruct their positions in 3D space. While these basic methods are not new, their presentation sets the stage for the more original Voronoı¨ developments of the next section. The following notations are adopted throughout the paper: upper-case letters denote 2D imaging quantities measured in pixels and ticks (intervals of time separating successive images); lower-case letters denote 3D world quantities measured in meters and seconds. In both cases, vectors are indicated in bold face.
2.1
Particle identification
The first step of the analysis consists in the localization of particle centroids on individual images. For each instant at which synchronized images are acquired under two viewpoints A and B, one seeks to identify sets of particle positions Rð ÞiA
n o
and Rð ÞjB
n o
on the two digital frames. Particle images show up as white blobs of a certain size against a dark background (see Fig. 2a). An image neighborhood associated with a particle can be approxi-mated by a Gaussian gray-level function centered on the particle centroid, with a diameter D that scales with the pixel diameter of the particle. Such bell-like regions are highlighted by convoluting the image with a Laplacian-of-Gaussian (Mexican hat) filter of width D (Ja¨hne 1995), as shown in Fig. 2b and c. Local brightness maxima of the highlighted images are then identified iteratively using a ‘‘dish-clearing’’ algorithm: a global maximum is found, then a Gaussian bell of diameter D is subtracted from the neighborhood gray-level values; a new global maximum is found, and so on. The position of each maximum is finally refined to subpixel accuracy by fitting a second-degree interpolation surface around the discrete pixel position. The resulting set of particle centers is shown in Fig. 2d. The expected rms accuracy on the X and Y image coor-dinates obtained in this fashion is of the order of 0.25 pixel Fig. 1a, b. Stereoscopic imaging of near-wall dispersion of particles:
ausing two synchronized sensors; b using a single sensor with system of mirrors. (1) transparent side-wall; (2) CCD camera; (3) light source; (4) mirrors
(Veber et al. 1997). More details about the procedure are given in Capart et al. (2002).
2.2
Transformation between world and camera coordinates A transformation is now required to relate the set of 2D image coordinates of any point P to its world coordinates. Define rP=[xP yPzP]Tas the world coordinates of point P
(Fig. 3), and Rð ÞPA ¼ X A ð Þ P Y A ð Þ P h iT
as the 2D image coordi-nates of point P associated with the camera viewpoint A (point P¢ in left image plane F). The transformation is obtained by modeling the image formation as a central projection from a virtual camera focal point A onto the image plane F (Tsai 1987; Jain et al. 1995). Conserving the alignment of points, this projective transformation constitutes a reasonable approximation of the image formation process and facilitates tremendously the ray tracing and matching operations (Trucco and Verri 1998; Faugeras 1999). For each viewpoint A, one can then specify a matrix [A(A)] and a vector b(A)such that a Xð ÞA Yð ÞA 1 2 4 3 5 ¼ Ah ð ÞAi xy z 2 4 3 5 þ bð ÞA ð1Þ
where a is a scalar parameter best interpreted in the context of Eq. (5) below. Matrix [A(A)] and vector b(A), on
the other hand, must be calibrated from a set of points Pk
for which both the world coordinates [xkykzk]Tand the
image coordinates Xkð ÞAY A ð Þ k
h iT
are known. Eliminating parameter a from (1), one obtains for each calibration point Pktwo linear equations in the unknown coefficients
að ÞijA and bð ÞjA of [A(A)] and b(A): xkaðAÞ11 þ ykaðAÞ12 þ zkaðAÞ13 xkXkðAÞa
ðAÞ 31 ykXðAÞk a ðAÞ 32 zkXðAÞk aðAÞ33 þ b ðAÞ 1 X ðAÞ k b ðAÞ 3 ¼ 0; ð2Þ
xkaðAÞ21 þ ykaðAÞ22 þ zkaðAÞ23 xkYkðAÞaðAÞ31 ykYkðAÞaðAÞ32
zkYkðAÞa ðAÞ 33 þ b ðAÞ 2 Y ðAÞ k b ðAÞ 3 ¼ 0: ð3Þ
A minimum of six calibration points Pk(12 equations)
are thus required to determine the 12 unknown coeffi-cients. Since the system is homogeneous, one needs to append one more equation, e.g.,
bð Þ1A ¼ 1; ð4Þ
to avoid the trivial solution with all zero coefficients. The system is further over-determined if more than six calibration points are used, hence a least-square procedure is needed to obtain an optimal solution. In practice, it is recommended to use a large number of calibration points with positions distributed evenly inside the viewing volume.
Once calibrated, Eq. (1) can be used to obtain image coordinates from known world coordinates. Conversely, a point P having its projection P¢ with image coordinates Rð ÞPA ¼ X A ð Þ P Y A ð Þ P h iT
under viewpoint A is known to belong to a ray AP (or AP¢ as seen in Fig. 3) defined by parametric equation
rAPð Þ ¼ ra Aþ asAP; ð5Þ
where a is a free parameter, and vectors rAand sAPare
given by rA¼ Að ÞA h i1 bð ÞA; ð6Þ sAP¼ Að ÞA h i1 X A ð Þ P YPð ÞA 1 2 6 4 3 7 5: ð7Þ
Fig. 3. Imaging geometry: physical point P and its image projec-tions P¢ and P¢¢ on the two stereo views. Rays emanate from focal points A and B of the two central projections. The trace of the epipolar plane comprising P is shown in dashed lines. Inset: due to imperfections in the imaging geometry, the two rays may not perfectly intersect
Fig. 2a–d. Particle-identification procedure: aimage fragment; b image a after low pass filtering; c image b after high pass filtering; dparticle positions at brightness maxima of c
The vector rAis the position of the focal point A of the
projection, and lies at the intersection of all rays associ-ated with viewpoint A. It is obtained by setting a to zero in (1). The vector sAP, on the other hand, specifies the
direction of the particular ray which pierces the image plane at P¢. It is obtained as the difference rP–rAor by
setting a to unity in (1).
When the imaged scene is immersed in a liquid and seen from the outside through a transparent wall, the projection is altered by refraction effects. In the simple case of a liquid-bathed scene observed through a plane wall, and under some restrictions relative to imaging configuration, however, the affine character of the pro-jection can be preserved to a very good approximation, provided that the calibration step is performed in refrac-tion condirefrac-tions identical to those of the actual experiments (calibration target immersed in the fluid). This can be verified through a detailed analysis of ray refraction, and has been checked empirically to lead to negligible errors in De Backer (2001).
2.3
Ray intersection and epipolar constraint
Consider now a physical particle P having unknown po-sition rPin 3D space. Let projections of the particle center
on the left (focal point A) and right (focal point B) views have known image positions Rð ÞPA (P¢ on plane F) and Rð ÞPB (P¢¢ on plane Y). Referring to (5), the corresponding rays are given by parametric equations
rAPð Þ ¼ ra Aþ asAP; ð8Þ
rBPð Þ ¼ rb Bþ bsBP; ð9Þ
where a and b are the two free parameters. The 3D posi-tion of the particle can now be retrieved by finding the intersection of the two rays AP¢ and BP¢¢ (Fig. 3). Due to imperfections of the imaging process, this intersection cannot be exact and instead we look for the point of closest encounter between the two rays. Let M and N be the points on both rays separated by the smallest inter-distance. Their positions
rM¼ rAþ aMsAP; ð10Þ
rN¼ rBþ bNsBP ð11Þ
are found by solving the following linear system in the two unknown parameters aMand bN:
sTAPsAP sTAPsBP sTBPsAP sTBPsBP aM bN ¼ ½rB rATsAP ½rB rATsBP : ð12Þ
The midpoint of the segment joining M and N constitutes an approximation of the true particle position P. The length ‘ of the segment, on the other hand, measures the distance of closest encounter between the two rays and provides an estimate of the quality of the approximation. They are given by
rP¼ 1 2½rMþ rN; ð13Þ ‘¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rN rM ½ T½rN rM: q ð14Þ Suppose now that the only available information about a particle P is the position Rð ÞPA of its image P¢ on the left
view (plane F). Can anything be said about projection Rð ÞPB
of the same physical point on the right view (plane Y)? Physical point P is known to belong to the plane ABP¢ containing both focal points A and B, and image P¢ on the left view. Its parametric equation is
r a;bð Þ ¼ rAþ asAPþ b r½ B rA ð15Þ
where a and b are again two free parameters. This plane AP¢B is called the epipolar plane (see Fig. 3). Its inter-section with the right image plane is a straight line , called the epipolar line, which constitutes the locus of all possible projections, in the right image plane, of a physical point having projection Rð ÞPA in the left image
plane (P¢ on F). Parametric equation R(B)(G) of the epi-polar line is easily found by projecting two points of the left ray AP¢ on plane Y, via (1).
The true projection Rð ÞPB (P¢¢ on Y) is known to lie
along this line, or, due to inevitable imperfections of the imaging process, to fall somewhere close to the epipolar line, i.e.,
Rð ÞPB Rð ÞB ð Þ;C ð16Þ
for a certain unknown value of parameter G. This restriction is known as the epipolar constraint, and can be used advantageously to guide the search for a correspondence between two stereoscopic views. It will be exploited below to accelerate the stereoscopic matching step.
Fig. 4a–c. Voronoı¨ diagrams (2D and 3D): a planar Voronoı¨ diagram {Vi}; b Voronoı¨ vertex star Si(thick lines) associated with cell Vi(thin lines) surrounding point Pi; c a 3D Voronoı¨ cell vkaround point pk;
to each natural neighbor (d) corresponds a particular facet of the Voronoı¨ polyhedron vk, perpendicular to a segment of the related vertex star Sk
3
Voronoı¨ imaging methods
With reference to Fig. 4a, consider a set of feature points {Pi} placed at positions {Ri} in the 2D plane or {ri} in 3D
space. The Voronoı¨ diagram is a geometrical construction that divides the space into a set of polytopes, or cells, surrounding each feature point. Each Voronoı¨ cell Vi(in
2D) or Vi(in 3D) encompasses the region that lies closer to
Ri(respectively ri) than to any other feature point of the
set. The construction presents many useful properties, discussed in a general context in Okabe et al. (1992). In the present imaging context, three main properties are of in-terest:
1. The diagram organizes the dispersion of points into a network of neighborhood relations that can be used to speed up nearest-neighbor queries (Preparata and Shamos 1985). This will prove very useful to accelerate the stereoscopic matching step, which involves a repeated search for points of the image plane located in the vicinity of a given epipolar line.
2. The diagram defines a natural partition of space into non-overlapping regions, each containing a single particle. The inverse area or volume of the cell thus defines a local concentration estimate at the smallest possible scale. This can be used to sample statistical distributions in a finer and more robust way than by binning observations into an arbitrary partition (Bernardeau and van de Weygaert 1996). The property is exploited below to estimate particle concentrations. 3. The Voronoı¨ diagram defines local patterns of
neigh-boring feature points, which can be used as match templates (Ahuja 1982; Song et al. 1999). For moving dispersions of points, these patterns remain stable over a certain duration of time, a property exploited for 2D particle tracking by Capart et al. (2002). The present work extends this idea to 3D.
Stereoscopic matching, volumetric sampling, and particle-tracking operations are now outlined in three separate subsections.
3.1
Stereoscopic matching 3.1.1
Basic procedure
For a single particle seen on two stereo views, the 3D particle position can easily be retrieved using the approximate ray intersection procedure of Sect. 2.3. To perform this operation when the imaged scene features a large number of particles, it is necessary first to solve the correspondence problem: find which particle image on one view corresponds to which one on the other. As required when all particles are identical, it can be solved on the basis of purely geometrical information by minimizing ray intersection discrepancies.
Consider two stereo images of a large number of identical particles (Fig. 1). From the particle-identification step, sets of 2D particle positions Rð ÞiA
n o
and Rð ÞjB
n o
have
been located on the left and right views, respectively. Using the inverse projective transformation of Sect. 2.2, 3D rays can then be traced back through each of these image points. Ray bundlesnrð ÞiAð Þaoand rð ÞJBð Þb
n o
are thus obtained for the two views, and solving the correspon-dence problem amounts to finding a pairing j(i) subject to the constraint that no more than one ray of one view can be paired with any one ray of the other (some rays can be left out, however, as a result of partial occlusion effects and differing viewing prisms). Once pairing j(i) has been formed, the set of measured 3D positions ^rrf g is easilyi
obtained from (13).
The pairing itself can be found by an optimization procedure, minimizing the objective function
X
i
‘ij ið Þ ð17Þ
where the distance of closest encounter ‘ijgiven by (14)
is taken as a ‘‘goodness-of-match’’ between any two rays ri(a) and rj(b). This constitutes a standard bipartite graph
optimization problem (e.g., Kim and Kak 1991), which is difficult (and computationally expensive) to solve thoroughly for large numbers of points. An approximate solution can however be found using the Vogel algorithm. It consists in considering for each ray the best match and the second best match, then constructing a reason-able global optimum by picking ray pairs in the order of maximum difference between first and second best choices.
3.1.2
Voronoı¨ epipolar screening
For large numbers of particles, the procedure above becomes prohibitively expensive in terms of computational time and memory allocation. This is because it requires the computation of discrepancies ‘ijfor all possible pairings of
rays, and the handling of a full matrix in the Vogel opti-mization procedure. To limit those expenses, it is advan-tageous to conduct a first screening of possible match candidates by resorting to the epipolar constraint (16). This can be done very efficiently using the 2D Voronoı¨ diagram. With reference to Fig. 5a, consider a particle i viewed on the left image at position RðAÞi . In the right image plane (Fig. 5b), let us construct the 2D Voronoı¨ diagram Vn ð ÞjB o on the set of particle positions Rð ÞjB
n o
. The projection Rð Þj ið ÞB of particle i in the right image is known to lie in the vicinity of epipolar line Rð ÞiBð Þ. A simple way to screenC match candidates is then to retain only those particle images j that have their Voronoı¨ cells Vð ÞjB pierced by the epipolar line (see Fig. 5b, where 12 such candidates are marked). If particles are expected to lie within a limited depth range [ymin... ymax] (where y is the normal to
side-wall coordinate), parametric line Rð ÞiBð Þ becomes a lineC segment, and the set of match candidates can further be restricted (to the five candidates shown in black in the example of Fig. 5a, b).
If N denotes the total number of visible particles, the procedure reduces considerably the number of candidates
to be considered from N to pffiffiffiffiN (screening along a line instead of the full area). The 2D Voronoı¨ diagram must be computed only once at a cost of NlogN
operations (Okabe et al. 1992). Epipolar line piercing then requires only a nearest-neighbor search and neighbor-by-neighbor propagation within the Voronoı¨ diagram, for a combined cost of pffiffiffiffiN. Overall, the Voronoı¨ epipolar screening reduces the cost of stereo matching from N2to N3/2. For the experiments described in Sect. 4, this
amounts to a 10–30-times speed-up over the brute-force approach.
Figure 5c and d show the reconstructed 3D dispersion of particles resulting from application of the complete stereo matching procedure. Starting from sets of image positions Rð ÞiA
n o
and Rð ÞjB
n o
, the procedure is seen to be successful at pinpointing the 3D positions ^rrf g of mosti
of the visible particles. An indication of the robustness of the approach is that partial occlusion relationships (visible by eye on the image but not exploited in the analysis) are well conserved by the reconstruction. It is also seen in Fig. 5d that, for such a densely packed dispersion, occlusion effects prevent the retrieval of particle positions beyond a depth of a few grain diameters. This issue is addressed again in the next subsection.
3.2
Volumetric sampling 3.2.1
Binning and Voronoı¨ sampling estimates of volumetric concentration
Supposing one can pinpoint a full set of 3D particle po-sitions {rk} inside the viewing region, this set can be
exploited to yield further spatial information about the granular dispersion. The quantity of most interest is the local volumetric concentration /(r). The simplest way of obtaining estimates for / is to subdivide the viewing volume into an arbitrary spatial partition {Wj}, then to
count the number of particle centroids njfalling into each
box Wj of the partition. The local concentration in each
box is then estimated as: ^
/
/j¼ njVp volume Wj
ð18Þ
where VP¼16pd3 is the volume of a single solid particle.
This is the so-called ‘‘binning’’ procedure, widely used to estimate statistical distributions.
Instead of defining an arbitrary partition, an alternative approach makes use of the 3D Voronoı¨ diagram. By con-struction, the Voronoı¨ diagram partitions the region into a set of cells {Vk}, each one containing a single particle.
Local concentrations can then be estimated as: ^
/
/k¼ VP
volume Vð kÞ
: ð19Þ
Advantages of the approach are that Voronoı¨ partition {Vk} is more natural than arbitrary partition {Wj}, and that
concentrations are sampled at the locations {rk} where
particle velocities are defined. To apply such estimates in the case of a dispersion bounded by a plane side-wall, it is convenient to construct an extended Voronoı¨ diagram on the basis of the set {rk} complemented by set {r¢k} of its
mirror images on the other side of the wall (refer to Fig. 8a, discussed more in details in Sect. 3.2.3). This guarantees that near-wall cells will be bounded by facets belonging to the plane boundary, hence allowing concen-tration estimate (19) to take the wall into account. 3.2.2
Occlusion effects
Unfortunately, for a dense dispersion of opaque grains, the set of measured particle positions ^rrf g picked up by thei
stereo procedure is not the full physical set {rk}, but a
subset of this including only the particles which are visible under both camera viewpoints. Occlusion effects intervene, whereby particles in the front hide from view particles in the back. To get a feel for how occlusion affects concen-tration estimates, it is useful to analyze first a simple model, introduced by Capart et al. (2002) in a monocular context and extended here to stereo acquisition.
Let us assume a 3D dispersion of spherical grains of diameter d, characterized by constant volumetric con-centration / in the near-wall region. Disregarding excluded volume effects due to impenetrability of the solid particles and side-wall, one may assume that the Fig. 5a–d. Principle of epipolar stereoscopic matching. a Considered
particle center on the left view; b identification of the match candidates on image 2, whose Voronoı¨ cells are pierced by the corresponding epipolar line (reduced to the segment in solid line); creconstructed 3D positions of the particles. Color is assigned according to depth, bright for near-wall particles, darker for deeper ones; d side view of c, solid line indicating wall position
particle centroids are distributed inside the 3D volume according to a homogeneous Poisson process, with an average number of particles per unit volume equal to lP¼ /V1
P . Considering the idealized stereo imaging
configuration shown in Fig. 6a, a particle centroid will be visible on both stereo views only if no other particle center invades the two viewing cylinders shown as hatched surfaces. Given the depth of the particle ypand the viewing
angle w under which it is seen by each camera, the hatched volume can be approximated by
V yð Þ P pd 2 4 yPþ y2P tan w d if yP2 tan wd ; ð20Þ V yð Þ P 3pd 3 16 tan wþ pd2 2 yP d 2 tan w if yP2 tan wd ; ð21Þ in which angle w is assumed to be small, and where the complex sub-volume of intersection of the two cylinders is only roughly estimated. The probability that a particle is viewed under both viewpoints is then given by Poisson distribution (e.g., Adrian 1991; Okabe et al. 1992) P 0 particle in V½ ¼ðlPVÞ
0
0! exp lð PVÞ ¼ exp lð PVÞ: ð22Þ The average volumetric concentration of visible parti-cles ^//is thus only a fraction (22) of the actual one / and
can be expressed as the following function of depth y: ^ / /ð Þ ¼ /P 0 particle in V yy ½ ð Þ ¼ / exp /V yð Þ VP : ð23Þ As expected, the two concentrations are equivalent right next to the side-wall (y0). Away from the side-wall, the concentration of visible particle drops as a result of the rising probability of occlusion.
Figure 7 compares the predictions of (23) with mea-surements ^//j yj obtained from actual fluidization cell
experiments. The measured functions were obtained by binning actual stereo measurements into partition {Wj}
given by Wj rjyj 1 2dy< y <yjþ 1 2dy ; ð24Þ
i.e., layers of small thickness dy parallel to the side-wall, then applying (18). The measured profiles of Fig. 7b ex-hibit clear humps and troughs. This is not an artefact, but a consequence of excluded volume effects (neglected in the simple occlusion analysis). Physical particle centers are
Fig. 6a, b. Occlusion analysis. Shaded volumes are the regions that must be free of occluding grains for a given particle center to be seen under: a two stereo views; b a single monocular view
Fig. 7a, b. Distribution of visible particles centers in function of the depth (normal to wall) for two characteristic concentrations (thick line 57%, thin line 26%). a Theoretical distribution derived from a 3D Poisson process; b experimental distribution gathered from fluidiza-tion cell tests. Wall posifluidiza-tion corresponds to the left of the graphs; intervals between successive graduations on the horizontal axis correspond to one particle diameter
located preferentially at depths corresponding to certain multiples of particle diameter d, in a crystalline-like order induced by the side-wall boundary. These effects are es-pecially strong for the denser packing (/0.57). The decay curves, which underlie these quasi-periodic variations, on the other hand, show the effect of occlusion. This expo-nential-like decay is seen in Fig. 7 to be reasonably well captured by the simple occlusion model introduced above. In particular, it is observed on both sets of results that away from the side-wall, the concentration of visible parti-cles may actually be higher for the case of lower concen-tration (/0.26) than for the densely packed case (/0.57). This is entirely a result of occlusion, and explains why one must be careful in trying to derive actual concentrations / from measured concentrations of visible grains ^//.
Another signal that caution must be exercised comes from examining the consequences of the simple Pois-son occlusion model outlined above in the case of mon-ocular imaging (Fig. 6b). In that case, the viewing region, which must be free of occluding particles, simplifies to a cylinder, and the surface particle density g obtained by projecting the images of all the visible particles onto the side-wall (i.e., as seen by the monocular camera) is given by (Capart et al. 2002)
^ g g¼ 4
pd2: ð25Þ
The key feature of (25) is that the resulting estimate ^ggis independent of the actual particle concentration /. The result is surprising because it is very tempting to assume that one can estimate grain concentrations by simply counting the surface density of particles seen on monoc-ular images, then resorting to some ad hoc calibration to convert ^gginto /. This obvious estimate, however, is geo-metrically completely insensitive to /. Only because of attenuation in illumination with depth does some sensi-tivity appear, making the measurements hinge upon illu-mination conditions, which are usually not well controlled in the laboratory. Such questionable estimates are nev-ertheless found in the literature including, we must confess, some of the authors’ earlier work (Capart et al. 1997).
3.2.3
Surface Voronoı¨ sampling of the near-wall concentration Because of occlusion, one should not expect to be able to measure concentrations inside the bulk of the opaque disperse phase. What should be possible, however, on the basis of the stereo position results ^rrf g, is to estimate thei
near-wall volumetric concentration ^
/
/0ffi / y 0ð Þ: ð26Þ
As seen in Fig. 7, occlusion effects are so severe for large concentrations that some care must be exercised even to attain this more modest objective. One possibility examined in a preliminary work (Spinewine et al. 2000) is to resort to estimate ^ / /i¼ VP volume 1Vð iÞ ð27Þ
where the {Vi} are the 3D Voronoı¨ cells constructed on the
basis of the ^rrf g and their mirror images ^rri f g (Fig. 8a). Toi
limit the effects of occlusion, estimates (27) are then sampled only within the near-wall Voronoı¨ cells, i.e., the cells {Vk}, which share a face with the side-wall. This yields
reasonable results up to moderate concentrations, for which particles are visible till a depth of several diameters. However, for larger concentrations (beyond /0.4), some second- and third-row particles become unseen (Fig. 8b), leading to unrealistically large volume of front-row Vor-onoı¨ cells. As a consequence, the concentration for dense packing is significantly under-predicted by (27).
To remedy this problem, it is proposed to resort to the following modified surface particle-density estimate: ^ g g0;i¼ 1 area Að 0ð ÞVi Þ ð28Þ where A0(Vi) is the plane face A0which the 3D cell Vi
shares with the wall boundary (or equivalently shares with its mirror image). The advantage of measuring a surface density rather than volumetric density is that the frontal face of the near-wall Voronoı¨ cells is not affected by occlusion. Estimate (28) is further assumed to hold uniformly over face A0(Vi). When averaging over a
certain extent of the wall surface, with area A large with respect to the area of an individual cell face, we have that
^ g g0 h i n0
A ð29Þ
where the brackets Æ æ denote a spatial average, and n0
is the number of near-wall particles picked up within area A, i.e., the number of visible particles which have a Voronoı¨ cell sharing a face with the wall within that area.
Stereoscopic surface estimate (29) differs from mon-ocular surface estimate (25) in one fundamental regard. In the monocular case (25), all visible particles are counted regardless of the depth. In the stereoscopic case (29), by contrast, only the near-wall particles (defined in the Voronoı¨ sense) are counted. The stereo near-wall estimate is therefore biased neither by occlusion effects (which Fig. 8a, b. Volumetric sampling: a mirror image construction of a Voronoı¨ diagram in which the plane side-wall is embedded; bincomplete near-wall Voronoı¨ cell closures resulting from occlusion of back row particles. (Note: whereas 3D diagrams are actually used in the experiments, 2D analogues are shown here for illustration purposes)
affect the volumes of the 3D Voronoı¨ cells) nor by visible particles located away from the wall (which affect the monocular estimate).
3.2.4
Monte Carlo estimation of the stereological coefficient A final step is required, which is to convert surface estimate
^ g g0
h i into volumetric estimate ^//0. This relationship can
be predicted a priori if the 3D distribution of particle centers within the measurement volume is assumed to derive from a homogenous Poisson process, defining once again lpas the average number of particles per unit
vol-ume. In that case, it follows from dimensional consider-ations that the relconsider-ationship between Æg0æ and lpmust be of
the form g0
h i ¼ vl2=3p ð30Þ
where Æg0æ is the expected surface density of near-wall
particles and v is a non-dimensional constant: the so-called stereological coefficient.
While it is not impossible that an exact value for con-stant v could be derived theoretically, this is beyond our level of expertise. We thus resorted to Monte Carlo sim-ulations (e.g., Ross 1990) to obtain an approximate value. Repeated simulation runs were conducted, each run involving the following steps:
1. With the aid of a random number generator, a 3D Poisson process of given intensity lpis simulated
within a box of finite size;
2. The resulting dispersion of points is reflected on the other side of one of the box walls (= the side-wall); 3. The 3D Voronoı¨ diagram is computed and the near-wall
cells identified;
4 The surface density g0æ is estimated by counting
near-wall particles or averaging the areas of the cell faces A0(Vi) coinciding with the mirror plane.
For step 4, estimates are sampled in a central region only in order to avoid boundary effects. Except for step 1, the procedure closely emulates the one applied to the ac-tual stereo measurements. A batch of 1,000 Monte Carlo runs was carried out, each run involving 2,000 parti-cles. The result of these computations is v=0.92±0.01, where the bounds correspond to a 95% confidence interval (see Ross 1990).
Recalling that lp¼ /Vp1, relation (30) implies
^ / /0¼ VP ^ g g0 h i v 3=2 ; ð31Þ
which is our estimate for the near-wall volumetric con-centration ^//0 in terms of the averaged wall-sampled sur-face number density ^h i. Strictly speaking, relation (31)gg0 with constant v0.92 should be expected to hold only in the dilute limit. For dense dispersions, excluded volume effects will make the distribution of particle centers differ from the results of a Poisson process. As verified by the experiments of Sect. 4.1, it will turn out nevertheless to constitute a reasonable approximation over the entire range of concentrations. Such a favorable situation is en-countered in other stereology problems (Lorz 1990; Okabe
et al. 1992). Estimate (31) is remarkable in that it derives from purely geometrical considerations. Up to moderate concentrations, it requires no ad hoc calibration. This contrasts with the much less favorable situation faced when trying to extract concentration estimates from monocular images (Capart et al. 2002).
3.3
Pattern-based particle tracking
The third operation that exploits the properties of the Voronoı¨ diagram is the particle-tracking step. Sets of 3D particle positions at successive times are first acquired by repeated application of stereo matching to each frame of a movie sequence. Let {ri,m} and {ri,m+1} be two such sets of
particle positions sampled at successive times tm and
tm+1=tm+Dt. Particle velocities vi(tm+1/2) can be estimated
by expression vi;mþ1=2¼
rj ið Þ;mþ1 ri;m
Dt ; ð32Þ
provided one can first ‘‘connect the dots’’, and establish a pairing j(i) between positions ri,mand rj(i),m+1belonging to
one and the same physical particle. When dealing with a moving dispersion of many identical particles, the main problem consists in establishing this correspondence: finding which particle on one snapshot corresponds to which one on the next. The particle-tracking problem can thus be seen as a time-domain variant of the stereo matching problem addressed previously.
For dilute particle dispersions or slow motion, the correspondence problem can easily be solved simply by pairing together the particles on one frame and the next, which are nearest to each other (see, e.g., Guler et al. 1999). For dense dispersions or rapid motion, however, legiti-mate pairing candidates may travel further away and the minimum displacement criterion breaks down. An alter-native approach derives from the following observation: while individual particles are identical to each other, the local arrangements that they form with their neighbors are unique and may be preserved by the flow long enough to serve as a basis for tracking. Particle pairing can then be performed based on pattern similarity.
Following Capart et al. (2002), we again resort to the Voronoı¨ diagram to implement such pattern-based track-ing. Nearby particles are paired according to the geomet-rical similarity of their Voronoı¨ cells. This similarity is estimated as follows: ‘‘stars’’ are first constructed by con-necting a given particle center to its Voronoı¨ neighbors (i.e., the particles with which it shares a cell face), as il-lustrated in Fig. 4b and c for both the 2D and 3D cases. The stars belonging to two pairing candidates can then be compared for goodness-of-fit by making their centers co-incide and measuring the distances between the star extremities. Once these ‘‘goodness-of-fit’’ indicators are available for all possible pair candidates, a global optimum problem can be defined and solved in the same fashion as in the stereo matching case (see Sect. 3.1 and Eq. 17).
For two dimensions, the overall method is illustrated in Fig. 9 for a plane granular flow. While graphical repre-sentation is harder in 3D, the algorithms themselves gen-eralize straightforwardly to the third dimension. The
reader is referred to Capart et al. (2002) for a detailed presentation and comparison with alternative approaches. Three-dimensional applications and results are detailed further below.
4
Liquid–granular flow applications
The Voronoı¨ imaging methods presented in the last two sections are now applied to two different cases of immersed particulate flows, each aimed at highlighting particular aspects of the techniques. Fluidization cell experiments are chosen as a first test case, ideally suited for validation of concentration estimates. The steady uniform flow of a water–granular mixture down an in-clined open channel was selected as a second test case. Featuring strong gradients in both solid fraction and granular velocities, this second case constitutes a chal-lenging application for both velocity and concentration measurements.
4.1
Fluidization cell experiments 4.1.1
Principle and set-up
The proposed imaging techniques are now applied to fluidization cell tests. The principle of the tests is as follows. Subject to an ascending water current, a layer of loosely packed grains expands into a fluidized suspen-sion. The concentration adapts to the water flux until the mean drag balances the submerged weight of the grains. In this fluidized state, the suspended particles undergo weakly correlated fluctuating motions, exploring a variety of spatial arrangements. Most important in the present context, the average concentration of particles is spatially uniform throughout the fluidized layer. Different concen-trations can be obtained simply by tuning the fluid flow. Beyond their intrinsic interest, such homogeneous states of known concentration are ideal for testing concentration measurement methods.
The device used for the experiments is presented in Fig. 10. The cylindrical fluidization cell has a height of 25 cm and an inner diameter of 10 cm. A 5 cm deep layer of small lead spheres is placed at the bottom of the cell to diffuse the incoming water current and provide uniform fluid velocity throughout the cross-section. Above this heavy static layer comes the granular bed to be fluidized. It is composed of light spheres (artificial pearls) of relative density qs/qw=1.048 and diameter D=6.1 mm. To allow
visualization without optical distortion, a plane
rectangu-lar observation window of dimensions 5·10 cm is fitted to the cell wall.
The evolving granular dispersion is imaged by two synchronized CCD cameras placed in a stereoscopic ar-rangement. The sensors, each 256·256 pixels in size, are positioned approximately 35 cm away from the cell, and mounted ±25 cm apart from each other. The angle be-tween the two viewpoints is thus around 20. Both lenses have a focal distance of 16 mm. Particle diameters span around 10 pixels. The resulting positioning error on par-ticle centroids is estimated to be of the order of 0.2 mm.
The viewpoint calibration procedure sketched in Sect. 2.2 is carried out by imaging an ‘‘open-book’’-shaped dihedron placed within the cell filled with water. A total of 56 calibrations points is used for the least-square derivation of matrixes A(A), A(B)and vectors b(A), b(B) needed to determine the left (A) and right (B) viewpoints. The world coordinate system (x, y, z) is as follows: the x–z plane coincides with the cell wall (x taken horizontally), and y represents the out-of-plane horizontal coordinate (i.e., depth inside the cell).
Two independent series of tests were carried out under slightly different camera and lighting conditions. For each series, measurements were performed for a number of experimental runs corresponding to different solid con-centrations (see Table 1). The Richardson–Zaki empirical correlation (Richardson and Zaki 1954) was found to de-scribe well the relation between fluidization velocity and concentration.
In order to derive granular velocities and concentrations, the sets of simultaneous images are first preprocessed to
Fig. 9a–c. Pattern-based velocimetric track-ing: a original image fragment and particle positions; b Voronoı¨ diagrams built on sets of particle positions at two successive time instants; c velocity vectors built by matching of the local Voronoı¨ patterns. (Note: whereas 3D diagrams are actually used in the exper-iments, 2D analogues are shown here for illustration purposes)
Fig. 10. The fluidization cell: 1 bounding plates; 2 Acrylic glass cell shaft; 3 impermeable joints and fixation screws; 4 incoming water flow; 5 layer of lead spheres to diffuse water current; 6 fluidized bed 236
pinpoint particle positions in the two image planes. Ste-reoscopic ray matching and intersection is then used to retrieve the 3D positions of the grains. These procedures were described and illustrated in Sect. 2. We now focus on the concentration estimation process.
4.1.2
Concentration estimation
For every fluidization velocity, the height of the fluidized layer furnishes a direct and fairly accurate measure of the bulk concentration, which can serve as validation of the imaging estimates. The latter are extracted from uncorrelated snapshots acquired at a very low frame rate. Since the number of particles seen on each image is limited, this is needed to sample a variety of particle configurations and obtain good statistics by averaging concentration measurements over a number of images (16 and 30 images per run were used, respectively for the first and second series of tests).
Figure 11a presents a plot of the bulk concentration versus its imaging estimate. The latter is obtained by positioning the particle centers, constructing their 3D Voronoı¨ diagram, and sampling the surface density of near-wall particles using estimate (29). Coefficient v=0.92 derived from the Monte Carlo simulations is then used in (31) to convert surface number density ^h i into volumetricgg0
concentration ^//0. In Fig. 11a, the imaging results are seen
to be close to the line of perfect agreement up to moderate concentrations (/0.4). For higher concentrations, the theoretically based imaging estimate ^//0begins to
under-estimate the actual concentration /. This departure from perfect agreement most likely arises due to excluded vol-ume effects, which induce a quasi-crystalline arrangement of the grains at dense packing. The Poisson process as-sumption adopted to derive relation (31) and compute the value of stereological coefficient v then breaks down.
Remarkably, however, the discrepancies are not very severe. Without any adjustment, the predicted relation remains accurate to 10% over the entire range. Further-more, the data sets obtained from the two independent series of tests are very consistent with each other. Based on the fluidization data, an empirically adjusted estimate can be proposed in the form
^ / / 0¼1 jtanh 1 j ^// 0 ð33Þ where j is a dimensionless coefficient and value j=1.2 provides a good fit to the measurements. The functional form of (33), involving an inverse tangent, reduces to (31) in the dilute limit, and introduces a significant correction only at high concentrations (/>0.4) for which the Poisson process assumption breaks down. The data are plotted
again in Fig. 11b based on this empirically adjusted rela-tion. The apparent robustness and sensitivity of the ste-reoscopic imaging estimate contrasts with the much less favorable properties of the monocular imaging indicators examined in Capart et al. (2002).
4.1.3
Three-dimensional particle motions
The evolution of the particle arrangement in time can further be characterized by following the 3D motions of the imaged particles using the procedures of Sect. 3.3. In this application, the movement is relatively slow and the tracking procedure easy to carry out. A more challenging application for the tracking component of the algorithms is Table 1. Test conditions for the fluidization experiments
First series of tests
Run number 1 2 3 4 5 6 7 8 9 10 11
Fluidization velocity (cm/s) 1.97 1.05 2.46 1.78 2.20 2.54 0.93 0.72 1.18 1.84 2.95
Bulk concentration (%) 37.9 52.1 30.7 40.5 33.6 29.2 55.7 59.0 43.5 38.3 26.6
Second series of tests
Run number 12 13 14 15 16 17 18 19 20 21 22
Fluidization velocity (cm/s) 0.90 1.11 1.26 1.70 1.96 2.36 2.57 1.08 1.87 0.95 1.48
Bulk concentration (%) 55.7 50.7 47.2 41.2 37.1 32.4 29.3 51.4 38.4 54.3 43.6
Fig. 11a, b. Granular concentration: bulk measurements vs. imaging estimate. m first series of tests; n second series of tests; – line of perfect agreement. a imaging estimate ^//0(31); b empirically adjusted estimate ^// 0(33) with j=1.2
the rapidly sheared example presented in the next section. The key component here is the stereoscopic positioning step. Reliable 3D positions must be measured in order to record continuous trajectories. To damp out high-fre-quency noise due to the finite sensor resolution, particle displacements are averaged over nine successive frames (with a frame rate of 200 fps, this corresponds to an in-terval of 0.04 s, which was checked to be much less than the average time between collisions) to yield velocity and tra-jectory measurements. Sample particle paths obtained for a fluidization concentration /=0.38 are presented in Fig. 12. The fluctuating motions of neighboring grains are influ-enced by each other on a scale of a few particle diameters. This appears typical of hydrodynamic effects whereby conjugate motions of the embedding fluid transmit the influence of particle motions over a certain distance.
Mean squared velocity fluctuations along the three spatial directions x, y, z are plotted in Fig. 13 for different fluidization concentrations /. The observed velocity fluc-tuations are greater at lower concentrations, for which particles have more freedom to move around. Lower concentrations also correspond to a faster fluidization flux and greater transfer of kinetic energy to the fluctuating motions. Vertical fluctuations Æw¢2æ are found to be slightly larger than the horizontal fluctuations Æu¢2æ and Æv¢2æ.
Deviations from isotropy are not very large, however. A significant observation is that, despite the presence of the side-wall, the in-plane and out-of-plane horizontal fluc-tuations have similar magnitudes.
4.2
Uniform debris-flow experiments 4.2.1
Principle and set-up
The second test case is an open-channel flow of a highly concentrated liquid–granular mixture imaged through the
side-wall. The flow is obtained in a recirculatory flume developed at the University of Trento, Italy, for the study of torrential sediment transport and debris-flow processes (Armanini et al. 2000). Shown in Fig. 14, the device fea-tures a glass-walled flume of adjustable slope (from 0 up to 23 degrees) having the following dimensions: length=6 m; width=20 cm; wall height=40 cm. A fast conveyer belt is used to recirculate material collected at the flume outlet. This scheme achieves steady, longitudinally uniform flow conditions within the flume.
The solid grains used for the tests are PVC particles having a cylindrical shape and the following dimensions: diameter=3.2 mm, height=2.8 mm; equivalent spherical diameter=3.5 mm. The volumetric concentration in a well-packed static assembly of such grains is around /rcp0.69.
The material density is qs=1,540 kg/m3. Water is again
used as entraining fluid. The specific flow chosen for the present testing purposes is a mature debris-flow case for which the transport layer, many grains thick, fills the entire flow depth (i.e., there is no grain-free water layer in the upper part of the flow). This case was selected because it combines a wide range of concentrations, rapid veloci-ties, and intense shear, all desirable features to test the applicability and limits of the technique. Operating con-ditions for the run are: bottom slope angle=7.2; total discharge=14.5 l/s; //=delivered solid concentration=ratio of granular discharge to water discharge=49%. The latter two parameters are derived from bulk measurements performed at the end of the run by diverting the outlet flow to a trap.
To resolve individual grains in a flowing layer of such thickness, a reasonably high image resolution is needed. For flow speeds of the order of 1 m/s, a high image ac-quisition frequency and an adjustable shutter are further necessary to reliably sample the grain motions and observe blur-free particles. For this purpose a high speed CCD camera was operated at a resolution of 480·420 pixels, a frequency of 250 f/s, and an exposure time shuttered down to 1/500 s. Because of the availability of only one camera with these characteristics, stereo viewing requires a Fig. 12a–c. Example of granular trajectories obtained for a
fluidiza-tion cell experiment at a bulk concentrafluidiza-tion of 38% (run #1 of the first series). a Top view, in the x–y plane; b 3D view; c front view, in the x–z plane
Fig. 13. Velocity fluctuations vs. bulk concentration. h horizontal componentpffiffiffiffiffiffiffiffiffiffihu02i; n vertical componentpffiffiffiffiffiffiffiffiffiffihw02i; · normal-to-wall componentpffiffiffiffiffiffiffiffiffihv02i
special device composed of four mirrors arranged as shown in Fig. 14. The camera is positioned at a distance of about 1.5 m from the side-wall of the flume. The mirror set-up is then interposed midway between the camera and the flume.
4.2.2
Experimental results
Measurements obtained by application of the 3D Voronoı¨ imaging techniques described in the present paper are shown in Figs. 15 and 16. Figure 15 displays typical tra-jectories of individual grains reconstituted using the ste-reoscopic matching and velocimetry tracking algorithms. The results shown have been filtered over three successive frames (0.008 s.) to suppress high-frequency noise. Results in the x–z plane show the higher velocities attained by the grains in the upper part of the flow layer. Results in the y–z plane show that trajectories can be reconstituted farther away from the side-wall in the upper part of the flow. This is because concentration is lower there and occlusion effects are less severe. Despite the high con-centration and the irregular nature of the granular motions, long granular trajectories are successfully reconstructed by the proposed methods.
Figure 16 shows vertical profiles for the mean velocities, fluctuation velocities, and solid concentration averaged from a sequence of 512 images. In Fig. 16a, the mean ve-locities in the vertical and normal to wall directions are seen to be close to zero. The mean longitudinal velocity, on the other hand, varies from zero in the static bed layer to maximum speed at the free surface. An approximately constant shear rate ¶u/¶z is obtained in the upper part of the flow. Fluctuating motions keep a roughly constant
magnitude in this upper part, as shown quantitatively in Fig. 16b and qualitatively in Fig. 15. There is a slight peak in the fluctuations at the free surface. This appears to be due to a corner effect involving intermittent (stick-slip) motion of partly emerged grains, interacting with the side-wall because of surface tension.
In the lowermost part, the bed is virtually motionless and fluctuations should decrease to zero. Residual values there represent noise due to inaccuracies in particle po-sitioning. Small errors in position, uncorrelated from one
Fig. 14a, b. Recirculating flume of the Universita` degli Studi di Trento for the experimental study of steady uniform debris flows: a photograph; b plane view; 1 flume, 2 recirculating belt, 3 camera, 4 mirror device, 5 lighting system
Fig. 15a, b. Granular trajectories tracked over a period of 0.04 s. a front view, b ‘‘side’’ view obtained by projection onto the
y–z plane
frame to the next, translate into apparent velocity fluctu-ations. This effect is small for the in-plane components of the mean squared velocity fluctuations Æu¢2æ and Æw¢2æ.
Near the bed, a much higher magnitude is observed for the normal to wall component Æv¢2æ. Not unexpectedly, this indicates that the stereo measurements of 3D positions achieve a lower accuracy in the depth-wise direction than in directions more closely parallel to the camera image planes. The resulting noise is seen in Fig. 16b to be sig-nificant. Encouragingly, however, it does not drown out the physical signal.
Figure 16c shows the measured concentration profile obtained using adjusted estimate (33). The measured solid concentration exhibits a maximum value in the static lower layer and decreases to a minimum near the free surface. Again, corner effects perturb the measurements right at the free surface. A slightly jagged profile is ob-served in the lower part of the flow. This occurs because the bed is virtually motionless there: the particle configu-ration does not change over time hence averaging over many frames does not improve the statistics for the con-centration measurements. Overall, the shape of the profile is quite regular and features reasonable concentration values. Values close to static packing are observed in the motionless bed. The following integration over depth further yields an imaging estimate for delivered concentration // / /¼ Rz0þh z0 /huidz Rz0þh z0 h idzu ; ð34Þ
in which z0 is the elevation of the rigid bed and h is the
flow depth. The resulting value is //¼ 0:46, comparable to the value of //¼ 0:49 derived from bulk measurements at the outlet. Even if the side-wall measurements were absolutely accurate, the two values would not be expected to perfectly coincide because the flow departs somewhat from uniformity in the transverse direction. Nonetheless, similar magnitudes should be obtained, as is the case with the present stereoscopic methods. Significantly, such agreement is obtained without any adjustment to the stereological relation derived from the fluidization tests, itself very close to the theoretically predicted relation. This was far from the case for the monocular imaging methods examined in Capart et al. (2002), for which concentration estimates had to be recalibrated in an ad
hoc fashion when going from one type of experiments to another.
5
Conclusions
Three-dimensional imaging techniques were developed for the measurement of near-wall particulate flows. The methods include special matching and tracking algo-rithms, which exploit the properties of projective geometry and Voronoı¨ diagrams to reconstruct 3D particle positions and velocities. A novel estimate based on the surface density of near-wall particles was also proposed to mea-sure volumetric solid concentrations. Because they can handle occlusion effects and resolve position and motion ambiguities, the methods are particularly well suited for applications involving highly concentrated, rapidly sheared dispersions of many identical particles. Fluidiza-tion cell tests and an open-channel solid–liquid flow ex-periment were used to demonstrate the potential of the proposed techniques, as well as point out some of their limitations. It is hoped that they will prove valuable tools for the measurement of particulate flows of scientific and industrial interest.
References
Adrian RJ (1991) Particle-imaging techniques for experimental fluid mechanics. Annu Rev Fluid Mech 23:261–304
Ahuja N (1982) Dot pattern processing using Voronoı¨ neighborhoods. IEEE PAMI 4:336–343
Armanini A, Fraccarollo L, Guarino L, Martino R, Bin Y (2000) Ex-perimental analysis of the general features of uniform debris flow over a loose bed. In: Wieczorek GF, Naeser ND (eds) Proceedings of the 2nd International Conference on Debris-flow hazard miti-gation, Taipei, August 2000. Balkema, Rotterdam, pp 327–334 Azanza E, Chevoir F, Moucheront P (1999) Experimental study of
collisional granular flows down an inclined plane. J Fluid Mech 400:199–227
Bernardeau F, van de Weygaert R (1996) A new method for accurate estimation of velocity field statistics. Mon Not R Astron Soc 279:693–711
Campbell CS (1990) Rapid granular flows. Annu Rev Fluid Mech 22:57–92
Capart H (2000) Dam-break induced geomorphic flows and the transition from solid- to fluid-like behaviour across evolving interfaces. PhD thesis, Universite´ catholique de Louvain, Belgium Capart H, Liu HH, Van Crombrugghe X, Young DL (1997) Digital
imaging characterization of the kinematics of water–sediment interaction. Water Air Soil Pollut 99:173–177
Capart H, Young DL, Zech Y (2002) Voronoı¨ imaging methods for the measurement of granular flows. Exp Fluids 32:121–135
Fig. 16a–c. Vertical profiles: amean granular velocities in m/s; bvelocity fluctuations in m/s; longitudinal (bold line), vertical (solid line) and normal-to-wall (dashed line) components; cestimated granular concentra-tion. Marks on the vertical axis represent distances of 2 cm. The free surface of the flow is located at the top of the graphs
Cui MM, Adrian RJ (1997) Refractive index matching and marking methods for highly concentrated solid–liquid flows. Exp Fluids 22: 261–264
De Backer A (2001) De´veloppement et application de techniques d’imagerie digitale pour l’analyse tridimensionnelle d’un e´coule-ment complexe. BEng thesis, Haute e´cole Leonardo da Vinci (ECAM), Belgium (in French)
Douxchamps D, Devriendt D, Capart H, Craeye C, Macq B, Zech Y (2000) Three-dimensional reconstruction of the oscillatory free-surface of a flow over antidunes: stereoscopic and velocimetric techniques. In: Proceedings of the Oceans 2000 Conference, Rhode Island, USA, September 2000. IEEE, Piscataway, N.J.
Drake TG (1991) Granular flow: physical experiments and their im-plications for microstructural theories. J Fluid Mech 225:121–152 Fairhurst PG, Barigou M, Fryer PJ, Pain JP, Parker DJ (2001) Using
positron emission particle tracking (PEPT) to study nearly neu-trally buoyant particles in high solid fraction pipe flow. Int J Multiphase Flow 27:1881–1901
Faugeras O (1999) Three-dimensional computer vision: a geometric viewpoint. MIT Press, Cambridge, Mass.
Fortes AF, Joseph DD, Lundgren TS (1987) Nonlinear mechanics of fluidization of beds of spherical particles. J Fluid Mech 177:467–483 Guler M, Edil TB, Bosscher PJ (1999) Measurement of particle
movement in granular soils using image analysis. J Comp Civ Eng 13:116–122
Hsiau SS, Jang HW (1998) Measurements of velocity fluctuations of granular materials in a shear cell. Exp Therm Fluid Sci 17:202–209 Hryciw RD, Raschke SA, Ghalib AM, Horner DA, Peters JF (1997)
Video tracking for experimental validation of discrete elements simulations of large discontinuous deformations. Comput Geotech 21:235–253
Ja¨hne B (1995) Digital image processing. Springer-Verlag, Heidelberg Berlin New York
Jain R, Kasturi R, Schunck BG (1995) Machine vision. McGraw-Hill, New York
Kang SY, Sangani AS, Tsao HK, Koch DL (1997) Rheology of dense bubble suspensions. Phys Fluids 9:1540–1561
Kim W, Kak A (1991) 3D object recognition using bipartite matching embedded in discrete relaxation. IEEE PAMI 13:224–251 Lorz U (1990) Cell-area distributions of planar sections of spatial
Voronoı¨ mosaics. Mater Characteriz 25:297–309
Malik NA, Dracos T, Papantoniou D (1993) Particle tracking veloci-metry in three dimensional flows – II: Particle tracking. Exp Fluids 15:279–294
Menon N, Durian DJ (1997) Diffusing-wave spectroscopy of dynamics in a three-dimensional granular flow. Science 275:1920–1922 Murai Y, Matsumoto Y, Yamamoto F (2001) Three-dimensional
measurement of void fraction in a bubble plume using statistic stereoscopic image processing. Exp. Fluids 30:11–21
Nakagawa M, Altobelli SA, Caprihan A, Fukushima E, Jeong EK (1993) Noninvasive measurements of granular flows by magnetic-resonance-imaging. Exp Fluids 16:54–60
Natarajan VVR, Hunt ML, Taylor ED (1995) Local measurements of velocity fluctuations and diffusion coefficients for a granular material flow. J Fluid Mech 304:1–25
Okabe A, Boots B, Sugihara K (1992) Spatial tessellations: concepts and applications of Voronoı¨ diagrams. Wiley, Chichester Phillips RJ, Armstrong RC, Brown RA (1992) A constitutive equation
for concentrated suspensions that accounts for shear-induced particle migration. Phys Fluids A 4:30–40
Preparata FP, Shamos MI (1985) Computational geometry: an intro-duction. Springer-Verlag, Heidelberg Berlin New York
Richardson JF, Zaki WN (1954) Sedimentation and fluidization: part I. Trans Instn Chem Engrs 32:35–53
Ross SM (1990) A course in simulation. Macmillan, Basingstoke Rouyer F, Lhuillier D, Martin J, Salin D (2000) Structure, density, and
velocity fluctuations in quasi-two-dimensional non-Brownian suspensions of spheres. Phys Fluids 12:958–963
Savage SB, Dai R (1993) Studies of granular shear flows, wall slip velocities, ‘‘layering’’ and self-diffusion. Mech Mater 16:225– 238
Seymour JD, Caprihan A, Altobelli SA, Fukushima E (2000) Pulsed gradient spin echo nuclear magnetic resonance imaging of diffu-sion in granular flow. Phys Rev Lett 84:266–269
Song XQ, Yamamoto F, Iguchi M, Murai Y (1999) A new tracking algorithm of PIV and removal of spurious vectors using Delaunay tessellation. Exp Fluids 26:371–380
Spinewine B, Capart H, Zech Y (2000) Stereoscopic imaging measurement of solid concentrations in dense granular flows. In: Proceeding of the 10th International Conference on Transport and sedimentation of solid particles, Wroclaw, September 2000. Akademia Rolnicza we Wroclawiu, Wroclaw, Poland,
pp 497–509
Trucco E, Verri A (1998) Introductory techniques for 3-D computer vision. Prentice Hall, Upper Saddle River, N.J.
Tsai RY (1987) A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses. IEEE J Robotics Auto 3:323–344
Ushijima S, Tanaka N (1996) Three-dimensional particle tracking velocimetry with laser-light sheet scannings. J Fluid Eng 118:352– 357
Veber P, Dahl J, Hermansson R (1997) Study of the phenomena affecting the accuracy of a video-based particle tracking veloci-metry technique. Exp Fluids 22:482–488
Virant M, Dracos T (1997) 3D PTV and its application on Lagrangian motion. Meas Sci Technol 8:1539–1552
Wildman RD, Huntley JM, Parker DJ (2001) Granular temperature profiles in three-dimensional vibrofluidized granular beds. Phys Rev E 63:1311–1320
Wildman RD, Huntley JM, Hansen JP (1999) Self-diffusion of grains in a two-dimensional vibrofluidized bed. Phys Rev E 60:7066–7075
Zhang DZ, Prosperetti A (1994) Averaged equations for inviscid disperse two-phase flow. J Fluid Mech 267:185–219
