2. Dynamics of small inertial spheroids, time scales, and numerical methods

We focus on spheroidal particles, which are ellipsoids of revolution with two equal semi-axes $a=b$ and a principal axis $c = \lambda a$. The aspect ratio $\lambda$ characterises their shape: oblate particles have $\lambda < 1$, spheres $\lambda = 1$ and prolate particles $\lambda >1$ (see figure 1a,b). We consider the case where such particles are suspended in a developed, incompressible turbulent velocity field $\boldsymbol {u}(\boldsymbol {x},t)$ and are much smaller than the associated Kolmogorov dissipative scale, namely $a,b,c \ll \eta = \nu ^{3/4}/\varepsilon ^{1/4}$, where $\nu$ is the fluid kinematic viscosity and $\varepsilon$ is the mean dissipation rate of kinetic energy. Moreover, we assume that the particles are much heavier than the surrounding fluid, i.e. their mass density $\rho _{p}$ is much larger than the fluid density $\rho _{f}$, and that their velocity relative to the fluid, together with their size, defines an infinitely small Reynolds number. Finally, the particles are presumed to be sufficiently dilute to neglect their feedback on the flow and turbulence is presumed to be sufficiently intense to ignore the effects of gravity.

Figure 1.

(a) An oblate particle, $\lambda < 1$, together with a representative trajectory on which green arrows stand for particle orientation $\boldsymbol {e}_{\hat {z}}$ and black ones for its acceleration ${\rm d} \boldsymbol {V}/{\rm d} t$. (b) Same for a prolate particle, $\lambda > 1$. (c) Various response times (see text for definitions), normalised by the equal-mass time $\bar {\tau }=2\rho _{p}\lambda ^{2/3} a^2 /(9\rho _{f} \nu )$, as a function of the particle aspect ratio.

The dynamics of non-spherical particles involves both translational and rotational motions. The translation is determined by the position $\boldsymbol {X}$ and velocity $\boldsymbol {V}$ of the particle's centre of mass in the inertial frame of reference $\mathcal {R} = (\boldsymbol {e}_{x}, \boldsymbol {e}_{y}, \boldsymbol {e}_{z})$ of the fluid flow. Under the above assumptions, they are given by the linear momentum equations (Brenner Reference Brenner1964)

(2.1a,b)\begin{equation} \frac{{\rm d} \boldsymbol{X}}{{\rm d} t} = \boldsymbol{V}, \quad \frac{{\rm d}\boldsymbol{V}}{{\rm d} t} ={-}\frac{1}{\bar{\tau}}\boldsymbol{\mathsf{A}}^{\mathrm{T}}\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{A}} \left[\boldsymbol{V} - \boldsymbol{u}(\boldsymbol{X},t)\right], \end{equation}

where we have introduced the response time $\bar {\tau }=2\rho _{p}\lambda ^{2/3} a^2 /(9\rho _{f} \nu )$ associated with a spherical particle with the same mass as the spheroid. Here $\boldsymbol{\mathsf{A}}$ denotes the rotation matrix that maps $\mathcal {R}$ to the reference frame $\widehat {\mathcal {R}} = ({\boldsymbol {e}}_{\hat {x}}, {\boldsymbol {e}}_{\hat {y}}, {\boldsymbol {e}}_{\hat {z}})$ of the particle, ${\boldsymbol {e}}_{\hat {z}}$ being along its revolution axis (see figure 1a,b). The drag tensor is expressed in frame $\widehat {\mathcal {R}}$ where it is diagonal, $\boldsymbol{\mathsf{D}}= \mathrm {diag}[{\mathsf{D}}_\perp, {\mathsf{D}}_\perp, {\mathsf{D}}_\parallel ]$, with

(2.2a,b)\begin{equation} {\mathsf{D}}_\perp= \frac{8 (\lambda ^2-1)}{3\lambda^{1/3}[\chi(2\lambda^2-3)+\lambda]},\quad {\mathsf{D}}_\parallel=\frac{4 (\lambda ^2-1)}{3\lambda^{1/3}[\chi(2 \lambda ^2-1)-\lambda]} \end{equation}

and

(2.2c)\begin{equation} \chi= \frac{\log(\lambda+\sqrt{\lambda^2-1})}{\sqrt{\lambda^2-1}}. \end{equation}

For detailed derivations of these drag coefficients, we refer the reader to Brenner (Reference Brenner1964). Note that when $\lambda <1$, the factor $\chi$ involves the logarithm of a unit-modulus complex number and reduces to $\chi =\arctan (\sqrt {1-\lambda ^2}/\lambda )/\sqrt {1-\lambda ^2}$, as usually stipulated for oblate particles.

