Measurements of the Solid-body Rotation of Anisotropic Particles in 3D Turbulence

We introduce a new method to measure Lagrangian vorticity and the rotational dynamics of anisotropic particles in a turbulent fluid flow. We use 3D printing technology to fabricate crosses (two perpendicular rods) and jacks (three mutually perpendicular rods). Time-resolved measurements of their orientation and solid-body rotation rate are obtained from stereoscopic video images of their motion in a turbulent flow between oscillating grids with $R_\lambda$=$91$. The advected particles have a largest dimension of 6 times the Kolmogorov length, making them a good approximation to anisotropic tracer particles. Crosses rotate like disks and jacks rotate like spheres, so these measurements, combined with previous measurements of tracer rods, allow experimental study of ellipsoids across the full range of aspect ratios. The measured mean square tumbling rate, $\langle \dot{p}_i \dot{p}_i \rangle$, confirms previous direct numerical simulations that indicate that disks tumble much more rapidly than rods. Measurements of the alignment of crosses with the direction of the solid-body rotation rate vector provide the first direct observation of the alignment of anisotropic particles by the velocity gradients of the flow.


Introduction
The motion of particulate matter in fluid flows has long been a central problem in both fundamental and applied fluid mechanics. The dynamics of spherical particles was a natural starting point for Stokes, who considered this problem in the 19th century [1]. A long list of researchers continued work to describe the motion of spherical particles at low Reynolds number, leading to the Maxey-Riley-Gatignol equations [2,3]. The development of numerical simulations [4,5] and experimental tools [6,7,8] capable of revealing the motion of spherical particles in complex flows has led to a recent resurgence of work on particle dynamics [9].
Extending the problem to non-spherical particles in general is important for many applications ranging from icy clouds [10], to bio-locomotion [11], to suspension flows in industrial settings [12]. In 1922, Jeffery analyzed the motion of axisymmetric ellipsoidal particles in Stokes flow [13]. Their tumbling rate can be expressed as: where p is a unit vector along the symmetry axis of the ellipsoid, α is the aspect ratio, Ω is the anti-symmetric part of the velocity gradient tensor (vorticity), and S is the symmetric part of the velocity gradient tensor (strain rate). A wide range of studies has explored deviations from Jeffery dynamics due to particle inertia, particle shape, and Reynolds number, among other factors [14]. Analytic work [15,16] and numerical simulations [17,18,19,20] have made significant progress extending studies of non-spherical particles in complex flows, but experimental measurements of the dynamics of anisotropic particles have lagged far behind due to the difficulty of measuring their time-dependent orientation in three dimensions. Recently, methods have been developed for measuring time-resolved orientation and position of thin rods in 3D turbulence with stereoscopic optical imaging [21]. Another technique uses large transparent anisotropic particles with tracer particles inside and measures their rotations with particle image velocimetry [22]. In this paper, we introduce a new way to measure the orientations and rotational dynamics of anisotropic particles that behave like ellipsoids. Bretherton has shown that many anisotropic particles have equivalent ellipsoids so that their tumbling rate follows Equation 1 with an effective aspect ratio [23]. We have identified particles of this class whose orientations can also be directly measured from stereoscopic imaging. Three-dimensional printing allows us to create anisotropic particles made of mutually perpendicular thin rods. Two perpendicular thin rods form a cross, while three perpendicular rods form a jack. Arguments in Bretherton [23] and resistive force theory calculations in Appendix A show that a cross rotates like a disk, an ellipsoid with aspect ratio α 1. Similarly, a jack rotates like a sphere, which is an ellipsoid with α = 1. We have used stereoscopic video imaging to directly measure the orientation of these objects as a function of time. The methods developed in this paper thus allow us to obtain Lagrangian measurements of the full solid-body rotation rate.
Since a jack rotates like a sphere, its solid-body rotation rate couples only to the vorticity. Thus, measurements of the rotation of neutrally buoyant Kolmogorov-scale jacks are direct Lagrangian vorticity measurements-a long sought goal of fundamental fluid mechanics. Several other methods for measuring Lagrangian vorticity have also been developed. For example, the vorticity optical probe measures the reflections from planar mirrors embedded in spherical particles [24]. Stereoscopic imaging has been used to measure the rotation of large spheres by applying patterns to the sphere surface [25,26] or embedding small fluorescent tracers in transparent particles [27]. Our technique of measuring the solid-body rotation of small jacks allows us to accurately measure the fluid vorticity using straightforward imaging methods.

