Translational and rotational dynamics of a large buoyant sphere in turbulence

We report experimental measurements of the translational and rotational dynamics of a large buoyant sphere in isotropic turbulence. We introduce an efficient method to simultaneously determine the position and (absolute) orientation of a spherical body from visual observation. The method employs a minimization algorithm to obtain the orientation from the 2D projection of a specific pattern drawn onto the surface of the sphere. This has the advantages that it does not require a database of reference images, is easily scalable using parallel processing, and enables accurate absolute orientation reference. Analysis of the sphere's translational dynamics reveals clear differences between the streamwise and transverse directions. The translational auto-correlations and PDFs provide evidence for periodicity in the particle's dynamics even under turbulent conditions. The angular autocorrelations show weak periodicity. The angular accelerations exhibit wide tails, however without a directional dependence.


Introduction
Particles suspensions in turbulent flows can be found in a wide range of natural and industrial settings. The behavior of these particles depends on several parameters including their size and density ratio, the particle Reynolds number, the particle Stokes number and the Taylor Reynolds number of the carrier turbulent flow. Understanding how the dynamics of particles is influenced by these parameters is crucial to make predictions on global phenomena of interest such as pollutant * : v.mathai@utwente.nl † : chaosun@tsinghua.edu.cn transport, cloud formation, and mixing in industrial processes (Toschi and Bodenschatz [25]). In many natural settings, particles can have a large size, and their density can be different from the carrier fluid. Theoretical studies often model such objects as passive, finite-sized particles advected by the flow (Biferale et al. [2], Calzavarini et al. [4], Maxey and Riley [17]). These are applicable in the limit of vanishing particle and shear Reynolds numbers. However, in most practical situations, the density-mismatch of particles with the fluid results in finite drift velocities, which leads to finite particle Reynolds numbers (Jimenez [12], MacKenzie and Leggett [15], Skyllingstad et al. [24], Mathai et al. [16] etc). In such situations, particle dynamics can be strongly influenced by the unsteady wakeinduced forces. For instance, a large buoyant sphere freely rising through quiescent fluid displays rich variability in translational dynamics (Ern et al. [6], Horowitz and Williamson [9,10]). The forcing responsible for such varied dynamics is linked to the vorticity shed in the wake of the particle (Achenbach [1], Govardhan and Williamson [7]). Since these forces may not act along the geometric centre of the particle, it is possible that these could as well induce torques on the body. Little is known about the resulting translational and rotational dynamics, particularly for the case of buoyant particles in turbulence, which motivates us to develop a reliable measurement technique for studying these issues.
In three dimensions, an object's location and orientation can be fully described by six independent variables. Many physical experiments rely on image analysis to obtain these parameters from experimental data (Bovik [3], Meyer et al. [19]). Most of these systems capture translation and retrieve orientation from relative motion of translating nodes (Klein et al. [13]). In the more elementary forms, the translation of nodes could be used to determine the velocities of the particles. However, these methods were not accurate enough when higher derivatives of orientation were to be determined. Zimmermann et al. [31] introduced a method based on the identification of possible orientation candidates at each time step using projections of a pattern painted on the surface of a sphere. They found surprisingly intermittent behavior in the acceleration statistics of a neutrally buoyant sphere of diameter of the order of the integral scale.
In the present work, we introduce buoyancy to the problem of spherical particle dynamics in turbulence. We recover the translation and rotation as a function of time. The core of the method, which is to compare experimental images to synthetic ones is the same as proposed by Zimmermann et al. [31]. The novelty here lies in the way the pattern is generated (both physically on the surface of the particle and numerically for the comparison with actual images to match the orientation). Hence, the synthetic images for any given orientation are analytically known and do not need to be determined from static images. Furthermore, the method is easily scalable using parallel processing, which enables an accurate absolute orientation reference. These aspects are elucidated in section 3. In addition, we describe a smoothing spline based roughness limiting technique that enables accurate representation of the higher derivatives of experimental data. Finally, in section 5, we present the main results of our investigation on translational and rotational dynamics.