The components of $\boldsymbol{\mathsf{D}}$ define two shape-dependent time scales, $\tau _\perp =\bar {\tau }/{\mathsf{D}}_\perp$ and $\tau _\parallel =\bar {\tau }/{\mathsf{D}}_\parallel$, associated with the drag in the directions perpendicular and parallel to the particle axis of symmetry, respectively. Shapiro & Goldenberg (Reference Shapiro and Goldenberg1993) introduced an effective response time $\hat {\tau }$ by performing an isotropic average of the mobility tensor $\bar {\tau } \boldsymbol{\mathsf{D}}^{-1}$ over all orientations:

(2.3)\begin{equation} \hat{\tau} = \frac{1}{3}\,\mathrm{tr}(\boldsymbol{\mathsf{D}}^{{-}1})\bar{\tau} = \frac{2\tau_\perp+\tau_\parallel}{3} = \bar{\tau}\lambda^{1/3}\chi = \frac{2\rho_{p}\lambda a^2}{9\rho_{f}\nu}\,\frac{\log(\lambda +\sqrt{\lambda^2-1})}{\sqrt{\lambda^2-1}}. \end{equation}

Fan & Ahmadi (Reference Fan and Ahmadi1995) considered an alternative effective response time $\check {\tau }$ obtained by averaging the drag tensor rather than the mobility matrix. This corresponds to the harmonic mean of the orientation-dependent response times, namely

(2.4)\begin{align} \check{\tau}=\left(\frac{1}{3}\,\mathrm{tr}(\boldsymbol{\mathsf{D}})\right)^{{-}1}\bar{\tau} = \left(\frac{2\tau_\perp^{{-}1}+\tau_\parallel^{{-}1}}{3}\right)^{{-}1} = \frac{\rho_{p}\lambda a^2}{2\rho_{f} \nu}\,\frac{4\chi^2(\lambda^2-1)^2 - (\lambda-\chi)^2}{(\lambda ^2-1)[10 \chi(\lambda^2-1)-3(\lambda-\chi)]}. \end{align}

The dependence of the two response times $\hat {\tau }$ and $\check {\tau }$ upon the aspect ratio $\lambda$ is shown in figure 1(c). Surprisingly, their difference hardly exceeds a couple of per cent. The relative discrepancy $\varDelta =(\hat {\tau }-\check {\tau })/\bar {\tau }$ between these two time scales can be obtained by asymptotic expansion in different limits. Specifically, we find $\varDelta \simeq (2/225)(\lambda -1)^2$ for nearly spherical particles ($\lambda \simeq 1$), $\varDelta \simeq ({\rm \pi} /56)\lambda ^{1/3}$ for thin disks ($\lambda \ll 1$), and $\varDelta \simeq (1/10)(\log \lambda )\lambda ^{-2/3}$ for slender fibres ($\lambda \gg 1$). Such tiny differences make it almost impossible to select the most relevant of these two response times from numerical or experimental data. Nevertheless, as observed in figure 1(c), the difference between these two time scales and the perpendicular and transverse response times becomes more evident and exceeds $10\,\%$ when there is a significant deviation of $\lambda$ from 1. This substantial dependence on $\lambda$ and $\bar {\tau }$ enables us to differentiate $\hat {\tau }$ and $\check {\tau }$ from other possible definitions of an effective particle response time.

Turning now to rotational dynamics, the orientation matrix $\boldsymbol{\mathsf{A}}$ evolves with an angular velocity $\boldsymbol {\varOmega }$, which is more conveniently expressed in the particle reference frame $\widehat {\mathcal {R}}$, so that

(2.5a,b)\begin{equation} \frac{{\rm d} \boldsymbol{\mathsf{A}}}{{\rm d} t} ={-}\begin{bmatrix} \,0 & -\varOmega_{\hat{z}} & \phantom{-}\varOmega_{\hat{y}}\\ \phantom{-}\varOmega_{\hat{z}} & \,0 & -\varOmega_{\hat{x}}\\ -\varOmega_{\hat{y}} & \phantom{-}\varOmega_{\hat{x}} & \,0 \end{bmatrix} \boldsymbol{\mathsf{A}}, \quad \frac{{\rm d} (\boldsymbol{\mathsf{I}} \boldsymbol{\varOmega})}{{\rm d} t} + \boldsymbol{\varOmega}\times(\boldsymbol{\mathsf{I}} \boldsymbol{\varOmega})= \boldsymbol{T}, \end{equation}

where $\boldsymbol{\mathsf{I}}= (4/15){\rm \pi} \rho _{p}\lambda a^5 \, \mathrm {diag}[1+\lambda ^2,1+\lambda ^2,2]$ is the spheroid's moment of inertia about its principal axis and $\boldsymbol {T}$ the hydrodynamic torque acting on the particle. The latter reads (Jeffery Reference Jeffery1922)