Printing 3D Particles
We use 3D printing to fabricate anisotropic particles in the shapes of crosses and jacks in order to measure the dynamics of axisymmetric ellipsoids across the range of aspect ratios. To print at both high resolution and in high quantity, we used a Connex 500 to make 10,000 of each particle shape. Arm lengths were 3 mm, which is 6 times the Kolmogorov length scale in our flow. The diameter of any given arm on a rod, cross, or jack is 300µm, the smallest we could achieve while maintaining the structural integrity of the particles.
To print particles with such small cylindrical arms, it was important to ensure that none of the arms were along the build axes. Arms that lay along the vertical build axis had defects and would often break off, while arms lying in the horizontal plane tended to flatten. Another important difficulty in printing O(10 4 ) particles is removing the support material. Connex printers use a different material for the support structure and the printed particles. The support material can be partially dissolved using a strong base solution (e.g., NaOH) without affecting the particles themselves. We found that using an ultrasonic bath made the removal process much more efficient, and the particles could be filtered out of the solution with almost no loss.
In order to have neutrally buoyant particles, the density of our fluid is matched to the particle density. The density is adjusted by adding calcium chloride to water [21]. We used the print material VeroClear, whose bulk density was quoted at 1.17 g/cm 3 . However, we found that the manufacturer's quote differed significantly from the density of the fluid in which the particles were neutrally buoyant. After the particles were immersed in the fluid for several hours, we found that different particles were neutrally buoyant at slightly different fluid densities ranging from 1.21 to 1.23 g/cm 3 . We chose the density of the fluid for the experiment as 1.22 g/cm 3 based on the population average. More work is needed to understand the mass density distribution inside 3D printed objects, but our particles are sufficiently density matched that their rotations should accurately represent the neutrally buoyant case [28]. To make the particles fluorescent, they were placed in a high concentration Rhodamine solution at elevated temperature (60-80 • C) for several hours.

Turbulent flow between oscillating grids
The turbulence is generated in a 1 × 1 × 1.5 m 3 octagonal prism using two parallel 8 cm-mesh grids oscillating in phase (see Figure 1). The measurements presented here were performed at a grid frequency of 1 Hz, which produced a flow with a root mean squared (rms) velocity U ≡ u i u i /3 = 20 mm/s and energy injection length scale, L = U 3 /ε = 60 mm so that R λ ≡ (15U L/ν) 1/2 = 91. The energy dissipation rate, ε = 133 mm 2 s −3 , was measured from the mean square tumbling rate of jacks as described  [29]). In the octagonal flow between oscillating grids, a central viewing volume in the focus of the four stereoscopically arranged video cameras is illuminated by a green Nd:YAG laser.
in section 3. The kinematic viscosity of the CaCl 2 solution that is approximately density matched to the particles is ν = 2.17 mm 2 /s. The Kolmogorov length scale is η = (ν 3 /ε) 1/4 = 526 µm, and the Kolmogorov time scale is τ η = (ν/ε) 1/2 = 128 ms. We chose this low Reynolds number to ensure that our particles were only 6η, where the mean square rotation shows only a small deviation from the tracer limit [17,30].
We use four stereoscopic cameras (1280×1024 at 450 Hz) with a custom real-time image compression system [31,29] to allow continuous imaging. A green Nd:YAG laser with 50 W average power illuminates a detection volume of roughly 3 × 3 × 3 cm 3 at the center of the tank, where the flow is quite homogeneous [32]. We used two perpendicular expanded beams in the horizontal plane, each reflected back upon itself, in order to provide 4 directions of illumination and minimize shadowing of some arms of a particle by other arms. With an average particle number density of only 5×10 −3 cm −3 , there was a particle in view less than 20% of the time, which made the image compression system essential for acquiring enough trajectories to converge statistics [31].