Experimental Setup
The experiments were conducted in the Twente Water Tunnel facility (TWT), designed to study particle-laden flows (see figure 1 and Rensen et al. [23]). The measurement section has dimensions 0.45 × 0.45 × 2 m 3 , with three glass walls providing optical access to perform particle tracking experiments. The setup houses an active grid above the measurement section, consisting of 24 independently rotating motors, which produce nearly homogeneous and isotropic turbulence with Re λ up to 300 in the downstream section of the water tunnel. The flow in the measurement section was characterized using a cylindrical hot film probe (Dantec 55R11) by following the same methodology as reported in Mercado et al. [18]. The experiment reported here was performed at Re λ ≈ 300. The dissipation rate = 505 × 10 −6 m 2 /s 3 , and the dissipation length and times scales were approximately 211 µm and 44 ms respectively. The sphere used in the study has a diameter, d p = 25 mm, with an effective mass ratio, m * ≈ 0.82, where m* is the ratio of mass of sphere to mass of the sphere's volume of water (see Govardhan and Williamson [7]). The mass ratio was chosen such that the mean rise velocity of the sphere matched the mean downward flow velocity in the measurement section. This was necessary for obtaining sufficiently long sphere trajectories for well converged Lagrangian statistics. The spheres were designed as spherical shells with the MP300 resin, with a bulk density ≈ 1089 kg/m 3 . The angular inertia of the sphere was comparable to that of a homogeneous sphere of water, as was the case in the study by Zimmermann et al. [31]. The spherical shells were made using a 3D printing technique, and the surface roughness was within 50 microns. The printing resolution depended on two factors. First, the resolution of the 3D printer used for making the stencil. In the present case, we could produce complex stencil designs with dimensions as small as a millimetre on a 25 mm sphere. The second limiting factor was the painting procedure itself. Once the stencil was fixed on the sphere, the pattern had to be spray painted. Here again, it was practical to spray paint patterns of 2 mm smaller dimension.
The recordings were made with two Photron PCI-1024 high-speed cameras at 500 fps and megapixel resolution (see figure 2(a) & (b)). The cameras were positioned at a 90 • angle between them and focused at the center of the test section on a 150×150 mm 2 area, which resulted in a spatial resolution of 150 µm/pixel. The images showed that perspective effects were negligible. This was done by placing a sphere of known orientation at the corners of the field of view. The retrieved orientation varied less than 3 • . The measurement volume was lit by eight 20 W LED lamps from the sides. The flow velocity was tuned to ensure that the spheres stayed in the viewing window for considerable duration.