(2.6)\begin{equation} \boldsymbol{T} = \left[ \begin{array}{@{}l@{}} T_{\hat{x}}\\ T_{\hat{y}}\\ T_{\hat{z}} \end{array}\right] =\frac{16}{3}{\rm \pi}\rho_{f}\,\nu\,a^3 \left[ \begin{array}{@{}c@{}} -\dfrac{\lambda^4-1}{\chi(2\lambda^2-1)-\lambda} \left(\varOmega_{\hat{x}}-\dfrac{1}{2}\omega_{\hat{x}}+\dfrac{\lambda ^2-1}{\lambda^2+1}{\mathsf{S}}_{\hat{z}\hat{y}}\right)\\ -\dfrac{\lambda^4-1}{\chi(2\lambda^2-1)-\lambda} \left(\varOmega_{\hat{y}}-\dfrac{1}{2}\omega_{\hat{y}}-\dfrac{\lambda ^2-1}{\lambda^2+1}{\mathsf{S}}_{\hat{x}\hat{z}}\right)\\ -\dfrac{\lambda ^2-1}{\lambda -\chi}\left(\varOmega_{\hat{z}}-\dfrac{1}{2}\omega_{\hat{z}}\right) \end{array}\right], \end{equation}

and depends on the fluid vorticity $\boldsymbol {\omega } = \boldsymbol {\nabla }\times \boldsymbol {u}$ and strain tensor $\boldsymbol{\mathsf{S}} =(\boldsymbol {\nabla } \boldsymbol {u} +\boldsymbol {\nabla } \boldsymbol {u}^{\mathrm {T}})/2$, both evaluated at the particle position and in the rotating reference frame $\widehat {\mathcal {R}}$. Note that the fluid velocity gradient sampled by the particle is itself determined by the dynamics (2.1a,b), resulting in an intricate coupling between translational and rotational motions through the fluid flow.

