1. Introduction
Particles carried by turbulent flows in nature, such as ice crystals in clouds (Korolev & Isaac Reference Korolev and Isaac2003; Shaw Reference Shaw2003; Pinsky, Khain & Shapiro Reference Pinsky, Khain and Shapiro2007) or plankton in the oceans (Karp-Boss, Azevedo & Boss Reference Karp-Boss, Azevedo and Boss2007; Jumars et al. Reference Jumars, Trowbridge, Boss and Lee2009; Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012), are rarely spherical; instead, they often have non-trivial, anisotropic shapes that may influence their motion as they are carried by the flow. The simplest shape to consider after a sphere is that of an ellipsoid. The tumbling rate of an axi-symmetric ellipsoid in Stokes flow is determined by the particle orientation and the velocity gradient tensor (Jeffery Reference Jeffery1922):
$$\begin{eqnarray}{\dot{p}}_{i}={\rm\Omega}_{ij}p_{j}+\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}(\unicode[STIX]{x1D61A}_{ij}p_{j}-p_{i}p_{k}\unicode[STIX]{x1D61A}_{kl}p_{l}),\end{eqnarray}$$
where
$\hat{\boldsymbol{p}}$
is a unit director along the symmetric axis of the particle,
${\it\alpha}$
is the aspect ratio of the ellipsoid, and
$\unicode[STIX]{x1D61A}_{ij}$
and
${\rm\Omega}_{ij}$
are the symmetric and antisymmetric parts of the velocity gradient tensor, respectively.
In a turbulent flow, a small ellipsoidal particle rotates in response to the velocity gradients along its Lagrangian trajectory. Because these Lagrangian velocity gradients are controlled by the small scales, they are similar in many different turbulent flows and have been the focus of extensive study (Meneveau Reference Meneveau2011). To understand the dynamics of ellipsoidal particles in turbulence, we need to extend our understanding of the Lagrangian statistics of the velocity gradient tensor to include the orientational dynamics that result from integrating equation (1.1) along particle trajectories. This is a challenging problem, both because of the complexity of statistically quantifying the particle orientation with respect to the velocity gradient tensor, and due to the difficulty of measuring the dynamics of anisotropic particles simultaneously with the velocity gradient tensors.
The study of anisotropic particles in fluid flows has a long history because of the many relevant applications. Leal (Reference Leal1980) provides a review of the older literature, and a wide range of work has followed, for example: Koch & Shaqfeh (Reference Koch and Shaqfeh1989), Szeri & Leal (Reference Szeri and Leal1993), Herzhaft et al. (Reference Herzhaft, Guazzelli, Mackaplow and Shaqfeh1996), Olson & Kerekes (Reference Olson and Kerekes1998), Parsa et al. (Reference Parsa, Guasto, Kishore, Ouellette, Gollub and Voth2011), Andersson & Soldati (Reference Andersson and Soldati2013) and Rosen, Lundell & Aidun (Reference Rosen, Lundell and Aidun2014). Turbulent flows advecting anisotropic particles provide a compelling test case, both because of the many applications and because of the nearly universal statistics of the velocity gradients experienced by small particles in many turbulent flows at large Reynolds number. However, the difficulty of accessing particle and fluid variables in turbulent flows has hindered work in this area. Rod-shaped particles (
${\it\alpha}\gg 1$
) were the first to be studied. Shin & Koch (Reference Shin and Koch2005) provided an extensive numerical study of rotational diffusion and the tumbling rate of rods in turbulence. They observed that the tumbling rate of rods is much smaller than that predicted for randomly oriented rods. Pumir & Wilkinson (Reference Pumir and Wilkinson2011) showed from numerical simulations that this suppression of the tumbling rate is caused by rods aligning with the vorticity vector. Parsa et al. (Reference Parsa, Calzavarini, Toschi and Voth2012) extended numerical study of the tumbling rate across the full range of aspect ratios of axi-symmetric ellipsoids and found that preferential alignment decreases the tumbling rate for almost all shapes. They also provided the first time-resolved experimental measurements of tumbling of rods in turbulence. Chevillard & Meneveau (Reference Chevillard and Meneveau2013) studied the full parameter space of tri-axial ellipsoids in numerical simulations and showed that the tumbling of rods is a challenging test case for stochastic models of the velocity gradient tensor in turbulence. Gustavsson, Einarsson & Mehlig (Reference Gustavsson, Einarsson and Mehlig2014) used analytical and numerical methods to show that the differences in tumbling between rods and disks can be understood using Lagrangian three-point correlations of the velocity gradient tensor. Parsa & Voth (Reference Parsa and Voth2014) made experimental measurements of the rotation of rods with lengths extending into the inertial range of turbulence and proposed that rotations of long rods should show inertial range scaling. Ni, Ouellette & Voth (Reference Ni, Ouellette and Voth2014) showed how Lagrangian stretching aligns the long axis of a particle with the vorticity in turbulence.
A question that remains unanswered is how to quantify the preferential orientation of particles that decreases their tumbling rate in turbulence. From numerical simulations (Pumir & Wilkinson Reference Pumir and Wilkinson2011; Gustavsson et al.
Reference Gustavsson, Einarsson and Mehlig2014) we know the probability distributions of the projection of
$\hat{\boldsymbol{p}}$
onto the vorticity and the eigenvectors of the strain rate. But predicting the tumbling rate requires the full joint probability distribution of the velocity gradient tensor and the particle orientation. This could be obtained from numerical simulations, although a thorough study has not yet been published. In this paper we show that this joint probability distribution is now accessible to experimental measurements.
One of the main challenges in simultaneous measurement of the orientations of particles and the velocity gradient tensor along the particle trajectory is that measuring the Lagrangian velocity gradient tensor is difficult and typically has large experimental uncertainties. This is primarily because the Kolmogorov length scale
${\it\eta}$
over which the flow is roughly linear is small, typically tens or hundreds of microns at high Reynolds numbers. Achieving this spatial resolution in an experiment is highly non-trivial. Multi-sensor hot-wire probes (Wallace & Vukoslavc̆ević Reference Wallace and Vukoslavc̆ević2010) can provide the required resolution, but acquire only single-point Eulerian information and typically require the use of Taylor’s hypothesis and a strong mean velocity. Non-intrusive optical methods, such as laser-induced fluorescence (LIF), particle image velocimetry (PIV), and particle tracking, are potentially viable alternatives. Particle tracking is best suited for our purposes, as it directly provides Lagrangian information without requiring interpolation or integration of velocity fields. But since particle tracking follows the motion of individual tracer particles, it has typically been restricted to fields that are sampled too sparsely to resolve the velocity gradient. At small Reynolds numbers, however, Lüthi, Tsinober & Kinzelbach (Reference Lüthi, Tsinober and Kinzelbach2005) successfully used particle tracking to measure the velocity gradient; thus, if the seeding density can be made large enough, the gradient can be resolved. One promising method for increasing the seeding density in particle tracking is to illuminate not the entire measurement volume at once but rather to section it into successively illuminated slabs by scanning a laser through it (Hoyer et al.
Reference Hoyer, Holzner, Lüthi, Guala, Liberzon and Kinzelbach2005). This scanning particle tracking opens the door for using particle tracking to measure the velocity gradient at larger Reynolds numbers, and has been used, for example, along with LIF to study the joint evolution of the velocity gradient and the density field in turbulent gravity currents (Krug et al.
Reference Krug, Holzner, Lüthi, Wolf, Tsinober and Kinzelbach2014).
In this paper, we report experimental results on the joint dynamics of rods and the velocity gradient tensor in turbulence using a scanning particle tracking system. In § 2, we show an example of a rod trajectory along with the local velocity gradient tensor to give a qualitative idea of the different ways the gradient contributes to the rod tumbling rate. In § 3, we present in detail our experimental setup and our scanning particle tracking system; an analysis of our experimental uncertainties in measuring the gradient is reserved for the appendices. Details of the data analysis are addressed in § 4. In § 5, we discuss our experimental results concerning the dependence of the rod tumbling rate on the velocity gradient; in particular, we describe quantitatively the relative contributions from the vorticity and the strain rate and show that both are necessary for understanding the orientational dynamics of rods.
2. Example
Before describing our techniques for measuring rod motion and the velocity gradient tensor in detail, we show in figure 1 an example of a measured rod trajectory along with vectors characterizing the local velocity gradients. The ribbon shows the full trajectory with a solid rod plotted only once every 25 time steps. The colour indicates the magnitude of rod tumbling rate
$|\dot{\boldsymbol{p}}|$
, which depends on the straining and swirling motion of the surrounding flow. The tumbling rate due to strain is
$$\begin{eqnarray}{\dot{p}}_{i}^{S}=\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}(\unicode[STIX]{x1D61A}_{ij}p_{j}-p_{i}p_{k}\unicode[STIX]{x1D61A}_{kl}p_{l})\end{eqnarray}$$
and this tends to align the rod with the most extensional direction of the local flow (red arrow), which is given by the eigenvector
$\hat{\boldsymbol{e}}_{1}$
of the strain-rate tensor
$\unicode[STIX]{x1D61A}_{ij}$
that corresponds to its largest eigenvalue
${\it\lambda}_{1}$
. Local swirling is characterized by the rotation-rate tensor
${\rm\Omega}_{ij}$
, which tends to rotate the rods about the local vorticity direction
$\hat{{\bf\omega}}$
(blue arrow) at a rate of
$$\begin{eqnarray}{\dot{p}}_{i}^{{\rm\Omega}}={\rm\Omega}_{ij}p_{j}.\end{eqnarray}$$
Reconstructed three-dimensional trajectory of a rod (
$700~{\rm\mu}\text{m}$
in length) over 1 s (527 frames). The colourmap shows the tumbling rate of a rod. The blue and the red arrows represent the vorticity vector
${\bf\omega}$
and the largest stretching direction
$\hat{\boldsymbol{e}}_{1}$
, respectively. The length of the two arrows indicates the magnitude of
${\bf\omega}$
and the eigenvalue
${\it\lambda}_{1}$
of strain-rate tensor corresponding to the largest stretching direction
$\hat{\boldsymbol{e}}_{1}$
.
For this rod trajectory, the magnitude of the total tumbling rate
$\dot{\boldsymbol{p}}$
as well as its two component contributions
$\dot{\boldsymbol{p}}^{{\rm\Omega}}$
and
$\dot{\boldsymbol{p}}^{S}$
are shown in figure 2(a,b). The total tumbling rate can be computed in two ways: by differentiating the rod orientation signal or by using measurements of the velocity gradient tensor and Jeffery’s equation
$\dot{\boldsymbol{p}}^{J}=\dot{\boldsymbol{p}}^{{\rm\Omega}}+\dot{\boldsymbol{p}}^{S}$
. As shown in figure 2(a), the two measurements agree well with each other, indicating both that our measurement of the velocity gradient tensor is accurate and that the rod is small enough that Jeffery’s equation holds. In figure 2(c), we plot the cosine of the angle between
$\dot{\boldsymbol{p}}^{{\rm\Omega}}$
and
$\dot{\boldsymbol{p}}^{S}$
. When this quantity is negative, the contribution to the rod tumbling due to strain works against that due to rotation. But when it is positive, the two contributions work cooperatively and lead to large tumbling rates, as can be seen at 0.16 and 0.56 s in figure 2 (corresponding to the red regions in figure 1).
The time series of the rod tumbling rate for the trajectory shown in figure 1. (a) The magnitude of total tumbling rate from two different measurements:
$\dot{\boldsymbol{p}}$
determined by differentiating the rod orientation (curve) and
$\dot{\boldsymbol{p}}^{J}$
calculated from Jeffery’s equation using the velocity gradient measurements around the rod (circles). (b) Two contributions to tumbling rate from vorticity
$|\dot{\boldsymbol{p}}^{{\rm\Omega}}|$
and from strain
$|\dot{\boldsymbol{p}}^{S}|$
. (c) The cosine of the angle between two vectors
$\dot{\boldsymbol{p}}^{{\rm\Omega}}$
and
$\dot{\boldsymbol{p}}^{S}$
.
3. Experimental apparatus and setup
We generated a turbulent flow in an octagonal Plexiglas water tank measuring
$1\times 1\times 1.5~\text{m}^{3}$
. Two grids with a mesh size of 8 cm were oscillated in phase with an amplitude of 12 cm. Details about the apparatus are given in Blum et al. (Reference Blum, Kunwar, Johnson and Voth2010). In this paper, all experiments were performed with a grid frequency of 1 Hz which creates turbulence in the centre of the tank with Reynolds number
$R_{{\it\lambda}}=140$
. The Kolmogorov length scale is
${\it\eta}=({\it\nu}^{3}/{\it\epsilon})^{1/4}=0.31$
mm and the Kolmogorov time scale is
${\it\tau}_{{\it\eta}}=({\it\nu}/{\it\epsilon})^{1/2}=0.093$
s, where
${\it\nu}$
is the kinematic viscosity and
${\it\epsilon}$
is the energy dissipation rate per unit mass whose measurement is described in appendix B. The temperature of the water was almost constant at
$21.8\pm 0.2\,^{\circ }\text{C}$
, giving a kinematic viscosity of
${\it\nu}=(9.61\pm 0.05)\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$
.
To measure the velocity gradient tensor simultaneously with rod motion, we need two kinds of particles: the rods themselves and small, spherical tracer particles. A key to this experiment is obtaining suitable particles. Fluorescent particles are convenient, as they will produce much better images, since some small residue from aluminium top and bottom plates and the bearings may be carried into the measurement volume. For tracer particles, we therefore used internally dyed polystyrene divinylbenzene (PS-DVB) particles with
$30~{\rm\mu}\text{m}$
diameter and density
$1.05~\text{g}~\text{cm}^{-3}$
, purchased from Thermo Scientific. The rods were nylon fibre (density
$1.10{-}1.13~\text{g}~\text{cm}^{-3}$
) from DonJer Corp., with major and minor axes of roughly 700 and
$30~{\rm\mu}\text{m}$
, respectively, giving an aspect ratio of
$23.3$
. We dyed the rods with Rhodamine-B, so that they absorbed green light (wavelength
${\it\lambda}=532$
nm) from our laser and emitted red light at the same wavelength as the tracer particles. A typical image of these fluorescent particles, captured through a Schneider B
$+$
W MRC Orange 550 band-pass filter, is shown in figure 3. The shapes of the two types of particles are well-defined and distinct from each other, and thus the particles can be easily separated using automated image analysis. Note that the brightness of the two kinds of particles is very similar, another important factor in the experiment: if the brightness of two types of particles were very different, it would be difficult to determine their positions accurately at the same time. Even though the tracer particles are much smaller than the rods, the internal dyeing emitted a much stronger fluorescent signal than the rods, which we dyed ourselves.
In our subsequent analysis, we make the assumption that the rods and the tracers do not exhibit inertial effects caused by their finite size or density difference with respect to the fluid. To characterize the validity of this assumption, we consider the lengths and Stokes numbers of the particles. The length of a rod is roughly 2.3
${\it\eta}$
. Studies of long neutrally buoyant rods show that the tumbling of a rod of this length is very close to the short rod limit (Shin & Koch Reference Shin and Koch2005; Parsa & Voth Reference Parsa and Voth2014). The tracer particles have diameters of 0.1
${\it\eta}$
, which is clearly in the small sphere limit (Voth et al.
Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002). The Stokes number
$St={\it\tau}_{p}/{\it\tau}_{{\it\eta}}$
is defined as the ratio between the time scale of the Stokes viscous drag
${\it\tau}_{p}=r^{2}/(3{\it\beta}{\it\nu})$
and the Kolmogorov time scale
${\it\tau}_{{\it\eta}}$
. Here,
$r$
is the radius of the particle and
${\it\beta}=3{\it\rho}_{f}/(2{\it\rho}_{p}+{\it\rho}_{f})$
is a coefficient capturing the effect of a density difference between the fluid (with density
${\it\rho}_{f}$
) and the particle (with density
${\it\rho}_{p}$
). If the particle response time
${\it\tau}_{p}$
is much smaller than the smallest scale
${\it\tau}_{{\it\eta}}$
, the inertial effect of the particle is negligibly small, and the particles can be safely treated as tracers. For our spherical particles,
$St=8.6\times 10^{-5}\ll 1$
. For a rod, the relaxation time is given by
${\it\tau}_{r}=2{\it\alpha}{\it\rho}_{p}r^{2}\ln ({\it\alpha}+\sqrt{{\it\alpha}^{2}-1})/9{\it\nu}\sqrt{{\it\alpha}^{2}-1}{\it\rho}_{f}$
, where
$r$
is the semi-minor axis (Shapiro & Goldenberg Reference Shapiro and Goldenberg1993; Zhang et al.
Reference Zhang, Ahmadi, Fan and McLaughlin2001). Using this expression, the Stokes number for the rods is
${\it\tau}_{r}/{\it\tau}_{p}=8\times 10^{-3}$
, which is also much less than unity. Thus, both the rods and the spherical particles are in the tracer limit, and will follow the fluid motion accurately.
The full resolution (a) and zoomed-in (b) images captured by one of the three cameras. The circle in (b) is centred at the centroid of a rod with 2 mm (6.7
${\it\eta}$
) radius. The tracer particles that fall in such a sphere will be used to calculate the velocity gradient tensor.
Schematic of the experimental apparatus. Lens group 1 controlled the thickness of the laser beam and needed to be placed before the piezo-electric mirror. Lens group 2 (
$\,f=15$
mm and
$f=150$
mm) expanded the beam in the vertical direction to form the final laser sheet. The laser sheet stayed nearly parallel across the whole scanning range (divergence is exaggerated in the figure). Detection volume is not to scale. Camera 3 is located below camera 1.
3.1. Scanning system
Once the particles are chosen, we need to determine the proper seeding density for the experiment. Generally speaking, we want a very high density of spherical particles and relatively low density of rods. The density of spherical particles is directly related to the spatial resolution, which is crucial for measuring the velocity gradient tensor. The maximum density is limited, however, if we illuminate the entire measurement volume, because the images of individual particles will overlap with each other when the particle seeding density is high. We therefore use a scanning particle tracking system for our measurements (Hoyer et al. Reference Hoyer, Holzner, Lüthi, Guala, Liberzon and Kinzelbach2005). The basic principle of this technique is a sacrifice of temporal resolution for improved spatial resolution. By subdividing the measurement volume into 10 slabs, for example, and successively scanning the illumination through them, we can increase the particle seeding density by a factor of 10 (in the ideal case), albeit at the cost of requiring a factor of 10 increase in camera frame rate and a decrease in total recording time.
A schematic of the experimental setup is shown in figure 4. The beam from a pulsed Nd:YAG laser with an average output power of 50 W was stretched independently in height and in width by two sets of lenses to create an illumination slab measuring 50 mm
$\times$
3 mm. To ensure a relatively uniform slab depth over a height of 30 mm, we did not image the top and bottom 10 mm of the laser slab. The scanning of the light slab was controlled by a piezo-electric driven mirror (S-224, PI Inc.) with a diameter of 15 mm and a maximum deflection angle of 2.2 mrad with sub-
${\rm\mu}\text{rad}$
resolution. The small deflection angle is magnified through the optical system to give a 12 mm scanning range in the centre of the flow chamber. For reasons we will explain in § 4.2, each slab overlaps with the previous one by
${\sim}50$
%. Compared to previous designs using a rotating prism to scan the illumination slab (Hoyer et al.
Reference Hoyer, Holzner, Lüthi, Guala, Liberzon and Kinzelbach2005), a piezo-electric driven mirror or an acousto-optic deflector has the potential to generate faster scanning rates, which are more suitable for turbulence with even higher Reynolds number.
Three Photron FASTCAM SA5 cameras with a resolution of
$1024\times 1024$
pixels were used to image the particles in a small volume measuring approximately
$3\times 3\times 4~\text{cm}^{3}$
in the centre of the tank. To resolve this small volume more than a half metre away from the sidewalls of the tank, each camera was fitted with a Nikkor 200 mm macro lens and a Kenko 1.6 teleconverter. The cameras were mounted on a custom-built frame on an optical table uncoupled from the turbulence tank to minimize camera vibration. Two of the cameras (labelled as ‘top cameras’ in figure 4) were mounted in the same lateral plane with a
$90^{\circ }$
angular separation, and looked down into the detection volume at an angle of
$16^{\circ }$
. The third camera was aimed up toward the same volume at an angle of
$21^{\circ }$
. Situating the three cameras in different planes helps to increase the stereomatching accuracy (Ouellette, Xu & Bodenschatz Reference Ouellette, Xu and Bodenschatz2006).
The piezo-electric mirror was driven with an adjusted saw-tooth signal at a frequency of 500 Hz, rising linearly for 80 % of each cycle and falling for the remaining 20 %. This driving produced a nearly linear scanning motion of the laser beam through the measurement volume followed by a quick return to the initial position. The cameras recorded images at a frame rate of 5000 Hz, so that each cycle of the mirror resulted in 10 captured sets of images, with eight of these linearly positioned through the measurement volume. We set the cameras to store only eight out of every 10 frames to conserve on-board memory, which could hold a total of 5456 images. Thus, approximately 1.4 s (14.7
${\it\tau}_{{\it\eta}}$
or 0.44 large-eddy turnover time) of data is held in the camera memory before being transferred to the computer hard drive. Our experimental protocol thus consisted of several steps. First, the two grids were driven at 3 Hz for 30 s to stir up the fluid and particles. The grids were then slowed to a 1 Hz oscillation rate, and the flow was given 1 min (approximately 17.5 large-eddy turnover times) to stabilize. The cameras then recorded images until their memory was full (for approximately 1 s), and the system subsequently rested until the data transfer to the hard drive was complete. The timing of this system was automatically controlled with Labview scripts. The results reported here come from a full day of measurements, giving roughly 300 data sets.
3.2. Simulation
To compare with the experimental data, one dataset from a direct numerical simulation (DNS) of homogeneous isotropic turbulence at Reynolds number
$R_{{\it\lambda}}=180$
was used. There are total of
$N^{3}=512^{3}$
collocation points for the entire volume. A total of
$7\times 10^{4}$
Lagrangian trajectories of the velocity gradient tensor were followed for
$O(1)$
large-eddy turnover times, and the time step for integrating the Navier–Stokes equations and tracking Lagrangian points was
$O(10^{-2}{\it\tau}_{{\it\eta}})$
. Along each Lagrangian trajectory, the orientation of a virtual infinitesimal rod with an aspect ratio of
$20$
was computed by integrating Jeffery’s equation (1.1) using a fourth-order Runge–Kutta method (Parsa et al.
Reference Parsa, Calzavarini, Toschi and Voth2012). The details of this simulation are given by Benzi et al. (Reference Benzi, Biferale, Calzavarini, Lohse and Toschi2009).
4. Data analysis
4.1. Image processing
To determine particle positions, the digital image from each camera is first segmented into groups of bright pixels, representing both tracers particles and rods. Typically, the image of one spherical tracer particle contains
$4{-}9$
pixels, corresponding to 2–3 pixels in diameter. The number of pixels in a rod image depends on its relative orientation to the camera. The maximum size is
${\sim}320$
pixels when its long axis is perpendicular to the optical path of the camera, in which case it is also highly anisotropic. But the image of a rod becomes quasi-spherical with an area of
${\sim}16$
when it points directly toward the camera. Ideally, the images of rods and spherical particles could be separated solely by their area. However, sometimes a spherical particle looks larger if it is out of focus. So we add another criterion for separating the images of spheres and rods: eccentricity, which is defined as the degree of geometrical deviation from a circle. It is determined by finding the best-fit ellipse to the pixel cluster and calculating the ratio of the distance between the foci and the major axis of the ellipse. We find that the eccentricity ranges from 0 to 0.6 for spherical particles and from 0.8 to 1 for rods. In practice, a particle image is considered to be a rod if its eccentricity is larger than 0.9 and its area is larger than 30 pixels. These criteria separate almost all rods from spherical tracer particles. Sometimes a rod that points directly toward one camera will be mistakenly identified as a tracer in that camera; but because we have other cameras viewing the same rod from different directions, where it will be distinctly elongated, it can be correctly identified by combining the information from all the cameras.
Diagram of one particle passing from one slab
$n$
(box in solid line) to the next one
$n+1$
(box in dashed line). The positions of the same particle in slab
$n$
and
$n+1$
are represented by solid circles and crosses, respectively. The time between two consecutive volume scans is
${\rm\Delta}t$
, and the time between two slabs is
${\it\delta}t$
. In the overlapping region, there are three pairs of joins, which can be used to connect two segments of trajectories in different slabs together.
(a) The trajectories that pass through at least five laser slabs in one typical movie with 5457 frames. Particle positions that belong to different laser slabs are shown by a different colour. (b) The histogram of
$y$
positions of the particles in each slab. Only particles from the centre of the detection volume (
$-5~\text{mm}<x<5~\text{mm}$
and
$-5~\text{mm}<z<5~\text{mm}$
) were used to ensure that the histogram represents the slab width and overlap.
4.2. Particle tracking
Once we know the particle type, its centre is determined by averaging the positions of its pixels weighted by their brightness values. This procedure is used for both types of particles. From the two-dimensional positions of the particles as measured by each camera, the three-dimensional positions can be found by stereoscopic matching, and subsequently connected together in time from one volume scan to the next (Ouellette et al. Reference Ouellette, Xu and Bodenschatz2006), although some modifications of common tracking algorithms are required.
In traditional particle tracking, the entire measurement volume is illuminated and is imaged at a constant rate of
$1/{\rm\Delta}t$
, where
${\rm\Delta}t$
is the time interval between two frames. Thus, finding candidate particles to extend a given trajectory requires searching for potential matches only at a time
${\rm\Delta}t$
in the future. But for the scanning system, as is sketched in figure 5, there are two relevant time intervals:
${\rm\Delta}t=1/500$
s, the time between two full scans of the volume; and
${\it\delta}t=1/5000$
s, the time between the illumination of two neighbouring slabs of the volume. Since a particle may or may not pass from one slab to another over an interval of
${\it\delta}t$
, the time at which a particle corresponding to the continuation of a trajectory may be found is not obvious: it may be found at
${\it\delta}t$
in the future, for example, if the particle moved between two successive slabs; at
${\rm\Delta}t$
in the future if it remained in the same slab; at
${\rm\Delta}t-{\it\delta}t$
in the future if it moved to a previous slab; or even potentially at other times. A previous scanning tracking system by Hoyer et al. (Reference Hoyer, Holzner, Lüthi, Guala, Liberzon and Kinzelbach2005) performed the search for the continuation of a particle trajectory in slabs
$n$
,
$n-1$
, and
$n+1$
. To handle this ambiguity, we first track the particles found in each individual slab separately, and then subsequently merge the short segments that belong to different slabs but that refer to the same particle. Within each slab, particles are tracked using a standard predictive tracking algorithm (Ouellette et al.
Reference Ouellette, Xu and Bodenschatz2006), with a recording time of roughly
${\sim}d/{\rm\Delta}t\tilde{u} =75$
frames, where
$d=3$
mm and
$\tilde{u} \sim 2~\text{cm}~\text{s}^{-1}$
are the thickness of one laser slab and the root-mean-square velocity of the flow, respectively. In figure 6(b), it is seen that each pair of neighbouring slabs overlaps by
$d_{ov}\sim 1.5$
mm. Thus, since the velocity of a particle is bounded, a particle that moves from one slab to the next will inevitably pass through the overlapping region and will be recorded twice during each volume scan. We refer to each doubly recorded position as a ‘join’ between the two trajectory segments. In figure 5, where we schematically demonstrate our tracking method, there are four joins. In principle, two joins are enough to merge two trajectory segments into one longer track. Estimating the number of joins in the overlapping region using simple kinematics gives
$d_{ov}/{\rm\Delta}t\tilde{u} =32$
. Thus, even a particle with a speed of
$10\tilde{u}$
would still have more than two joins in the overlapping region. We apply this merging procedure consecutively for all neighbouring slabs, and thus link all segments in different slabs together into longer trajectories. In figure 6(a), we show those trajectories that pass through at least five different slabs, which evenly cover the entire volume, to demonstrate qualitatively that our measurements are robust and that the overlapping regions are large enough for splicing.
The stereomatching and tracking procedures for rods are almost the same as for the tracers. The primary difference is that we need to keep track of the orientation of a rod in addition to its position. In two dimensions, the orientation, defined as the angle between the major axis of the rod and horizontal axis of the image, can be extracted from all three cameras. The three-dimensional orientation of rods can then be uniquely determined from these three two-dimensional angles and the viewing directions of the cameras (Parsa Reference Parsa2013). Errors in determining the orientation arise mainly from the finite aspect ratio of the rods. In the experiments reported here, the aspect ratio of the rods is 20, roughly four times that in our previous experiments (Parsa et al. Reference Parsa, Calzavarini, Toschi and Voth2012). Thus, the uncertainty in orientation in these data is smaller than it was for our previous results (Parsa Reference Parsa2013).
4.3. Interpolation and differentiation
The next step in processing the data is to filter and differentiate the trajectories to obtain the velocities of the tracers and tumbling rates of the rods. Common methods to accomplish this task include convolving the trajectory with appropriate kernels (Mordant, Crawford & Bodenschatz Reference Mordant, Crawford and Bodenschatz2004) and polynomial fitting (Voth et al.
Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002). In the scanning system, the convolution method is difficult to apply because the points along the trajectories are not evenly spaced in time, leading to difficulties in accurate numerical integration. We therefore use polynomial fitting. In addition to smoothing and differentiating, fitting also allows us to accurately interpolate the measured positions and velocities in time so that all data throughout the measurement volume are contemporaneous. That is, we can use the fits to extract the kinematics of the flow field not at the measured space–time positions of the particles but rather at the times
$t_{n}=4{\it\delta}t+n{\rm\Delta}t$
(
$n=0,1,2,3,\dots$
), so that we acquire one full velocity field for each scan of the volume, measured at the time corresponding to the illumination of slab 5 (halfway through the volume scan). To accomplish this interpolation, for each
$n$
, we fit a polynomial to all measured data points in the range
$[t_{n}-({\it\tau}_{f}-1)/2,t_{n}+({\it\tau}_{f}-1)/2]$
along each track, where
${\it\tau}_{f}=9{\rm\Delta}t$
is the temporal length of the fit. This choice minimizes noise without unduly affecting the signal (Voth et al.
Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002). From the polynomial fits, we extract smoothed positions, velocities, and accelerations of the tracers at
$t_{n}$
spread over the entire volume. We apply a similar process to the rods to extract their orientation and tumbling rate. Note, however, that we define the orientation of the rods by a unit vector along their major axis. Smoothing the orientation signal decreases random error in the measurement of the orientation direction, but may change the magnitude of this vector. Therefore, we must re-normalize the orientation vector for each rod after smoothing.
4.4. Seeding density and shape factor
To determine how large we could make the seeding density of the spherical tracer particles and still obtain good measurements, we tested the scanning system by slowly increasing the number of tracer particles. When the particle concentration is low, the ratio between the number of successfully stereomatched particles and the number of particles detected in each two-dimensional image is almost constant. We can determine the maximum seeding density at which we can still resolve the particles by locating the point at which this ratio begins to decrease. In our experiments, this point corresponds to roughly 500 stereomatched particles in each slab. After accounting for doubly imaged particles in the overlapping regions and trajectories shorter than
${\it\tau}_{f}$
, for which we cannot measure velocities or accelerations, we can reliably measure about 2000 velocity vectors in each volume scan. This number varies somewhat over time due to sedimentation of particles to the sidewalls and bottom plate, so we add particles over the course of an experimental run to maintain a roughly constant particle concentration. The seeding density of rods is kept low to avoid interactions between them. We have roughly 10 rods in the measurement volume at any given time, so the non-dimensional concentration is roughly
$n\ell ^{3}\sim 10^{-3}$
, where
$n$
is the number density and
$\ell$
is the length of a rod. This is far below the concentration at which rod–rod interactions become important.
After tracking, smoothing, and differentiation, we are left with trajectories of rod orientations along with the velocities of many tracers surrounding them. Measuring the velocity gradient tensor around each rod requires us to estimate the spatial gradient from multiple velocity vectors inside a small volume. To do this, we first locate all tracers within a 2 mm (
${\sim}6{\it\eta}$
) radius of the centre of a rod. This radius is chosen to have sufficient tracer particles (
$6{-}10$
, given our tracer-particle seeding density) surrounding the rod to estimate the velocity gradient well. This number is comparable to what has been used in previous experiments, where a radius of nearly
$8{\it\eta}$
was found to be sufficient so that viscous effects dominate and the velocity field is close to linear (Lüthi et al.
Reference Lüthi, Tsinober and Kinzelbach2005).
To estimate the gradient, consider
$N$
tracers at positions
$\boldsymbol{x}^{n}(t)$
and with velocities
$\boldsymbol{u}^{n}(t)$
(
$n=1,2,3,\dots ,N$
) that are sufficiently close to a rod. Their relative position and velocity with respect to the centre of mass of the particle cloud are
$\boldsymbol{x}^{\prime n}(t)=\boldsymbol{x}^{n}(t)-\sum _{n}\boldsymbol{x}^{n}(t)/N$
and
$\boldsymbol{u}^{\prime n}(t)=\boldsymbol{u}^{n}(t)-\sum _{n}\boldsymbol{u}^{n}(t)/N$
, respectively. Determining the velocity gradient tensor
$\unicode[STIX]{x1D608}_{ij}$
can then be formulated as a least-squares problem by finding the minimum value of the squared residuals
$S=\sum _{n}[u_{i}^{\prime n}(t)-\unicode[STIX]{x1D608}_{ij}x_{j}^{\prime n}(t)]^{2}$
(Pumir, Bodenschatz & Xu Reference Pumir, Bodenschatz and Xu2013). In general, however, the tracer particles in the cloud surrounding a rod are randomly distributed in space, and will sometimes lie in almost the same plane. Such cases will introduce a large error in the determination of the out-of-plane components of the velocity gradient. To exclude these cases, we use the inertia tensor
$\unicode[STIX]{x1D644}=\unicode[STIX]{x1D65C}/\text{tr}(\unicode[STIX]{x1D65C})$
with
$\unicode[STIX]{x1D628}_{ij}=\sum _{n}x_{i}^{n}x_{j}^{n}$
to characterize the shape of the tracer-particle cloud.
$\unicode[STIX]{x1D644}$
can be diagonalized in an orthogonal basis with three eigenvalues
$I_{1}\leqslant I_{2}\leqslant I_{3}$
. For a symmetric object,
$I_{1}=I_{2}=I_{3}=1/3$
, while for a co-planar particle cloud
$I_{3}\approx 0$
. In our experiment, we require
$I_{3}>0.1$
to rule out cases that will have large errors. Empirically, we find that our results are not sensitive to the choice of this threshold when it is in the range from 0.07 to 0.2. In addition to the overall shape of the particle cloud, the distribution of the tracer particles inside the cloud may also affect the estimate of the velocity gradient: sometimes, the particles will all be concentrated in a corner, for example, and the centre of mass of the tracer cloud will be far away from the rod centre. This situation will lead to a biased estimate of the velocity gradient at the rod centre. To avoid this case, we require that the distance between the rod centre and centre of the tracer cloud be no more than one third of the radius of gyration of the tracer cloud.
4.5. Alignment between vorticity and strain-rate tensor
Details of the errors in our measurements and a comparison between the dissipation rate estimated from the measured velocity gradient tensor and the velocity structure functions are reported in appendices A and B. To demonstrate the quality of our measurements, however, we briefly consider some of the well-known geometric properties of the velocity gradient – in particular, the alignment between the vorticity and the eigenvectors of the rate of strain. The rate of strain
$\unicode[STIX]{x1D61A}_{ij}$
is a symmetric second-rank tensor, and can be described by its three eigenvectors
$\hat{\boldsymbol{e}}_{i}$
(
$i=1,2,3$
), which correspond to the eigenvalues
${\it\lambda}_{i}$
with
${\it\lambda}_{1}\geqslant {\it\lambda}_{2}\geqslant {\it\lambda}_{3}$
. Intuitively, one would expect that the vorticity
$\hat{{\bf\omega}}$
would tend to align with the most extensional eigenvector,
$\hat{\boldsymbol{e}}_{1}$
. But, in general, vorticity is preferentially aligned with the intermediate eigenvector
$\hat{\boldsymbol{e}}_{2}$
(Ashurst et al.
Reference Ashurst, Kerstein, Kerr and Gibson1987). We show the probability distribution functions (PDFs) of the cosine of the angle between
$\hat{{\bf\omega}}$
and the
$\hat{\boldsymbol{e}}_{i}$
in figure 7 for our measured velocity gradients and for the direct numerical simulation. We choose to use the cosine of the angle instead of the angle itself because two vectors that are randomly oriented with respect to each other in three-dimensional space will have a uniform distribution for the cosine of the angle between them. In figure 7, the overall trend of the PDFs is very similar for the experiments and the simulation: as expected, the vorticity is best aligned with
$\hat{\boldsymbol{e}}_{2}$
. The trends, however, are more pronounced in the simulation due to experimental inaccuracies in the measurement of the velocity gradient, and potentially from the coarse-graining introduced by estimating the velocity gradient from a finite-sized tracer cloud.
Probability distribution function of the cosine of the angle between vorticity
$\hat{{\bf\omega}}$
and eigenvectors of the strain-rate tensor
$\hat{\boldsymbol{e}}_{i}$
: (a) experimental measurements, (b) simulation results.
5. Results and discussion
Figure 8 shows the alignment of rods with respect to the vorticity,
$\hat{{\bf\omega}}$
, and the three eigenvectors of the strain rate,
$\hat{\boldsymbol{e}}_{i}$
, for both our experiments and the simulation. For both cases, rods are much more strongly aligned with
$\hat{{\bf\omega}}$
than they are with
$\hat{\boldsymbol{e}}_{i}$
. This is because
$\hat{{\bf\omega}}$
and
$\hat{\boldsymbol{p}}$
each independently tend to align with the strongest Lagrangian stretching direction, as defined by the maximum eigenvector of the left Cauchy–Green strain tensor (Ni et al.
Reference Ni, Ouellette and Voth2014). The alignments of rods with the strain-rate eigenvectors are weaker because rods align with the Lagrangian stretching integrated over several Kolmogorov times, and this direction is usually quite different from the instantaneous stretching direction defined by the strain-rate eigenvectors.
Probability distribution function of the cosine of the angle between the rod orientation
$\hat{\boldsymbol{p}}$
and both vorticity
$\hat{{\bf\omega}}$
and eigenvectors of the strain-rate tensor
$\hat{\boldsymbol{e}}_{i}$
: (a) experimental measurements, (b) simulation results.
The alignment distributions from the simulation in figure 8 are in excellent agreement with previous simulations (Pumir & Wilkinson Reference Pumir and Wilkinson2011). Our experiment agrees quite well with the simulations, except that the experiments show that rods are slightly better aligned with
$\hat{\boldsymbol{e}}_{1}$
than they are with
$\hat{\boldsymbol{e}}_{2}$
, while the simulations show the opposite. Prior work on material lines (Wan Reference Wan2008) has shown that the relative strength of their alignment with
$\hat{\boldsymbol{e}}_{1}$
and
$\hat{\boldsymbol{e}}_{2}$
is sensitive to Reynolds number and perhaps to the particular driving of the flow. Material lines rotate just like thin rods with large aspect ratio
${\it\alpha}\gg 1$
. Thus, the small difference between the experiment and the simulation may result from differences in Reynolds number and the forcing mechanism.
In turbulence, the tumbling rate of rods has been found to be much smaller than it would be if the rods were randomly oriented (Shin & Koch Reference Shin and Koch2005; Parsa et al.
Reference Parsa, Calzavarini, Toschi and Voth2012). This decrease has been qualitatively understood to be a result of the alignment of rods with the vorticity vector. Given the alignment of rods with the vorticity, rods will preferentially rotate about their long axis, and therefore the contribution of
${\rm\Omega}_{ij}$
to the rod tumbling rate
$\dot{\boldsymbol{p}}$
will be weak.
To understand the rod tumbling rate in more detail, we can start with Jeffery’s equation to estimate the parameter space and evaluate the complexity of the problem. Recall that the velocity gradient tensor
$\unicode[STIX]{x1D63C}$
is a second-rank tensor with eight independent components, assuming incompressibility. For turbulent flow with isotropic small scales, the overall orientation in space does not matter so five scalars are sufficient to characterize
$\unicode[STIX]{x1D63C}$
(Meneveau Reference Meneveau2011). There are many possible choices for these five scalars, such as two of the eigenvalues of
$\unicode[STIX]{x1D64E}$
and the three components of the vorticity in the
$\unicode[STIX]{x1D64E}$
eigenframe, or five scalars constructed from moments of the velocity gradient tensor such as the well-known
$R=-(\unicode[STIX]{x1D608}_{im}\unicode[STIX]{x1D608}_{mn}\unicode[STIX]{x1D608}_{ni})/3$
and
$Q=-(\unicode[STIX]{x1D608}_{im}\unicode[STIX]{x1D608}_{mi})/2$
(Cantwell Reference Cantwell1992). Specifying the relative orientation of the rod in the
$\unicode[STIX]{x1D64E}$
eigenframe requires two additional independent angles. Thus, together with the five scalars required to determine
$\unicode[STIX]{x1D63C}$
, fully characterizing the rod tumbling rates requires a seven-dimensional parameter space. Denoting the generalized coordinates in this space by
$\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C})$
and the corresponding probability density function as
$\mathscr{P}(\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C}))$
, the mean-squared rod tumbling rate can be expressed as
$$\begin{eqnarray}\displaystyle \langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle & = & \displaystyle \int \text{d}\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C})\mathscr{P}(\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C}))\left[{\rm\Omega}_{ij}p_{j}+\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}(\unicode[STIX]{x1D61A}_{ij}p_{j}-p_{i}p_{k}\unicode[STIX]{x1D61A}_{kl}p_{l})\right]\nonumber\\ \displaystyle & & \displaystyle \times \,\left[{\rm\Omega}_{im}p_{m}+\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}(\unicode[STIX]{x1D61A}_{im}p_{m}-p_{i}p_{q}\unicode[STIX]{x1D61A}_{qn}p_{n})\right]\nonumber\\ \displaystyle & = & \displaystyle \int \left[{\rm\Omega}_{ij}{\rm\Omega}_{im}p_{j}p_{m}+\left(\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}\right)^{2}(\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{im}p_{j}p_{m}+p_{i}p_{k}\unicode[STIX]{x1D61A}_{kl}p_{l}p_{i}p_{q}\unicode[STIX]{x1D61A}_{qn}p_{n})+T\right]\nonumber\\ \displaystyle & & \displaystyle \times \,\text{d}\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C})\mathscr{P}(\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C})),\end{eqnarray}$$
where
$T$
represents six cross-terms between strain and rotation. If we assume that the rods are randomly oriented and are uncorrelated with the velocity gradient tensor, this equation simplifies considerably because some of the averages can be taken independently; for example,
$\langle \unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{im}p_{j}p_{m}\rangle =\langle \unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{im}\rangle \langle p_{j}p_{m}\rangle$
. This assumption leads to
$$\begin{eqnarray}\langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle =\frac{\langle {\rm\Omega}_{ij}{\rm\Omega}_{ij}\rangle }{3}+\frac{1}{5}\left(\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}\right)^{2}\langle \unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle .\end{eqnarray}$$
Given that, in isotropic turbulence,
$\langle \unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle =\langle {\rm\Omega}_{ij}{\rm\Omega}_{ij}\rangle =\langle {\it\epsilon}\rangle /{\it\nu}$
, the normalized rod tumbling rate would depend only on its eccentricity:
$$\begin{eqnarray}\frac{\langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle }{\langle {\it\epsilon}\rangle /{\it\nu}}=\frac{1}{6}+\frac{1}{10}\left(\frac{{\it\alpha}^{2}-1}{{\it\alpha}^{2}+1}\right)^{2}.\end{eqnarray}$$
However, in turbulence, rods are not randomly oriented, but are coupled with the velocity gradient tensor. Thus, characterizing their tumbling requires understanding of the full seven-dimensional PDF
$\mathscr{P}(\boldsymbol{X}(\,\boldsymbol{p},\unicode[STIX]{x1D63C}))$
. Since the Lagrangian velocity gradients are similar in many different turbulent flows, this PDF should be approximately universal. Experimentally, we do not have enough samples to obtain this PDF at a reasonable bin size, and it would be difficult to present even if we obtained it. The typical way to approach this complicated PDF would be to try to find a suitable low-dimensional projection to simplify the problem. We will show below that although the two most important dimensions are the magnitude of the strain and the enstrophy, the relative orientation of a rod in turbulence cannot be neglected in determining its tumbling rate.
Mean-squared tumbling rate of rods conditioned on the alignment between rod and vorticity
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle$
. The black dash-dotted line represents the mean-squared tumbling rate and the dashed line is from simulation results. Open symbols show the tumbling rate from differentiating the rod trajectories of orientation, and closed symbols show the calculation from Jeffery’s equation applied to the measured velocity gradient tensor.
Mean-squared tumbling rate of rods conditioned on (a) normalized enstrophy and (b) normalized dissipation rate. The symbols are the same as in figure 9.
Mean-squared tumbling rate of rods conditioned on the alignment between the rod and eigenvectors of the strain-rate tensor: (a)
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle$
; (b)
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{2}\rangle$
; (c)
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle$
. The symbols are the same as in figure 9.
In the rest of this paper, we will use the conditional average of the rod tumbling rate along different dimensions to characterize the importance of each dimension and to provide new experimental insights into rod tumbling dynamics. Recall that, as shown in figure 2, we have two different ways of experimentally determining the rod tumbling rate: a direct measurement,
${\dot{p}}_{i}$
, and a measurement inferred from Jeffery’s equation,
${\dot{p}}_{i}^{J}$
. To check for any potential systematic offset between these different measurements, we computed the mean-squared tumbling rates for each, finding
$\langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle =0.10\langle {\it\epsilon}\rangle /{\it\nu}$
and
$\langle {\dot{p}}_{i}^{J}{\dot{p}}_{i}^{J}\rangle =0.09\langle {\it\epsilon}\rangle /{\it\nu}$
for experiments and
$\langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle =\langle {\dot{p}}_{i}^{J}{\dot{p}}_{i}^{J}\rangle =0.09\langle {\it\epsilon}\rangle /{\it\nu}$
for simulations. Since
${\dot{p}}_{i}^{J}={\dot{p}}_{i}^{S}+{\dot{p}}_{i}^{{\rm\Omega}}$
, the mean-squared tumbling rate has three contributions,
$\langle {\dot{p}}_{i}^{S}{\dot{p}}_{i}^{S}\rangle$
,
$\langle {\dot{p}}_{i}^{{\rm\Omega}}{\dot{p}}_{i}^{{\rm\Omega}}\rangle$
and
$2\langle {\dot{p}}_{i}^{S}{\dot{p}}_{i}^{{\rm\Omega}}\rangle$
, which are respectively equal to 0.053, 0.069 and
$-0.039$
for experiments and 0.062, 0.073, and
$-0.045$
for simulations. To compare the experiment with simulation in detail, we systematically shifted the values of
${\dot{p}}_{i}^{J}$
and the tumbling rates computed from the simulation so that their mean-squared values were
$0.10\langle {\it\epsilon}\rangle /{\it\nu}$
. The shift was done by multiplying all data points in
${\dot{p}}_{i}^{J}$
and simulation results by 1.11, which will only move the curve up without changing its trend.
In figure 9, we plot the mean-squared rod tumbling rate conditioned on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle$
for
${\dot{p}}_{i}$
,
${\dot{p}}_{i}^{J}$
, and the simulation. All three curves nearly collapse, and show that the rod tumbling rate monotonically decreases by more than 50 % as its alignment with the vorticity increases. For
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle =0$
, the rod is almost perpendicular to the vorticity vector and its tumbling rate has a larger contribution from the vorticity than from the strain. If only the vorticity contributed, the mean-squared tumbling rate should be
$\langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle \approx \langle {\rm\Omega}_{ij}{\rm\Omega}_{ij}\rangle /3$
. Since
$\langle {\rm\Omega}_{ij}{\rm\Omega}_{ij}\rangle =\langle {\it\epsilon}\rangle /2{\it\nu}$
,
$\langle {\dot{p}}_{i}{\dot{p}}_{i}\rangle /(\langle {\it\epsilon}\rangle /{\it\nu})$
at
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle =0$
should be
$1/6$
, which is close to but slightly larger than the measured result of
${\sim}0.14$
. This discrepancy indicates that the strain contribution for
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle =0$
is on average in the direction opposed to the vortical motion, just as it would be for a rod in Jeffery orbits in a uniform shear flow.
In the other limit, for
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle =1$
, the rod is almost perfectly aligned with the vorticity, so its tumbling rate has no contribution from vortical motion. The tumbling rate does not vanish, however, due to the coupling to the strain field. The contribution from the strain
$\unicode[STIX]{x1D64E}$
can be investigated by considering the alignment of its three orthogonal eigenvectors
$\hat{\boldsymbol{e}}_{i}$
with
$\hat{\boldsymbol{p}}$
. As shown in figure 8, both simulation and experiment suggest that the rod tends to be aligned with
$\hat{\boldsymbol{e}}_{1}$
and
$\hat{\boldsymbol{e}}_{2}$
and perpendicular to
$\hat{\boldsymbol{e}}_{3}$
, indicating that the rod primarily lies in the plane formed by
$\hat{\boldsymbol{e}}_{1}$
and
$\hat{\boldsymbol{e}}_{2}$
.
Another way to assess the relative contribution to the rod tumbling rate of vortical motion and strain is by conditioning on the magnitude of
${\rm\Omega}_{ij}$
and
$\unicode[STIX]{x1D61A}_{ij}$
, which can be represented by the enstrophy
${\it\omega}^{2}={\it\omega}_{i}{\it\omega}_{i}$
and the dissipation rate
${\it\epsilon}=2{\it\nu}\langle \unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$
. Figure 10 shows the mean-squared rod tumbling rate conditioned on these two quantities. The tumbling rate monotonically increases by nearly two orders of magnitude with a similar dependence on both enstrophy and dissipation rate, indicating that the vorticity and the strain equally contribute to the tumbling rate of the rod on average. This two-decade change of the conditional rod tumbling rate is larger than other dimensions, which also suggests that the enstrophy and the dissipation rate are the most important two dimensions.
In addition to the effect of the strain magnitude, we also investigated the effect of the orientation between the rod and the eigenframe of the strain-rate tensor on the tumbling rate, as shown in figure 11. Recall that, in turbulence, the vorticity tends to be aligned with
$\hat{\boldsymbol{e}}_{2}$
(Ashurst et al.
Reference Ashurst, Kerstein, Kerr and Gibson1987). It would then be expected that the conditional average of the tumbling rate on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{2}\rangle$
(figure 11
b) would share some similarity with the conditional average on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle$
(figure 9). In both cases, the tumbling rate monotonically decreases with increasing alignment, but the trend is much steeper for the case of
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle$
, because the vorticity vector plays a stronger role in determining the rod tumbling. The conditional average of the tumbling rate on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle$
(figure 11
a) is more complicated to explain. If the rod were perfectly parallel or perpendicular to
$\hat{\boldsymbol{e}}_{1}$
, the strain would act to stretch or compress the rod rather than to rotate it, and the tumbling rate due to strain would be zero in both limits. But if a rod were oriented at 45
$^{\circ }$
with respect to
$\hat{\boldsymbol{e}}_{1}$
, the strain contribution to its tumbling rate would be maximized. So the squared tumbling rate should be small for both limits of alignment at
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle =0$
and
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle =1$
, but should have a peak near
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle =\cos (45^{\circ })=0.7$
. Our results in figure 11(a) are consistent with this expectation. The same argument can be made for the conditional average on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle$
, as shown in figure 11(c). In this case, however, the peak is much higher and the data are skewed towards
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle =1$
, as compared with the data conditioned on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle$
. Two factors are important here. First, rods are rarely aligned with
$\hat{\boldsymbol{e}}_{3}$
, as seen in figure 8. Events where this occurs are likely to be events with large vorticity, so the rotation rate will be large. In addition, the vorticity is preferentially perpendicular to
$\hat{\boldsymbol{e}}_{3}$
, and so a rod that is better aligned with
$\hat{\boldsymbol{e}}_{3}$
gains a larger contribution to its tumbling from vortical motion. Both of these effects move the peak up and towards
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle =1$
.
We have now described the dependence of rod tumbling rate on five independent dimensions: the magnitudes of the strain and the vorticity, captured by the dissipation rate and the enstrophy, respectively, and three independent characterizations of the rod orientation with respect to the velocity gradient tensor. All of them are important in determining the rod tumbling rate, and none can be neglected. The additional two dimensions necessary for determining the problem concern the orientation of vorticity in the strain eigenframe, and can be taken to be, for example,
$\langle \hat{{\bf\omega}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle$
and
$\langle \hat{{\bf\omega}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{2}\rangle$
. We find that the conditional mean-square tumbling rate has only a very weak dependence on these two dimensions. This may seem surprising since the relative orientation of vorticity and strain determine whether they reinforce or cancel each other as seen in figure 2. For prediction of the mean-square tumbling rate, the parameter space may be able to be simplified to five dimensions, but it seems that a complete picture of the rotation of rods in turbulence will require specifying the full seven-dimensional PDF of rod orientation and the velocity gradient tensor.
6. Conclusion
We have presented an experimental investigation of the tumbling rate of rods in turbulence. To assess the relative importance of various factors on the tumbling rate and to thereby estimate the effective dimensionality of the problem, we also simultaneously measured the full velocity gradient tensor near the rods. We obtained the gradient by implementing a scanning particle tracking system, allowing us to image a high concentration of tracer particles. The quality of our velocity gradient measurements is comparable to previous experiments, and we were able to measure the trajectories of anisotropic particles at the same time.
We carefully explored the mean-squared tumbling rate conditioned on several different variables, including
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle$
,
${\it\epsilon}$
,
${\it\omega}^{2}$
,
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle$
,
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{2}\rangle$
and
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle$
. These variables were chosen to give a framework within which the relative contributions to the rod tumbling rate are easily interpretable. We found that the mean-squared tumbling rate was dependent on all five of these dimensions to some degree. As these variables can be separated into two classes, those that depend on the strain rate and those that depend on the vorticity, our results suggest that the rod tumbling rate depends approximately equally on
$\unicode[STIX]{x1D61A}_{ij}$
and
${\rm\Omega}_{ij}$
.
We provide experimental evidence that rods are preferentially aligned with the vorticity, thus diminishing the potential contribution from
${\rm\Omega}_{ij}$
and partially explaining why the measured rod tumbling rates are smaller than they would be if the rod orientation were random. We also found that the mean-squared rod tumbling rate monotonically decreases with increasing
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{2}\rangle$
, where
$\hat{\boldsymbol{e}}_{2}$
is the intermediate strain-rate eigenvector, just as it does with
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{{\bf\omega}}\rangle$
. The contribution from the local strain is small when the rod is either aligned with or orthogonal to
$\hat{\boldsymbol{e}}_{1}$
, but is largest when the rod is oriented at roughly
$45^{\circ }$
with respect to
$\hat{\boldsymbol{e}}_{1}$
. The PDF of the rod tumbling rate conditioned on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle$
is similar to that conditioned on
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{1}\rangle$
; the dependence is stronger, however, for
$\langle \hat{\boldsymbol{p}}\boldsymbol{\cdot }\hat{\boldsymbol{e}}_{3}\rangle$
because the rod feels the full effect of both the strain and the vorticity simultaneously when it is aligned with
$\hat{\boldsymbol{e}}_{3}$
. This result suggests that the intermittent tumbling rate of rods may be linked back to the geometric information contained in the relative orientation of a rod with the local velocity gradient.
Joint probability distribution function of one of the diagonal components of the velocity gradient tensor versus the summation of the other two. The dashed line represents the case for incompressible fluid, i.e.
$\partial u_{i}/\partial x_{i}=0$
.
Finally, the new experimental technique we describe enables us to extract the flow motion near particles in a fully turbulent three-dimensional system, and has the potential to be applied to many other problems such as larger particles where larger particle Reynolds number makes simulations much more difficult or the flow field generated by active anisotropic particles such as bacteria.
Acknowledgements
We thank F. Toschi and E. Calzavarini for providing us with the DNS data. We acknowledge support from US National Science Foundation grants DMR-1206399 and DMS-1211952 to Yale University and DMR-1208990 to Wesleyan University and COST Actions MP0806 and FP1005.
Appendix A. Quality of experimental measurements
Here, we examine the quality of our experimental measurement of the velocity gradient tensor in several ways. First, we test the incompressible flow condition that requires the trace of the gradient to vanish. Figure 12(a–c) shows the joint probability density function (jPDF) of one of the diagonal components of the velocity gradient tensor versus the sum of the other two. The shape of the contour lines for all panels looks very similar, indicating that both the flow and the error are nearly isotropic. For perfectly incompressible flow, the vertical axes should be equivalent to the horizontal axes, as given by the dashed line. Deviations from the dashed line are caused by experimental uncertainties in measuring the velocity gradient. These deviations can be quantified using the ratio between major and minor axes of the contour ellipses in the jPDF. The larger this ratio is, the more accurate the measurement of the velocity gradient tensor. For all three panels in figure 12, this ratio is roughly 3.3, which is larger than the value of 2.6 in the early experiments of Lüthi et al. (Reference Lüthi, Tsinober and Kinzelbach2005). More recent experiments using scanning particle tracking systems by Hoyer et al. (Reference Hoyer, Holzner, Lüthi, Guala, Liberzon and Kinzelbach2005) and Krug et al. (Reference Krug, Holzner, Lüthi, Wolf, Tsinober and Kinzelbach2014) have measured deviations from incompressibility that are comparable or slightly smaller than ours.
Joint probability distribution function of the tumbling rate
${\dot{p}}_{i}$
from rod orientation measurements versus the tumbling rate calculated from Jeffery’s equation
${\dot{p}}_{i}^{J}$
applied to the measured velocity gradient tensor. The dashed line represents the case when two tumbling rates agree with each other, i.e.
${\dot{p}}_{i}={\dot{p}}_{i}^{J}$
.
The energy dissipation rate
${\it\epsilon}$
calculated from different methods: (a) comparison between
${\it\epsilon}$
obtained from the dissipative range of the longitudinal (
$D_{LL}$
, dot-dashed line) and transverse (
$D_{NN},solidline$
) second-order structure function and the direct measurement of the velocity gradient tensor (black dashed line). (b) Dissipation rate obtained from compensated third-order structure function (solid line) and the direct measurement of the velocity gradient tensor (black dashed line).
A more stringent and relevant test of our experimental accuracy is comparing the directly measured rod tumbling rate
${\dot{p}}_{i}$
computed by differentiating the time-dependent rod orientation with that calculated from Jeffery’s equation,
${\dot{p}}_{i}^{J}$
, using the measured velocity gradient tensor. Figure 13(a–c) shows the jPDF of
${\dot{p}}_{i}^{J}$
and
${\dot{p}}_{i}$
. The contours of the jPDF again are ellipsoidal and deviate from the linear dashed line, just as they did in figure 12. The ratio between the major and minor axes of the contours is also close to 3.3, indicating that the uncertainty in the measurements of the velocity gradient tensor is not amplified through Jeffery’s equation.
Previously, it has been proposed that the incompressibility condition can be used as a weighting factor for smoothing trajectories (Lüthi et al. Reference Lüthi, Tsinober and Kinzelbach2005). The apparent incompressibility is characterized by
$$\begin{eqnarray}{\it\delta}=\frac{\left|\displaystyle \frac{\partial u}{\partial x}+\displaystyle \frac{\partial v}{\partial y}+\displaystyle \frac{\partial w}{\partial z}\right|}{\left|\displaystyle \frac{\partial u}{\partial x}\right|+\left|\displaystyle \frac{\partial v}{\partial y}\right|+\left|\displaystyle \frac{\partial w}{\partial z}\right|}.\end{eqnarray}$$
The weighting factor is a function of this
${\it\delta}$
; it is designed to be 1 for small
${\it\delta}$
, and smoothly changes to zero for large
${\it\delta}$
. Thus if we smooth all the components of the velocity gradient using this weight, it nearly guarantees that the filtered velocity gradient has a small apparent compressibility. Note, however, that incompressibility only relies on the three diagonal components in the velocity gradient tensor, which are independent of the other six components; thus, compressibility may not be directly related to the accuracy of the full velocity gradient measurement. To test this argument, we plot in figure 14(a) the incompressibility test we used in figure 12 but only using samples with
${\it\delta}<0.15$
. This sampling is equivalent to applying a weighting function that equals 1 for
${\it\delta}<0.15$
and 0 for
${\it\delta}>0.15$
, similar to the one used before (Lüthi et al.
Reference Lüthi, Tsinober and Kinzelbach2005). After enforcing incompressibility, we would expect that the contours should lie closer to the linear dashed line. We also plot the jPDF of
${\dot{p}}_{i}^{J}$
and
${\dot{p}}_{i}$
in figure 14(b) as we did in figure 13, but again enforcing
${\it\delta}<0.15$
. In this case, the shape of the contours does not change much, suggesting that weighting by the measured compressibility does not significantly improve the overall accuracy of the full velocity gradient tensor measurements.
Appendix B. Dissipation rate measurement
From the velocity gradient tensor, we can directly obtain the mean energy dissipation rate
$\langle {\it\epsilon}\rangle =2{\it\nu}\langle \unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$
, one of the most important parameters in turbulence. For most experiments in homogeneous, isotropic turbulence without direct access to the velocity gradient tensor, the dissipation rate is extracted from the scaling of the second- or third-order Eulerian velocity structure functions in the inertial range. The velocity structure function of order
$n$
is defined to be the
$n\text{th}$
statistical moment of the velocity increment across a spatial scale
$r$
:
$\langle [\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})-\boldsymbol{u}(\boldsymbol{x})]^{n}\rangle$
. Applying Kolmogorov’s theory of isotropic turbulence (Kolmogorov Reference Kolmogorov1941), the longitudinal
$D_{LL}$
and transverse
$D_{NN}$
components of the second-order structure function (
$n=2$
) scale as
$D_{LL}=\langle {\it\epsilon}\rangle r^{2}/15{\it\nu}$
and
$D_{NN}=2\langle {\it\epsilon}\rangle r^{2}/15{\it\nu}$
in the dissipative range
$r<{\it\eta}$
. Figure 15(a) shows both
$D_{LL}$
and
$D_{NN}$
compensated with
$r^{2}/15{\it\nu}$
and
$2r^{2}/15{\it\nu}$
respectively. There is a very short plateau near
$r\approx {\it\eta}$
, and the two compensated structure functions collapse with each other well in that range. Thus, the value of that plateau tells us the dissipation rate, which in our experiment is
$\langle {\it\epsilon}\rangle \approx 9\times 10^{-5}~\text{m}^{2}~\text{s}^{-3}$
. In addition, figure 15(b) shows the third-order longitudinal structure function (
$D_{LLL}$
). The dissipation rate can be obtained from the Kolmogorov
$4/5$
law (
$D_{LLL}=-4\langle {\it\epsilon}\rangle r/5$
) in the inertial range
${\it\eta}\ll r\ll L$
, where
$L$
is integral length scale. Here, because of the finite Reynolds number and the small measurement volume, there is a sharp cutoff near the infrared end of the inertial range. There is therefore only a narrow peak roughly from 50
${\it\eta}$
to 60
${\it\eta}$
, but the estimate of
$\langle {\it\epsilon}\rangle$
from that peak is very close to the estimate from the dissipation-range scaling of the second-order structure functions. The dashed line in both panels shows the estimate from the velocity gradient measurement from tracer particles within a
${\sim}6{\it\eta}$
radius of the rod. It is clear that this value is slightly smaller than the estimate from both velocity structure function methods. This is because the velocity gradient is coarse-grained over a small volume, which removes contributions from some of the smallest scales to the energy dissipation rate.
