Stereoscopic and velocimetric reconstructions of the free surface topography of antidune flows

Intro Experiments Particle imaging Stereo Velocimetry Results Conclusions / references

5. Velocimetric method

In contrast with the relatively direct character of the stereoscopic method, the velocimetric technique is indirect in nature. The technique seeks to reconstruct the surface topography on the basis of measurements of the horizontal velocity field. For that purpose, the method exploits the Bernoulli relation applicable to free surface streamlines. This relation takes a particularly simple form if the flow can be approximated as a superposition of small perturbations upon a rapid mean current. To carry out the procedure, images from the top camera are first used to track particles on the water free surface and estimate their horizontal velocities. The elevation field is then obtained after suitable longitudinal and transverse averaging of the individual velocity vectors. These steps are detailed in the next paragraphs.

5.1 Elevation-from-velocity principle

The following two assumptions are taken to apply locally (on the scale of a few wavelengths) to antidune flows: 1) the time evolution of the loose sediment bed is sufficiently slow that it appears stationary to the rapidly flowing fluid, hence the flow can be taken as quasi-steady; 2) the flow can be considered inviscid (but it does not have to be assumed irrotational). The dynamics of the free surface can then be described by requiring conservation of mechanical energy along streamlines (see e.g. Batchelor 1967). Each surface streamline is subject to the following set of ordinary differential equations:

,
,
,
(15a,b,c)

,(16)

where ζ  is a curvilinear coordinate measured along the streamline projection in the horizontal plane (x, y), u = (u, v, w) is the water velocity, and η is the free surface elevation. Equations (15a,b,c) describe the streamline as a parametric curve , along which the Bernoulli equation (16) is written. Vertical velocity w can be eliminated from the Bernoulli equation using (15c). It follows that if the horizontal components of the surface velocity field (uv) of a steady flow are known (and if dissipation can be neglected), one can solve for the free surface elevation η along each streamline, up to a constant of integration. The constant can be determined by invoking a suitable boundary condition (point of known elevation along the streamline) or stationarity condition (known average elevation). This constitutes the basis of the velocimetric technique proposed here for the measurement of free surface topography. On a streamline-by-streamline basis, the above principle can be applied to a wide class of flows beyond the present antidune case: examples include flows over oscillatory bottoms, uniform currents past isolated accidents of topography, or shooting flows controlled upstream by a vane of known opening.

Before proceeding, it is useful to reduce the system to a dimensionless form. Given a characteristic scale [u] for the velocity, a characteristic wavenumber can be defined as (Reynolds 1965)

.(17)

Conversion to dimensionless variables can then be achieved by the transformation:

,
,
,
(18a,b,c)

,
,
,
,
(19a-d)

resulting in the following set of dimensionless equations:

,
,
,
(20a,b,c)

.(21)

5.2. Approximation for small amplitude oscillations over a rapid mean current

An especially simple situation arises if the free surface can be approximated as small amplitude oscillatory perturbations superposed upon a rapid mean flow. This is the case for the present antidune flows, for which one can further assume that: 1) the mean flow is oriented in the longitudinal x direction; 2) the mean flow is approximately uniform, with only a weak mean velocity gradient along the transverse y direction. In that case, it is possible to decompose velocities into mean and perturbation components and introduce a small parameter ε such that

,
,
(22a,b)