Image Analysis.
Orientation.-An example of a jack imaged by all four cameras is shown in Figure 2. When a particle is visible on all four cameras, we find the three dimensional position of the particle using stereomatching methods [33]. We developed a nonlinear fitting algorithm to determine the orientation of a particle from a set of stereoscopic images. Any orientation is specified by a rotation matrix, O, which can be parameterized by three Euler angles (φ, θ, ψ) [34]. From a measurement of the particle's center and an initial guess of its orientation, we construct a model of the particle and project it onto each of the four image planes using the calibrated camera parameters [35]. The total difference in intensity between a model image and an experimental image provides the residual that is minimized by a nonlinear search in Euler angle space.
To minimize computational cost, we project the endpoints of each arm onto the image plane of each camera, and then model the intensity distributions in two dimensions. The model image of an arm is formed by a Gaussian intensity distribution along the width of each arm (in 2D) and a Fermi intensity distribution across the length of the arm. For jacks and crosses, each arm in the model has an identical intensity distribution. As can be seen from Figure 2, the arms in the experimental images are not of uniform intensity. The observed intensity has a non-trivial dependence on the angles between the arms, the illumination, and the viewing direction [36], which has not yet been included in our model. However, we found the simple model adequate enough for our purposes.
The orientation-finding algorithm must be seeded with an initial guess for the orientation. Except for the first frame of a track, we use the previous frame as the initial guess, since the rotation between frames is less than 0.01 radians. The initial guess for the beginning of a particle track is reliably generated using an optical tomographic reconstruction algorithm [37]. We also compared tomographic reconstruction as a method for finding particle orientation, but we found it to be both less accurate and more computationally expensive than our own algorithm.
The Euler angles found for a jack give 1 of the 24 orientations related by symmetry. To see this, consider 1 of the 6 arms of a jack and define it by the vector p =ẑ; there are 4 symmetrically identical orientations obtained by rotations of π/2 about the z axis. There are 4 such orientations for each of the 6 arms, for a total for 24. A cross has 8 symmetrically identical orientations, and a rod has 2. We ensure that we have a The two different color sheets trace out the path of the particle through space and time. The length of the particle track is 336 frames, or 5.7 τ η . A cross is shown every 15 frames. (b) A reconstructed trajectory of a jack in three dimensional turbulence. The three different colors distinguish the arms of the jack and trace out their path as the particle rotates. The dark green line denotes the trajectory of the jack's center. The length of the particle track is 1025 frames, or 17.5 τ η . A jack is shown every 50 frames. (Note: neither the crosses nor the jacks shown above are drawn to scale.) consistent series of orientations along a single trajectory by comparing the orientation found in each frame with that from the previous frame, and choose the orientation of the particle that produces the smallest total rotation between frames. Figure 3(a) shows a representative track for a cross, which is 5.7 τ η long. It demonstrates the effectiveness of our algorithm at determining the full range of orientations. At various points along the track, the rotation of the cross about a vector nearly coplanar with its two arms is clearly visible. In Figure 3(b) we also show an example jack trajectory that is 17.5 τ η long.
Solid-body rotation rate.-Once we have measurements of orientation and position at each time step, the natural quantity to consider is the rate of change of the unit vector defining the particle's orientation, which we call the tumbling rate,ṗ. For crosses and jacks, we can also measure the full solid-body rotation rate vector, ω s , from a series of orientation measurements smoothed along the particle's trajectory. One method for doing this has been described in [26]. We take a different approach using the tools we already developed for least squares optimization in Euler angle space. The problem can be framed as finding the initial orientation matrix, O(t i ) and the rotation matrix over a single time step, R, that together give the particle orientation matrix as a function of time, where τ f is the period between images. A non-linear least squares fit is used to find the Euler angles for the matrices O(t i ) and R that best match the measured orientation matrices. Then, the rotation matrix R can be decomposed as a rotation by an angle Φ about the solid body rotation axis,ω s , in accordance with Euler's theorem [34], from which we obtain the magnitude of the solid body rotation rate, ω s = Φ/τ f . The solid-body rotation rate, ω s , is related to the tumbling rate byṗ = ω s × p. The difference between the two quantities is thatṗ does not depend on the vector component of ω s lying along p. We use this relationship to determine the tumbling rate from measurements of particle orientation and their solid-body rotation rate. In order to correct for the contributions from orientation measurement errors to a measurement of the mean square tumbling rate, one needs to measureṗ over a a range of fit-lengths, τ = t − t i . Figure 4 shows our measurement of ṗ iṗi as a function of τ . The solid lines are fits of the function, where the first term models the random error that dominates at small τ and the second term models the effect of filtering experimental measurements of orientation [28]. The fit parameter C gives an estimate of the function value at τ = 0 if there were no random errors.