Method
The position and orientation of the sphere are determined using image analysis methods. The background was painted gray for contrast between the black and white pattern on the sphere. The sphere is separated from the background by subtracting the absolute difference in intensity from the background image. Subtracting the absolute difference yields a full dark circle with The accuracy of CHT varied with angle of view. Therefore, choosing a good pattern is one of the most important steps. We found that a pattern that contained almost equal fraction of dark and bright areas for most angles of view improved the circle detection accuracy. Additionally, we used the outputs from previous frames to improve and speed up the detection process. The CHT also returns the detected circle diameter D c in pixels from the image. Since the sphere diameter (25 mm) is known, we use a calibration (25/D c ) mm/pixels for the recorded images. This ensures a correction for the minor magnification changes due to the sphere moving forward or backward in the measurement section.
There are several approaches to obtaining the three dimensional trajectory of a particle. One approach is to reconstruct the spatial positions using multiple cameras (Ouellette et al. [21]). Precise spatial reconstruction is needed when there are multiple particles in the image and identifying and separating particles can be a challenge. In the present work, there are only a few particles in the measurement volume. We use a different method using two orthogonal cameras to obtain the three dimensional trajectory. The camera arrangement is such that the magnification and field of view are comparable between the two cameras. The camera field of view covered the full width of the water tunnel. We first obtain two dimensional trajectories of particles from both cameras. The redundant data corresponding to the vertical motion for two cameras is used to compare individual trajectories. In order to avoid ambiguities, we cross-correlate the vertical acceleration time series from the two cameras. A match is said to be found when the cross correlation of two acceleration time series yields a coefficient greater than 98 percent, even for long trajectories such as the one shown in figure 2(b). This yields a simple yet robust three dimensional track of the particles without needing to resort to complex spatial reconstruction algorithms.
Obtaining the rotation of the sphere is a more complex task compared to position tracking. Hence, the remainder of this section will be used to describe the method that is used to obtain the rotation. The method to determine the rotation can be divided into four parts. Initially, a suitable boolean surface pattern is created. This pattern can be described as a piecewise constant analytic function F (θ, φ) such that, given a coordinate on the surface of the sphere, the function returns either a zero or a one depending on the color of its corresponding infinitesimally small surface element. Here, θ and φ are the azimuthal and the polar angles respectively, and F (θ, φ) is radius independent. Second, the pattern is drawn onto a physical sphere. This is realized using a 3D-printed painting stencil and an airbrush system. It is imperative that F (θ, φ) is painted as accurate as possible to decrease the introduced error in this step. Third, a synthetic 2D image is constructed from a projection of the surface of the sphere onto a plane. The synthetic image for any given orientation is analytically known and does not need to be determined from static images. The projection is a function of the angle of rotation of the sphere and can be conceptually understood as the analog of the projection of the physical sphere on the record-ing camera. Finally, the rotated and projected synthetic pattern is compared to an image of the physical sphere. The minimizing function can take different forms. For instance, a cross-correlation function between the synthetic pattern and the camera image could be used to find the best match. Alternately, a suitable cost function may be used to search for a match. Here we use a cost function, defined as the sum of the absolute difference between the binarized image pixels and the corresponding pixels in the synthetic image. The orientation for which this comparison yields the best match is then determined using a Nelder-Mead minimization algorithm. These steps are illustrated in a flowchart in figure 3. The Fortran90 code and the stl-file of the stencil used for painting the spheres have been included as supplemental material. The choice of painted pattern is an important step in the method of orientation detection. A necessary condition is that the projection of F (θ, φ) onto a plane must be unique for each orientation of the sphere. Additionally, it is desirable that the pattern contains a minimum of edges and corners. The latter criterion is essential for fast con-vergence of the algorithm to the global minimum. Here we use the axis-angle method to describe the orientation considering its straightforward and singularity-free definition. It allows for smooth and continuous rotation from any orientation including around the Euler-angle singularities as it does not suffer from gimbal lock problems.
The pattern F (θ, φ) represents a simply connected region, even then, some local minima may still arise due to the two-tone color limit. When the orientation is unknown, for instance in the first frame of a movie, a comparison of several initial orientation estimates solves this local minima problem. This yields the global minima. However, every initial estimate reduces the performance. Fortunately, given a sufficiently high ratio of frame rate over rotation rate, any frame in a sequence may use the outcome of the previous frame as initial estimate, reducing the number of initial estimates required to just one. This causal method is the preferred method for analyzing large sequential datasets. The computational performance of the method is approximately 30 frames per second on a contemporary computer, suggesting that processing of large datasets is straightforward.
To determine the numerical accuracy of the method, several sets of 1024 projections are created using pseudorandom generated orientations. These synthetic images are used as input to our algorithm. A probability density function of the difference between the actual orientation and the orientation as found by the algorithm is shown in figures 5(a) & (b), and the corresponding scaling of the standard deviation as function of the image width N (in pixels) is shown in figure 6. The standard-deviation of accuracy scales as O(N −2 ), and assuming a Gaussian error distribution, this means that less than 1% of the measured data showed an error larger than 1 degree for N ≥ 50. Also, figure 6 shows that noise decreases the accuracy to some extent but does not affect the reliability of the algorithm. In the experiment, most image artifacts arise due to shadows and glare, which may be reduced by using diffuse light-sources and a matte-paint finish on the surface of the sphere.
The axis-angle output, ρ ≡ (k x , k y , k z , α), is defined with respect to a reference orientation in the camera coordinate system (see leftmost image in figure 4 (a)). A typical output obtained from a buoyant sphere in a turbulent flow is given in figure 9(a). This output has little physical significance. Quantities of greater relevance to a rotating sphere are its rotational kinetic energy Iω 2 and the net torque exerted on it by the surrounding fluid Iα, where I is the moment of inertia of the sphere, and ω and α are the angular velocity and acceleration respectively. We adopt the following approach to compute ω and α in the lab coordinate system. Firstly, a high framing rate is used for accurate estimation of the time derivatives. In the present case, a frame rate of 500 fps ensures about 30 recordings in one Kolmogorov time. Therefore, within the inter-frame time interval ∆t = 1/500 sec, the particle's angular velocity may be assumed constant.
According to Euler's theorem, the angular velocity is linked to the inter-frame axis-angle, ρ ∆ ≡ (k ∆x , k ∆y , k ∆z , ∆α), by the relation Now ρ ∆ may be obtained from the reference orientation based axis-angle output (ρ) through a few transformations. Firstly, the axis angle output is converted to rotation matrix using Rodrigues formula. Let [R 0,i ] and [R 0,(i+1) ] represent the rotation matrices corresponding to rotations from the reference orientation to the orientations in i th and (i + 1) th images respectively. Then This follows from the identity that [R] −1 = [R] T for rotation matrices. From [R i,(i+1) ], ρ ∆ ≡ ρ i,(i+1) may be obtained, and eq. 1 gives the angular velocity ω(t) in the lab coordinate system.
The accuracy of the detection in a real experiment needs to be verified by other methods. We use a twocamera arrangement for this. If the orientations determined by the two cameras are comparable, then the method can be regarded accurate. Figure 9(b) shows the three orientation angles obtained from Camera 1 for a long particle trajectory. The same sphere was captured by Camera 2, which was placed at a 90 • angle with Camera 1. The difference in orientation prediction is plotted in figure 9(c). The maximum deviation is within 2.5 degrees. These differences are due to experimental noise and may be filtered out in the data smoothing step, to be explained in the following section. Figure 7: Roughness vs error frontier for a family of candidate smoothing splines. The kink marks the critical error tolerance of the optimal smoothing spline. Inset shows log(C(s)) vs. log(E(s)). The optimal fit corresponds to the maximum curvature point indicated by the red star symbol.