where is the mean current velocity, and measures the strength of the perturbation velocities u', v' (having characteristic scale [u']) relative to the mean current velocity (having characteristic scale [u]). One can then substitute expansions

(23)

(24)

in eqs (20)-(21) to obtain, at order 0,

,(25)

and at order 1,

(26)

where it was assumed that transverse variations in mean velocity are at most O(ε ). One can finally translate back to dimensional variables and obtain the key result

,(27)

which is particularly simple both to use and to interpret: superimposed upon a mean surface level , the free surface topography perturbations are proportional to the perturbations in longitudinal velocity u', but opposite in phase. For torrential flows, the water velocity is slowest at the crests and fastest at the troughs, an observation pointed out already by Leonardo Da Vinci in his notebooks several centuries ago (Da Vinci ca. 1500/1975).

In practice, the ratio of fluctuation velocities to mean velocities obtained in the present experiments is about 0.1, with mean velocities of about 1 m/s. This leads to errors on the perturbed elevation of the order of (1 m/s)2(0.1)2/(10 m/s2) ≈ 1 mm. More precise estimates obtained by calculating second order corrections are listed in Table 1. These second order corrections are comparable in amplitude to measurement errors. Because they are also highly sensitive to measurement noise, only the first order estimate is retained in what follows. Whatever the order of approximation, an obvious limitation of the velocimetric method is that only the topography perturbations are accessible, not the mean elevation. If needed, this has to be obtained by other means.

5.3. Horizontal velocities in world coordinates

Horizontal velocities are obtained as follows from the image sequences acquired with the top camera (viewpoint V = T). First, particle positions on the digital images are pinpointed using the algorithm of section 3.1. Particles are then tracked in the image plane using the Voronoï match algorithm of section 3.2. Index pairs j(i) are obtained, identifying matches , likely to correspond to the same physical particles on successive images.

To convert image coordinates into world coordinates , the image positions are first projected onto a plane corresponding to the mean water surface . This is carried out as follows. To each particle image viewed under camera viewpoint T is associated a ray

(28)

where and are given by equations (3a,b). On the other hand, plane Π0 can be specified by its Cartesian equation , where is the unit normal to the plane and its signed distance from the origin (with an arbitrary point belonging to the plane). Upon substitution, the parameter value corresponding to the point of intersection of the ray with the plane is given by

.(29)

In the experimental configuration, plane Π0 is close to horizontal with Cartesian equation . Once world coordinates have been approximated by their projections in plane Π0, the horizontal components of velocity vector ui,m are estimated from:

,
,
(30)

where Δt is the time interval separating successive frames.

The above procedure would be sufficient to obtain accurate measurements of the horizontal velocities if either the water free surface was perfectly planar, or if we had been able to place the camera far above the water surface with a precisely vertical optical axis. In practice however, the free surface is perturbed by the antidune oscillations and the top camera is oriented slightly obliquely (as illustrated on Fig. 2). As a result, a non-negligible bias results when the projection is performed onto a hypothetical non-perturbed plane. To correct for this bias, an iterative procedure is used: horizontal velocities derived from particle motions projected onto the mean water surface are used to construct a first estimate of the perturbed free surface topography. Particle motions are projected again onto the estimated perturbed topography, then a new horizontal velocity field is obtained, and so on. Since free surface slopes depart only moderately from horizontal, this procedure converges very quickly and two iterations were found sufficient.

5.4. Averaging of the velocity field

An averaged horizontal velocity field remains to be extracted from the full set of horizontal particle velocities defined at horizontal positions . Rather than combining right away the vectors corresponding to different times t into one single set to be indiscriminately binned and averaged, the regularisation procedure adopted here seeks to preserve the streamline structure which underlies the various individual velocity vectors. This structure presents important implication for the error formation process, and it is sought to account for it by splitting the averaging procedure into two steps: 1) along-trajectory averaging; 2) transverse averaging. For each of the two steps, a balance between attenuation and noise error can again be sought in order to minimise the combined error.

Introducing hats to distinguish measured estimates from hatless "true" quantities, consider measurements acquired at successive time instants tm along a single particle trajectory. We have the standard particle tracking velocimetry (PTV) estimate (e.g. Adrian 1991)

(31)

where measurement approximates the true particle position at time tm, and measurement approximates the true instantaneous velocity at intermediate time . A simple stochastic model for error formation can be derived by assuming (Wernet and Pline 1993) that measured positions correspond to the actual particle position to which is added Gaussian white noise, i.e.

(32)

where σx is the root mean square (rms) error on particle position, and where the γ m are independent Gaussian random variables of zero mean and unit standard deviation. Invoking (31) and (32) as well as Lagrangian relation x(t) =  u(t)dt, velocity estimate can be expressed in terms of the actual velocity history u(t) as

,(33)

where

.(34)