Results
From the measured solid-body rotation rates of many jacks and crosses, we can obtain the probability density function (PDF) of the squared tumbling rate, which is shown in Figure 5. Also shown is the PDF obtained from direct numerical simulations of spheres [21]. All three PDFs agree within experimental uncertainties. Numerical work has shown that there should be slightly more probability density in the tail of the PDF  Figure 4: Measurements of the mean square tumbling rate, ṗ iṗi , as a function of the fit length for both jacks (blue) and crosses (red). To determine the true value, we extrapolate to zero fit-length by fitting the data to a stretched exponential [28].
for rods and disks than there is for spheres [21]. For large tumbling rates, our data shows the opposite, with jacks having slightly higher probability than crosses, although the sources contributing to this small discrepancy are known. The error bars shown in Figure 5 account for random error as well as the systematic error that results from the fit-length dependence of the tumbling rate measurements. An additional source of error that is not included in the error bars is the self-shadowing of particles as they pass through certain orientations dependent on the camera configuration. The most dramatic cases are when an entire arm is missing from each of the four cameras, such that a jack, for example, will look like a cross along part of a given trajectory. The reduction in accuracy in determining these particular orientations occasionally leads to erroneously high measurements of the solid-body rotation rate, which pushes additional probability density towards the tail of the PDF. This effect is stronger for jacks than for crosses, which is why the jack PDF is slightly higher in our measurements at large tumbling rates. Previous work by Parsa et al. showed that the mean square tumbling rate for axisymmetric ellipsoids, p i p i , is strongly affected by alignment of anisotropic particles by turbulence [21]. They found agreement between simulations and experimental measurements of rods, but particles with a wider range of aspect ratios could not be measured experimentally with the tools available at that time. We have now measured the mean square tumbling rate of crosses and jacks using the techniques described in section 2. In Figure 6, we show the mean square tumbling rate normalized by the Kolmogorov timescale as a function of aspect ratio. Our measurements show crosses tumbling at a considerably higher rate than jacks or rods, but still more slowly than the randomly-oriented prediction. There is good agreement with simulations across the full range of aspect ratios.
In turbulence experiments, it is always a challenge to measure the energy dissipation rate, which appears in the normalization of the vertical axis in Figure 6. We attempted to make independent measurements using non-fluorescent tracer particles, but did not succeed because of reduced light scattering from tracers in the densitymatched CaCl 2 solution. However, the measurements of jack rotations provide a new way to measure the energy dissipation rate. Because they rotate like spheres, they give a direct measurement of the vorticity, and in isotropic turbulence their vorticity is directly related to the energy dissipation through Ω ij Ω ij = S ij S ij = ε/ν. This implies that for spheres ṗ iṗi = 1 6 ε/ν [21]. We use this to determine our energy dissipation rate. This makes our jack data at α = 1 match the simulations by definition, and the agreement of the cross data at aspect ratio α = 0.1 with numerical simulations is an independent result of the measurement.
The solid green curve in Figure 6 gives the mean square tumbling rate as a function of aspect ratio for randomly oriented ellipsoids. While the tumbling rate of jacks is unmodified from the randomly oriented case, both rods and crosses show Figure 6: Mean square tumbling rate as a function of aspect ratio. Solid squares are our data for crosses and jacks. The other data is from [21]. Plus symbols are numerical simulations. The triangle and circle are experimental measurements of rods at R λ = 214 and 160. The solid line is an analytic result for randomly oriented ellipsoids. smaller tumbling rates due to the effects of alignment by turbulence. Shin and Koch were the first to notice that the tumbling rate of rods that have correlated with a turbulent flow is reduced in comparison to that of randomly oriented rods [17]. More recent numerical work showed that the effects of alignment persist across the full range of aspect ratios [21]. The leading-order effects are contained in the Lagrangian threepoint correlations of the velocity gradients, which imply higher tumbling rates for disks than for rods [16]. A number of studies on the dynamics of rods show that they align with the vorticity, and this has been used to explain the slower tumbling rate of rods since the component of the vorticity along the rod axis does not contribute to its tumbling rate [19,20,38,39]. We find that alignment with vorticity is also responsible for the reduction of the tumbling rates for crosses.
We can directly measure the preferential alignment of a particle by measuring the angle between the orientation of the particle and its solid body rotation rate vector. In Figure 7, we plot the PDF of the magnitude of the cosine of this angle, |p ·ω s |, for both crosses and jacks, and compare them with numerical simulations. The peak near |p ·ω s | = 0 for crosses in Figure 7(a) shows that disks preferentially align with p perpendicular to ω s . Figure 7(b) confirms the same story as Figure 7(a). It shows that an arm of a cross is preferentially aligned with the solid body rotation rate vector. The height of the peak in Figure 7(b) is lower than that in Figure 7(a) because any vector in the plane of the disk is equally likely to align with ω s . In a turbulent flow, disks rotate like a spun coin rather than a tossed frisbee. For jacks in Figure 7(c), there is The observed alignment in Figure 7 seems natural if it is thought of as a result of the Lagrangian stretching of the fluid. Stretching will align both the vorticity and a long axis of a particle with the stretching direction [39]. For rods this means that p aligns with the vorticity, which decreases the tumbling rate. Because disks have their long axis perpendicular to p, they preferentially have p perpendicular to the vorticity, which makes p also preferentially perpendicular to the solid-body rotation rate. This creates a larger tumbling rate for disks sinceṗ = ω s × p, consistent with the data in Figure 6.
Although disks tumble faster than rods, they still tumble slightly slower than if they were randomly oriented, as seen in Figure 6. Since disks have a much smaller difference from the randomly oriented case, one might conclude that disks are less strongly aligned than rods by the turbulence. However, this is largely a result of the wayṗ is defined. We have measured the normalized mean square tumbling rate of a unit vector along one of the arms of a cross as 0.12 ± 0.02, which is smaller than the tumbling of either crosses or jacks and the same as measurements for rods within experimental error. This explicitly indicates that crosses (disks) are also strongly aligned by turbulence, and it is only the definition of p as perpendicular to the two arms that increases the mean square tumbling rate so much in comparison to rods. It is not easy to directly compare the degree of alignment of rods and disks because disks have an entire plane that can align with the stretching, while rods have only a director. The picture that emerges from our data and previous work is that anisotropic particles are aligning with the Lagrangian stretching direction, and this suppresses the tumbling rate for all particles except nearly spherical oblate ellipsoids. The amount of suppression is strongly dependent on the axis whose tumbling rate is considered.