Evaluating higher derivatives from experimental data
Evaluating the forces and torques is an important step in understanding any dynamical system. One way to do this is by using accelerometers. However, this may not always be convenient since a direct measurement of acceleration usually requires attaching devices to the body. Moreover, the method is intrusive in nature and often leads to variations in the mass and center of gravity of the system (Zimmermann et al. [32]). In such situations, it may be suitable to numerically compute the acceleration from the second derivative of position of the body. This approach is the basis of our present study, where we determine both the position and orientation of a large buoyant sphere in turbulence using recordings from a high speed camera. Estimating the derivatives from experimentally determined position and orientation data is a nontrivial task. The difficulty arises because data obtained from experiments will have inherent noise due to measurement uncertainties. Historically, there have been two popular methods for smoothing particle trajectories in turbulent flows. One method involves fitting parts of the particle trajectory to a polynomial of second order or higher. Voth et al. [29] used a second-order polynomial, Lüthi et al. [14] and Mercado et al. [18] used a third-order polynomial. Other researchers have used a Gaussian kernel for smoothing (Mordant et al. [20], Volk et al. [28]). Both the mentioned methods employ piecewise discontinuous fitting of analytical functions to smooth out the noise. The effectiveness of these methods in filtering out the experimental noise depends on the fitting parameters chosen. In spite of the extensive literature on smoothing methods, there still prevails a general lack of consensus on the method of finding the optimal fitting parameter. Existing guidelines for the choice of the fitting parameters, such as those suggested by Mordant et al. [20], require knowledge of the smallest time scales in the flow apriori, and therefore are likely to introduce a bias into the analysis.
In this paper, we explore an alternate method using smoothing splines to filter out the experimental noise. The method is based roughly on the work by Epps [5]. While this was originally used to estimate the translational accelerations experienced by a water-entering object (Truscott et al. [26], Truscott [27]), the method may be applied to particle trajectories in turbulence. Consider a general set of experimental data y(t) acquired at high temporal resolution. We use a smoothing spline based roughness limiting method to reduce the experimental noise in the data for obtaining higher derivatives. The function spaps in Matlab is defined by two fitting parameters: the order of the smoothing spline and the error tolerance, E(s) = |y(t) − s(t)| 2 dt. We use a quintic smoothing spline, which ensures that the second derivatives are properly represented. This leaves us with one fitting parameter, the error tolerance. We use a roughness estimate, R(s) = t N t1 | d 3 s dt 3 | 2 dt, to scan for the most suitable fit to the experimental data. This is based on the assumption that the true function does not have very large changes in acceleration, which typically are due to noise in the data. In figure 7, we demonstrate the step to determine the optimal fit from experimental data. Increasing the error tolerance beyond E crit does not reduce the roughness of the curve. For error tolerances below E crit , the roughness of the fits are increased significantly, indicating that the noise contained in the data is not properly removed. Therefore, the fit corresponding to the kink can be thought to best represent the true curve.
We test the sensitivity of the velocity and acceleration to the choice of the spline. For this, we use a sample analytical function defined as a combination of sine functions, y = a sin x + b sin 2x + c sin 4x + d sin 10x, where a, b, c and d were generated randomly in the [0 1] range. We introduce random noise to the signal with signalnoise-ratio ≈ 4. The true function and noisy input are shown in figure 8(a). In figures 8(c)-(h), we compare the velocity and acceleration estimates for critical (E crit ), sub-critical (E 1 ) and super-critical (E 2 ) error tolerances. The sub-critical error tolerance (figure 8(f)) yields very high accelerations, while the super-critical error tolerance over-smoothes the curve ( figure 8(h)). Clearly, the acceleration estimate for critical error tolerance was obtained without prior knowledge of the timescales of the  Another interesting feature of the present method is that the kink (figures 7 & 8(b)) flattens out for experimental data with low signal-noise-ratios. Thus, the method also serves as a check for the quality of data, which can be difficult to estimate for standard smoothing methods. Alternately, the spline fit guess for the particle acceleration from a few representative tracks could be used to determined the optimal parameters for Gaussian Kernel smoothing (Ouellette [22]). This can speed-up the fitting process when sampling large datasets, while ensuring that the fitting windows are properly chosen.