This constitutes an observation equation (Honerkamp 1998) for the PTV measurement process, i.e. an idealised description of the connection between the true signal and the measured one. From expression (33), two discrepancies are seen to arise between measurements and the actual instantaneous velocities u(t). The first discrepancy is a corrupting high frequency noise due to random position error. The second discrepancy results from averaging over interval Δt, which attenuates high frequency velocity variations.

Minimisation of the combined error can again be sought by balancing the effects of noise and attenuation. A filtered signal can be obtained by convoluting signal with a rectangular window of width (2n+1)Δt:

.(35)

(The same operation can also be applied to the transverse velocity measurements to yield along-trajectory filtered values .) Substitution of (35) into (33) yields an observation equation for the filtered velocities

(36)

where

.(37)

Relation (36) is equivalent to the earlier observation equation (33), except that it is applied to the longer time interval (2n+1)Δt. From (36), the root mean square (rms) error on the filtered velocity due to position error is given by

.(38)

Because measurements are obtained along a trajectory, position errors cancel each other along the way and only the endpoint errors matter. As a result, the rms error on velocity (38) due to position errors decreases in proportion to 1/(2n+1), a much faster decrease with filter width than the factor which would be obtained if errors on velocity were uncorrelated Gaussian random variables. It is to account for this feature of the error formation process that a two step procedure starting with along-trajectory filtering is selected in preference to indiscriminate binning of the whole set of velocity measurements. Because the present experiments feature variations of relatively long wavelength, one can approximate , yielding for the random position error the rough estimate

.(39)

The resulting values for the different runs are listed in Table 1. Converted back to pixel coordinates, the particle positioning error is around 0.1 pixel. Due to the limited resolution of the camera, such small errors in pixel units still translate into significant errors on velocities.

The corresponding attenuation error can be estimated by idealising the velocity history of a particle along its streamline as a harmonic signal, i.e.

(40)

where angular frequency ω ~ 2π[u]/[λ], amplitude , and characteristic scales [u], [λ], and [u'] can be estimated for a given run (see Table 1). Alternatively, equation (17) can be invoked to approximate . Convoluting signal (40) with a rectangular window of width (2n+1)Δt results in filtered signal

(41)

where the damped amplitude is given by (e.g. Jähne 1995)

(42)

in which again sinc(x) = sin(x)/x. Developing (42) by Taylor around the origin and setting , the rms attenuation error on the filtered signal can be approximated by

.(43)

Balancing position noise error (38) and attenuation error (43), one obtains for the conditions of the present experiments an optimal filter width (2+ 1) = 13. The estimated error on trajectory-filtered velocities due to either noise or attenuation is of the order of 0.01 m/s.

The last remaining operation is transverse averaging. This is performed in two steps. First, filtered along-trajectory velocities are sampled at the intersections of the trajectories with transverse planes placed at equal intervals x along the x-axis. This yields for each x a set of irregularly distributed data points . The second step of the procedure consists of converting these measurements into regularly spaced data at positions yj of a uniform grid. The averaging is performed by way of overlapping bins of equal width Δy, where the bin width Δy is again chosen to minimise the combined error resulting from noise and attenuation. Noise at this stage comes both from the imperfect PTV measurements and from fluctuations in time of the physical trajectories, and is assumed to be Gaussian. Using the same error balancing approach as before (see section 4.3), the optimal bin width is found to be Δy ≈ 35 mm. The estimated error due to either noise or attenuation at the transverse averaging stage is again of the order of 0.01 m/s. Noise and attenuation errors incurred at both the streamline and transverse averaging stages yield an overall rms velocity error of the order of 0.02 m/s. This amounts to errors of the order of 2 % relative to the mean velocity (see Table 1). Upon invoking the approximate Bernoulli equation (27), elevation errors due to velocity measurement errors are of the order of 2 mm.

To these must be added the neglected contribution from second order perturbations. Assuming the errors to be independent, an overall error is again estimated by adding the squares of the rms contributions and taking the square root. The resulting values for the various runs are given in Table 1. Reconstruction errors arising from interpolation onto a regular grid are not significant for the velocimetric method because interpolation is conducted to first order accuracy. The interpolants used are piecewise linear rather than piecewise constant. The latter approach, less accurate but more robust, is used for binning in the stereo method due to the presence of stereo mismatches.