Conclusions
We have developed a method for measuring the time-resolved Lagrangian orientation and solid body rotation rate of anisotropic particles in a turbulent flow. By measuring the rotation of 3D printed jacks and crosses we are able to extend previous measurements of rods to cover the full range of aspect ratios of axisymmetric ellipsoids. Moreover, we have provided a way to directly probe Lagrangian vorticity with a single particle measurement, which has potential for application in a wide range of flows and Reynolds numbers. We find that the mean square tumbling rate, ṗ iṗi , agrees with DNS data at points spanning the full range of aspect ratios. Our measurements show that crosses are preferentially aligned in turbulence with their orientation axis perpendicular to their solid body rotation rate vector. To the best of our knowledge, this is the first direct experimental measurement of the preferential alignment of particles with their rotational motion in a turbulent flow. Our results support a natural picture of alignment in turbulence where particles have their long axes aligned with the Lagrangian stretching direction of the fluid flow.

Appendix A. Resistive Force Theory
Resistive force theory is widely used to model the motion of slender bodies at small Reynolds number (see, e.g., [40]). We will apply it first to a small thin rod of width a in a general velocity field. Then we will generalize to multiple perpendicular rods in order to show that jacks rotate like spheres and crosses rotate like disks.
Resistive force theory assumes that the differential viscous drag forces acting on a small segment of a slender body are [41]: where df is parallel to the rod, while dF is perpendicular to the rod and responsible for producing all of the torque. The vector u is the relative velocity between the segment and the fluid, and its subscripts indicate the component parallel or perpendicular to the slender object. We can use these expressions to calculate the total torque due to fluid flow past the rod. Since the particles are on the order of the Kolmogorov scale, we linearize the velocity field: Taking the center of the particle to be at rest at the origin, the velocity can be written as: where p is the normalized orientation vector of the rod and r is the radial distance from the center of the rod. The drag force is determined by the relative velocity between the fluid and the rigid-body rotation of the rod: The toque comes from the component of velocity perpendicular to the rod, u ⊥ . We can write this by subtracting the parallel component from the full relative velocity. Entering u ⊥ into Equation A.1, we have: where we have written the velocity gradient in terms of its symmetric and antisymmetric parts, S ij and Ω ij , respectively. We now rearrange (A.4), noting that p k Ω kl p l = 0 since p k p l forms a symmetric tensor while Ω kl is antisymmetric, and compute the differential torque element: Integrating along the length of the rod, l, we have the total torque, Now, consider the term imn (ω s ) m p n = (ω s × p) i . This is nothing else than the rate at which the orientation of the rod is changing,ṗ. Thus, we have: For any orientation, p, (A.7) can be used to compute the qth component of the torque on a single rod. For Stokes flow: τ q = 0, for q=1,2, and 3. (A.8) When we impose this constraint on (A.7), we immediately recover Jeffery's equation (Equation 1) in the limit of infinite aspect ratio (i.e., for a rod): p i = Ω ij p j + S ij p j − p i p k S kl p l (A.9) With this established, we can now explore the rotational dynamics of objects that are composites of perpendicular rods with arbitrary lengths. From (A.6), we can read off the qth component of the torque on the αth arm, τ α q . The total qth component of the torque is then: where N ≤ 3 is the number of arms on the particle, l α is the length of the αth arm, and Γ α is the pre-factor for the αth arm, Γ α ≡ aC d l 3 α . From here, we can extend the calculation made for single rods to particles composed of N rods. The zero net torque condition of Stokes flow still holds true, but there will now be contributions from multiple arms because of the sum over α. If we consider axisymmetric particles, l 1 = l 2 = l 3 , and are careful to define p along the symmetry axis of the particle, Equation A.10 leads to an analog of Jeffery's equation: where α p ≡ l 3 /l 1 . We find that resistive force theory for the rotational motion of a set of three perpendicular rods with l 1 = l 2 = l 3 reproduces Jeffery's equation, except that the aspect ratio appears with cubic dependence rather than quadratic, such that the particles behave like ellipsoids with an effective aspect ratio [23]. The cases we study are jacks, l 1 = l 2 = l 3 with α p = 1, and crosses, l 1 = l 2 but l 3 = 0 with α p = 0 . These both happen to be cases where α 3 p = α 2 p , so, within the limits of resistive force theory, we can simply say that jacks rotate like spheres and crosses rotate like disks.