Results and discussion
We present results on the translational and rotational motion of a large buoyant sphere in nearly homogeneous and isotropic turbulence. We first address the question of how the particle's translational velocity and acceleration decorrelate in time. In Figure 10 (a), we plot the Lagrangian autocorrelation function for the horizontal (x & z directions) velocity and accelerations. The particle response scale is fairly well predicted from the relation τ v = d p /(St × U r ), where St−Strouhal number ≈ 0.2, and U r is the measured mean rise velocity relative to the flow. Interestingly, both velocity and acceleration decorrelate in the same time ∼ 0.5 sec, and they both display periodicity. This decorrelation behavior is fundamentally different from that of a passive finite-size particle in turbulence, for which the acceleration decorrelates in much shorter time than the velocity (Ishihara et al. [11]). Therefore, in the present case, the dominant velocities and accelerations originate from vortices with  the same timescale. Figure 10(b) shows the angular velocities and accelerations in the horizontal direction. The time axis is normalised with the particle-sized eddy time scale, τ d = (d 2 p / ) 1/3 ∼ 1 sec. The angular acceleration crosses the first minima in approximately 1 τ d , which is four times longer compared to the minima crossing time for translational acceleration. The angular velocity components decorrelate even slower than the angular acceleration, and they show only a weak periodicity. The angular velocities and translational accelerations are weakly correlated, with a correlation coefficient ≈ 0.05. Surprisingly, we do not find any preferential orientation of the angular velocity vector with the translational acceleration vector, as found by Zimmermann et al. [30]. Therefore, in the present case, particle rotation does not appear synchronized with its translational motion, and we find no strong evidence suggesting a lift force. More detailed studies are to be conducted to gain further insights into the underlying physics. In future, we also aim to look into the flow structure around the sphere along with its motion.
In figure 11(a) & (b), we show the probability density functions (PDFs) of the horizontal and vertical components of the velocity and acceleration. The horizontal velocity PDF shows a symmetric flat-head distribution with sub-Gaussian tails. The behaviour is notably dif-ferent from the nearly Gaussian velocity PDFs found for neutrally buoyant particles (Mercado et al. [18]). The flat top of the PDF may be explained from a typical time series of velocity of the sphere. We observe that the particle undergoes repeated cyclic motions, weakly disturbed by the carrier flow. The velocities in these cycles lie in the ± 1.5 a 2 1/2 range. More extreme accelerations are less frequent and contribute to the low probability tails of the PDFs. The horizontal acceleration PDF is Gaussian. The vertical velocity and accelerations have positively skewed distributions, clearly indicating a directional dependence. This strong anisotropy is not inherent in the carrier turbulent flow. The water tunnel flow was reported to be nearly isotropic in earlier studies (Mercado et al. [18]). Therefore, the anisotropy here is expected to originate from the up-down asymmetry set up by the vortex-shedding downstream of the sphere. The PDFs of the rotational quantities reveal a different story (Figure 11(c) & (d)). Both horizontal and vertical angular velocities follow a nearly gaussian distribution. The angular acceleration PDFs show symmetric distributions but with wide tails compared to the translational acceleration PDFs. It may be noted that Re λ and particle size ratio Ξ = d p /η are comparable to the experiment of Zimmermann et al. [31]. Thus, changing only the particle-fluid density ratio brings in these new effects, leading to observable differences in particle dynamics in turbulence. Future experiments will be aimed at tracking simultaneously the particle and the flow around it.