Equations (2.5a,b) and (2.6) define two rotational response times associated with the particle tumbling (rotation along the $\hat {x}$ and $\hat {y}$ axes) and with its spinning (rotation along the spheroid's axis of symmetry $\hat {z}$). They correspond to the ratio between the components of the moment of inertia and the shape-dependent coefficient appearing in Jeffery's torque and read

(2.7a,b)\begin{align} \tau_{tumb} = \frac{\rho_{p}\lambda a^2}{20\rho_{f}\nu}\,\frac{\chi(2 \lambda^2-1)-\lambda}{\lambda^2-1}= \frac{3}{10}\tau_\parallel, \quad \tau_{spin} = \frac{\rho_{p}\lambda a^2}{10\rho_{f} \nu}\, \frac{\lambda-\chi}{\lambda^2-1} =\frac{3}{10}(2\tau_\perp{-}\tau_\parallel). \end{align}

These two times are displayed in figure 1(c). As stressed by Zhao et al. (Reference Zhao, Challabotla, Andersson and Variano2015) and Marchioli, Zhao & Andersson (Reference Marchioli, Zhao and Andersson2016), they remain shorter than the translational time scales for all values of the aspect ratio. More precisely, one can actually show that both the tumbling and the spinning response times are smaller than $(9/20)\hat {\tau }$ for all $\lambda >0$, this bound being attained by $\tau _{tumb}$ for $\lambda \to \infty$ and by $\tau _{spin}$ for $\lambda \to 0$. The separation of time scales between the translational and rotational motions of the particle is hence clear, but one can question whether a factor ${\approx }2$ is enough to assume complete decoupling. To address this, we have recourse to DNS and investigate shape dependence in the statistics of heavy spheroids transported by a homogeneous isotropic turbulent flow.

The three-dimensional incompressible Navier–Stokes equations with large-scale forcing are integrated using the parallel pseudospectral solver LaTu with third-order Runge–Kutta time marching (see Homann, Dreher & Grauer Reference Homann, Dreher and Grauer2007). Relevant simulation parameters are summarised in table 1. Equations (2.1a,b) and (2.5a,b) for the dynamics of spheroidal particles are integrated numerically using an exponential integrator with the same time stepping as the fluid flow. The fluid velocity and its gradient are obtained at the particle position by tricubic interpolation of the corresponding fields from the Eulerian grid. The orientation of each individual particle is represented as a quaternion (see Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008; Siewert et al. Reference Siewert, Kunnen, Meinke and Schröder2014a) to ease numerical integration and stability. In order to span the particles’ parameter space, we have considered $90$ different families of $500\,000$ particles each. They combine nine different aspect ratios, ranging from $\lambda = 0.1$ to $10$, and ten response times, spanning ${\bar {\textit {St}}}= \bar {\tau }/\tau _\eta = 0.1$ to $6.4$, where $\tau _\eta =(\nu /\varepsilon )^{1/2}$ denotes the Kolmogorov dissipative time scale. The properties of spheroids are stored with a period ${\approx }6.8\tau _\eta$. Statistics are computed over approximately four large-eddy turnover times after a statistical steady state is reached.

Table 1.

The DNS parameters: $N^3$, number of collocation points; $\nu$, kinematic viscosity; $\Delta t$, time step; $\varepsilon$, average dissipation rate; $\eta = (\nu ^{3}/\varepsilon )^{1/4}$, Kolmogorov dissipative scale; $\tau _\eta = (\nu /\varepsilon )^{1/2}$, Kolmogorov time; $u_{rms}$, root-mean-square velocity; $L = u_{rms}^3/\varepsilon$, large scale; $\tau _L = L/ u_{rms}$, large-eddy turnover time; and ${Re}_\lambda = \sqrt {15}\, u_{rms}^2/(\nu \varepsilon )^{1/2}$, Taylor-scale Reynolds number.

4. Concluding remarks

In this study, we investigated the dynamics of heavy spheroidal particles transported by a homogeneous isotropic turbulent flow. Our numerical simulations provided evidence that the isotropically averaged Stokes number, as introduced by Shapiro & Goldenberg (Reference Shapiro and Goldenberg1993), and based on the effective response time obtained by averaging the particle mobility tensor over all possible orientations, captures the shape dependence in their translational dynamics. Whether this definition is more relevant than that proposed by Fan & Ahmadi (Reference Fan and Ahmadi1995), where the drag tensor is averaged instead, remains to be determined by future work with much more precise statistics. Nevertheless, our results demonstrate that the single-particle statistics of an inertial spheroid, including the probability distributions of velocity and acceleration and their time correlations, are similar to those of an equivalent spherical particle whose effective radius can be written as $a_{eff} = (\lambda \chi )^{1/2} a$, where the dependence of $\chi$ on the aspect ratio $\lambda$ is given in (2.2c).

This paves the way for the development of macroscopic models for turbulent transport of non-spherical particles, in which particles can be approximated as effective spheres, eliminating the need to consider the intricacies of their angular motion. This simplified approach holds significant potential for the various applications mentioned in the introduction. Our observations suggest that angular and translational dynamics are only weakly correlated, which can be explained by the fast time scales associated with particle spinning and tumbling. Specifically, the response to fluid-flow rotation is at least twice as fast as translational equilibration.

While our study has focused on statistically isotropic situations, it is worth noting that intricate relations between rotation and translation could arise when anisotropies are significant. A first instance is when particle gravitational settling is of the same order as turbulent motions. The preferential orientation of particles during their fall that sets their settling speed could compete with orientation fluctuations induced by turbulence. A second instance is flow with a mean shear or in the presence of boundaries. In these cases, prolate particles tend to orient in the direction of the flow, whereas oblate particles are more likely to align in the direction of its gradient. Turbulent structures in such flows also exhibit strong anisotropies, and the response of spheroids to turbulent fluctuations may depend non-trivially on their shape. Recent experimental results by Baker & Coletti (Reference Baker and Coletti2022) suggest that, near to the walls of a channel flow, rods tumble more frequently than disks, and the latter respond more slowly to fluid-velocity fluctuations. These observations imply that efficient macroscopic models for the transport of spheroids by anisotropic flow may need to weigh differently the components of the particle mobility tensor. Further work is needed to develop a comprehensive understanding of the effects of anisotropies on the turbulent transport of non-spherical particles.

Finally, many questions remain unanswered regarding two-particle statistics of spheroids in turbulence. While we observed a decoupling between translation and rotation at the single-particle level, this may not hold true for their relative motion. Previous studies have shown that, in the absence of inertia, spheroidal particles tend to align with the eigendirections of the Cauchy–Green tensor (Ni et al. Reference Ni, Ouellette and Voth2014), and we expect these correlations to persist for low-inertia particles. In this case, particles concentrate on dynamically evolving attractors with a fractal structure related to the stretching and compression directions of the Cauchy–Green tensor. This could lead to intricate relationships between clustering and alignment, possibly making it impossible to describe particle spatial patterns in terms of a single shape-dependent Stokes number. For particles with large inertia, the presence of caustics, where particles with very different histories come arbitrarily close to each other, competes with fractal clustering. Such particles are likely to have significantly different orientations, making it even more challenging to predict how their relative motion depends on their shape.