Conclusions
We have conducted experimental measurements of the dynamics of a large buoyant sphere in a turbulent flow and described the methods and procedures to track simultaneously its position and orientation. A fast, accurate and adaptable method that determines the absolute orientation of the sphere in 3D space is described and validated. The standard deviation of error for this method is approximately 0.2 degrees. Distortion using salt & pepper noise with SNR= 2 increases the standard deviation of error to 0.7 degrees, however without impacting the reliability of the method.
We followed the approach of using the experimentally determined data to calculate time derivatives (velocity and acceleration) of position and orientation. To this end, we employed a smoothing spline based roughness limiting technique, which enables an accurate estimation of higher derivatives. The method has the advantage that it yields an optimal fit that is continuous and differentiable.
Our results on the velocity and acceleration statistics reveal that buoyant spheres have very different dynamics from the well-explored class of neutrally buoyant particles in turbulence (Zimmermann et al. [31], Toschi and Bodenschatz [25], and Homann and Bec [8]). We detect the influence of trailing wake, resulting in periodicity in the Lagrangian auto-correlations and anisotropy in the translational dynamics. The rotational velocity and acceleration PDFs show wide tails, however without any observable skewness in the streamwise and transverse directions. The present measurements provide clues about how a buoyant particle interacts with a turbulent flow. The methods presented here open up a new direction in the exploration of particle dynamics in turbulence.