4.1.1. Angular dynamics near the channel centre

To qualitatively explain why the distributions in figure 3 have power-law tails we consider a two-dimensional toy model for the angular dynamics (Zhao et al. Reference Zhao, Gustavsson, Ni, Kramel, Voth, Andersson and Mehlig2019). In two dimensions, the left Cauchy–Green tensor has two eigenvectors, the expanding eigenvector $\hat {\boldsymbol {e}}_{{L}1}$ and the contracting eigenvector $\hat {\boldsymbol {e}}_{{L}3}$, but not the intermediate eigenvector $\hat {\boldsymbol {e}}_{{L}2}$. This means that the two-dimensional model may explain the dynamics of the angle $\alpha$ but not that of the angle $\beta$, see figure 1. Following Gustavsson & Mehlig (Reference Gustavsson and Mehlig2016) we model the homogeneous and isotropic fluid-velocity fluctuations as Gaussian random functions that are white in time, but have smooth spatial correlations. In two dimensions, (4.1) can be written in terms of two angles: $\phi$, the angle that the Lagrangian stretching direction makes with the $x$-axis in the channel coordinates, and $\alpha \approx |\boldsymbol {\delta n}|, |\boldsymbol {\delta n}| \ll 1$, the angular separation between the particle and the Lagrangian stretching direction. The angular separation between the Lagrangian stretching direction and the symmetry vector $\boldsymbol {n}$ of a particle with shape parameter $\Lambda = 1-\delta \Lambda$ is simply given by $\alpha = \phi (\Lambda =1)-\phi (\Lambda =1-\delta \Lambda )$. Here $\phi _{\Lambda }$ follows Jeffery's equation (2.2), for $\boldsymbol {n} = [\cos \phi _\Lambda , \sin \phi _\Lambda ]^{\textsf T}$. In the following we drop the subscript in $\phi _{\Lambda =1}$.

We start off by assuming the fluid velocity-gradient matrix to be a Gaussian random variable with time correlation $\tau$. The dimensional parameters of the problem are the strength of velocity fluctuations, $u_0$, the correlation length $\eta$ and the correlation time $\tau$. Thus there are two relevant time scales, the correlation time $\tau$ and the advection time $\eta /u_0$. Out of these one can make one dimensionless parameter ${Ku} = u_0 \tau / \eta$ (Duncan et al. Reference Duncan, Mehlig, Östlund and Wilkinson2005). We non-dimensionalise the fluid-velocity gradient as ${\boldsymbol{A}} = ({{Ku}^2}/{\tau }) {\boldsymbol{A}}^\prime$ and the time as $t= ({\tau }/{{Ku}^2}) t^\prime$ and drop the primes in ${\boldsymbol{A}}, t$ in the following. Assuming tracelessness, and isotropy for the fluid gradient matrix ${\boldsymbol{A}}$, one finds that ${\boldsymbol{A}}$ has three independent components $O_{12}, S_{11},S_{12}$. Then for the variables $\phi$ and $\alpha$, (4.1) can be written as

(4.2a)\begin{gather} \dot{\phi} = - O_{12} -\sin 2 \phi S_{11} + \cos 2 \phi S_{12} , \end{gather}

(4.2b)\begin{gather}\dot{\alpha} = -(2 \alpha \cos 2 \phi + \delta \Lambda \sin 2 \phi) S_{11} + (- 2 \alpha \sin 2 \phi + \delta \Lambda \cos 2 \phi) S_{12}. \end{gather}

The white noise limit is taken as ${Ku} \to 0$, which corresponds to assuming that the fluid velocity correlation time is the shortest time scale in the problem, $\tau \ll ({\eta }/{u_0})$. We find that the drift coefficients vanish. The diffusion coefficients are given by

(4.3)\begin{equation} \left.\begin{array}{c@{}} \mathcal{D}_{\alpha \phi} =\mathcal{D}_{\phi \alpha}= \dfrac{1}{2} \delta\Lambda , \\ \mathcal{D}_{\phi \phi} = \tfrac{3}{2}, \\ \mathcal{D}_{\alpha \alpha} = \dfrac{1}{2} (\delta \Lambda^2 + 4 \alpha^2). \end{array}\right\} \end{equation}

The stationary Fokker–Planck equation for the joint distribution $P(\alpha ,\phi )$ is

(4.4)\begin{equation} \left[ \frac{3}{2}\frac{\partial^2}{\partial \phi^2} + \delta\Lambda \frac{\partial^2}{\partial \phi \ \partial \alpha} +\frac{\partial^2}{\partial \alpha^2} \frac{\delta\Lambda^2 + 4 \alpha^2}{2} \right]P(\alpha,\phi) = 0. \end{equation}

We obtain the marginal distribution $P(\alpha )$ by integrating out $\phi$ from (4.4), requiring symmetry $P(\alpha ) = P(-\alpha )$, and that $P(\phi )$ is normalised to unity. The result is

(4.5)\begin{equation} P(\alpha) = \frac{\delta \Lambda }{[(\delta \Lambda) ^2+4 \alpha ^2 ] \tan ^{-1}\left(\dfrac{{\rm \pi} }{\delta \Lambda }\right)}. \end{equation}

Equation (4.5) captures the qualitative features mentioned before, namely the power-law form of the distribution and its cutoff at small angles, of the order of $\delta \Lambda$. In the regime $\delta \Lambda \ll |\alpha | < {\rm \pi}/2$, where $\alpha$ is much larger than $\delta \Lambda$, the relative-angle distribution exhibits power-law tails. By contrast, in the regime $\alpha <\delta \Lambda$, the relative-angle distribution shows a plateau. Equation (4.2b) shows that at small angles $\alpha < \delta \Lambda$, $\dot {\alpha }$ is dominated by the additive term, which is the term independent of $\alpha$. This gives rise to the plateau in the distribution for $\alpha$, indicating random uncorrelated motion of $\phi _\Lambda$ and the Lagrangian stretching direction.

Using (4.5) we find that $\langle \alpha ^2 \rangle =\frac 12 \delta \Lambda$ for $\delta \Lambda \ll {\rm \pi}$, so that the variance of the relative angle is proportional to the cutoff angle in the distribution of $\alpha$. The two-dimensional toy model predicts a linear dependence of the cutoff angle on $\delta \Lambda$. Figure 6(b) shows that the cutoff angle $\alpha _c$ depends linearly on $\delta \Lambda$ also in DNS of turbulent channel flow.

Figure 6.

(a) Mean time between tumbles $\langle \tau \rangle s$ as a function of $\delta \Lambda = 1-\Lambda$ for a set of trajectories in the layer $0<z^+ <10$, $s$ is the mean shear along the observed trajectories. Symbols correspond to simulation results, the dashed line corresponds to theory. (b) Here, $\alpha _{c}$ as a function of $\delta \Lambda$ for $z^+\sim 4$ (red triangles) and $z^+\sim 180$ (orange squares) where $\alpha _c$ describes the transition between the power law and the plateau in figure 3(a) in orange, and figure 4(a) in red. Blue circles and green diamonds correspond to $z^+ \sim 10, 20$, respectively. The black dashed line shows the reference slope $\delta \Lambda ^1$, as predicted by theory.

Geometrically, the power law in $\alpha$ can be understood as a consequence of the fact that the Lagrangian stretching direction $e_{L1}(\boldsymbol {x}(t),t)$ acts as an attractor for the orientation field of particles with shape parameter $\Lambda$, $\phi _{\Lambda }(\boldsymbol {x}(t),t)$ along the same Lagrangian trajectory $\boldsymbol {x}(t)$. This clustering of orientations is analogous to spatial clustering of particles in turbulence in the advective limit (Gustavsson et al. Reference Gustavsson, Berglund, Jonsson and Mehlig2016; Meibohm et al. Reference Meibohm, Pistone, Gustavsson and Mehlig2017), and correlated random walks in the inertia free limit (Dubey et al. Reference Dubey, Meibohm, Gustavsson and Mehlig2018). Moreover, just as in the case of advected particles in one dimension, the power law in $\alpha$ is a result of purely diffusive dynamics with a diffusion coefficient proportional to $\alpha ^2$, leading to the same power-law exponent.

In the model the power-law exponent equals $-2$, and it is independent of $\delta \Lambda$. The exponent is very nearly $-2$ for the DNS results but since the theoretical model we have considered is two-dimensional, we don't expect it to explain the exponent, but merely the mechanism. That the model agrees only qualitatively with DNS results is to be expected, the model neither describes the additional degree of freedom $\beta$, nor does it capture the persistent nature of the turbulent velocity-gradient fluctuations. But numerical simulations of a three-dimensional model show qualitatively similar results.

4.1.2. Angular dynamics near the channel wall

Near the channel wall the fluid-velocity gradient has a large shear component. Since the flow is anisotropic it is natural to express the angular dynamics in terms of the channel-fixed basis $\hat {\boldsymbol {x}}$, $\hat {\boldsymbol {y}}$ and $\hat {\boldsymbol {z}}$ with the corresponding Euler angles $\theta$ and $\phi$, figure 1(b).

We consider a simple model for a spheroidal particle experiencing simple shear and additive noise. Turitsyn (Reference Turitsyn2007) used this model to study the angular dynamics of single polymers in a chaotic flow with mean shear. While Turitsyn (Reference Turitsyn2007) assumed $\Lambda =1$, here we consider general values of the shape parameter $\Lambda$.

Decompose the fluid gradient matrix ${\boldsymbol{A}}$ as a sum of the mean $\ \bar {\!\!{\boldsymbol{A}}}$ and the fluctuations ${\boldsymbol{A}}^\prime$ as ${\boldsymbol{A}} = \,\bar {\!\!{\boldsymbol{A}}} + {\boldsymbol{A}}^\prime$. The only non-zero component of $\ \bar {\!\!{\boldsymbol{A}}}$ is $\bar {A}_{xz} = s$ (Challabotla et al. Reference Challabotla, Zhao and Andersson2015). Jeffery's equation (2.2) for the angles $\phi ,\theta$ in terms of $s$ and ${\boldsymbol{A}}^\prime$ is given by

(4.6)\begin{gather} \dot{\phi} = -\frac{s}{2} (1 -\Lambda \cos 2 \phi) + \eta_\phi, \end{gather}

(4.7)\begin{gather}\dot{\theta} = -\Lambda \frac{s}{4} \sin 2 \phi \sin 2 \theta + \eta_\theta. \end{gather}

Here $\eta _\phi$ and $\eta _\theta$ are fluctuating terms which in general depend on ${\boldsymbol{A}}^\prime$, the angles $\theta , \phi$ and the shape parameter $\Lambda$. In the vicinity of the channel wall the shear strength $s$ is much larger than the magnitude of fluctuations of the velocity-gradient matrix, and thus the particle spends long times close to $\phi = 0$. Since we have $\eta _\phi = A_{zx} + O(\phi )+ O(\theta )$ in the limit $\phi ,\theta \to 0$, one can approximate $\eta _\phi \approx A_{zx}$. Let the autocorrelation for $\eta _\phi$ be given by $\langle \eta _\phi (0)\eta _\phi (t)\rangle = C_0 \ f(|t|/\tau )$, so that $C_0$ quantifies the magnitude of the fluctuations and $\tau$ is the correlation time. Then, (4.6) has three time scales, $1/s, 1/\sqrt {C_0}$ and $\tau$. Near the wall the shear strength is much larger than the strength of the fluctuations, $s \gg \sqrt {C_0}$. This gives that the dynamics of $\phi$ is comprised of two regions, the deterministic fast region with time scale $1/s$ and the stochastic slow region with time scale $1/\sqrt {C_0}$. Observation of typical trajectories show that near the channel wall the fluctuations of $\eta _\phi$ are much faster than the time scale of the slow, stochastic $\phi$ dynamics. Thus we take the white noise limit as $\tau \sqrt {C_0} \to 0$, while holding $2 D= \int _{-\infty }^\infty \, \textrm {d} t \eta _\phi (0) \eta _\phi (t) \propto C_0 \tau$ constant.

Figure 2 appears to show a marked difference between the dynamics of the Lagrangian stretching direction and a rod-like particle with $\Lambda = 0.9963$. In particular one observes that near the wall, the Lagrangian stretching direction seems to align along $\hat {\boldsymbol {x}}$ whereas $\Lambda =0.9963$ seems to tumble along the same trajectory. This can be explained by calculating the mean time between tumbles, defined as the average time taken by the particle to travel from $\phi =+{\rm \pi} /2$ to $\phi =-{\rm \pi} /2$, figure 6(a). A theoretical calculation using (4.6), see appendix A, shows good quantitative agreement with simulations. This is expected, considering that the parameter $\sqrt {C_0} \tau \sim 10^{-3}$ for the trajectory set used to obtain figure 6(a), so that the white noise limit is a good approximation.

Next we analyse the dynamics for $\phi$. This is because the $\phi$ dynamics are independent of $\theta$, whereas the $\theta$ dynamics are slave to the process $\phi$.

Assuming $\eta _\phi$ to be a Gaussian random variable, white in time, with the intensity of fluctuations $2 D = \int _{-\infty }^\infty \, \textrm {d} t \eta _\phi (0) \eta _\phi (t)$, Turitsyn (Reference Turitsyn2007) obtained the Fokker–Planck equation for the distribution of $\phi$,

(4.8)\begin{equation} \frac{\partial P(\phi,t)}{\partial t} = \frac{s}{2}\frac{\partial }{\partial \phi}[ (1 -\Lambda \cos 2 \phi)P(\phi,t)] + D \frac{\partial^2}{\partial \phi^2} P(\phi,t). \end{equation}

Equation (4.8) differs from (9) in Turitsyn (Reference Turitsyn2007) because we consider general values of $\Lambda$ and not just $\Lambda =1$, and that the strength of fluctuations $D$ is defined slightly differently. Following Turitsyn (Reference Turitsyn2007), the steady-state distribution $P(\phi )$ is obtained as the time independent solution of (4.8),

(4.9)\begin{equation} P(\phi) = \mathcal{N} \int_0^{\rm \pi} \mathrm{d}x \exp-\frac{s}{2 D}(x -\Lambda \cos(2 \phi -x)\sin x). \end{equation}

Here one integration constant is determined by periodic boundary conditions, $P({\rm \pi} /2) = P(-{\rm \pi} /2)$, and $\mathcal {N}$ is a normalisation constant which can be computed using $\int _{-{\rm \pi} /2}^{{\rm \pi} /2} \mathrm {d} \phi P(\phi ) =1$.

In order to understand how the relative angle between a slender rod and the Lagrangian stretching direction behaves, next we rewrite (4.1) for the relative angle, $\delta \phi = \phi (\Lambda =1)-\phi (\Lambda =1-\delta \Lambda )$. The joint equations of motion for $\phi , \delta \phi$ up to the second order in $\delta \phi$ are given by,

(4.10)\begin{gather} \dot{\phi} = -\frac{s}{2} (1 -\cos 2 \phi) + \eta_\phi, \end{gather}

(4.11)\begin{gather}\dot{\delta\phi} = \frac{s}{2} \delta\Lambda \cos 2 \phi - (s \Lambda \sin 2\phi ) \delta\phi + (s \Lambda \cos 2\phi ) \delta\phi^2. \end{gather}

First assume $\phi$ takes its mean value, $\phi _0 \approx ({\sqrt {{\rm \pi} }}/{\Gamma (1/6)})(\frac {3}{2})^{{1}/{3}} (D/s)^{{1}/{3}} \ll 1$. Then we have $\dot {\delta \phi } \approx ({s}/{2}) \delta \Lambda - (s \Lambda 2\phi _0 ) \delta \phi + (s \Lambda ) \delta \phi ^2$ to the first order in $\phi _0$. This equation has fixed points $\phi _0 \pm \phi _0 \sqrt {1-{\delta \Lambda }/{2\Lambda \phi _0^2}}$. For $\delta \Lambda \ll 2 \Lambda \phi _0^2$ there is a stable fixed point at $\delta \phi _c \approx {\delta \Lambda }/{4 \Lambda \phi _0}$ and an unstable fixed point at $2 \phi _0 - {\delta \Lambda }/{4 \Lambda \phi _0}$. At the end of this section, we argue that near the channel wall the distributions of $\alpha$ and $\beta$ get large contributions from the distributions of $\phi$ and $\theta$. Thus the stable fixed point explains the linear dependence of $\alpha _c$ on $\delta \Lambda$ near the channel wall, figure 6(b).

When $\delta \Lambda = 2 \Lambda \phi _0^2$, the two fixed points merge to one, and for $\delta \Lambda > 2 \Lambda \phi _0^2$ there are no fixed points. For $\delta \Lambda > 2 \Lambda \phi _0^2$ the linear term in $\dot {\delta \phi }$ can be ignored. For the following we restrict ourselves to the regime $\delta \Lambda \ll 2 \Lambda \phi _0^2$ and ignore the $\delta \phi ^2$ term in (4.11):

(4.12)\begin{equation} \left.\begin{array}{c@{}} \dot{\phi} = -\dfrac{s}{2} (1 -\cos 2 \phi) + \eta_\phi, \\ \dot{\delta\phi} = \dfrac{s}{2} \delta\Lambda \cos 2 \phi - (s \Lambda \sin 2\phi ) \delta\phi. \end{array}\right\} \end{equation}

Chertkov et al. (Reference Chertkov, Kolokolov, Lebedev and Turitsyn2005) and Turitsyn (Reference Turitsyn2007) argued that the $\theta$ distribution must have power-law tails in the case of a single polymer in a shear flow, with different power-law exponents arising from deterministic and stochastic $\phi$ dynamics, but were unable to analytically estimate the value of the exponent for the stochastic region. Similarly, we find that the $\delta \phi$ distribution exhibits two regimes corresponding to deterministic and stochastic $\phi$ dynamics, both of which lead to power-law tails for the distribution of $\delta \phi$ with different exponents. In addition we calculate the power-law exponent both in the deterministic and the stochastic region. The deterministic regime $|\phi | \gg \phi _0$ leads to a power-law exponent $-{3}/{2\Lambda }$ for $P(\delta \phi )$ and the stochastic regime $|\phi | \ll \phi _0$ with power law exponent $-1-{1}/{\Lambda }$ for $P(\delta \phi )$.

First consider the deterministic regime, $|\phi | \gg \phi _0, \delta \phi > {\delta \Lambda }/{4 \Lambda \phi _0}$. In this regime the terms $\eta _\phi$ for $\dot {\phi }$ and $({s}/{2})\delta \Lambda \cos 2\phi$ for $\dot {\delta \phi }$ in (4.12) can be ignored. These equations can then be integrated to obtain $\delta \phi = C \sin ^{2\Lambda }\phi$ where $C$ is an integration constant. Using $P(\phi ) \sim \phi ^{-2}$, we obtain $P(\delta \phi ) \sim \delta \phi ^{-{3}/{2\Lambda }}$ by a change of variables.

Next consider the stochastic regime $|\phi | \ll \phi _0, \delta \phi > {\delta \Lambda }/{4 \Lambda \phi _0}$. Then the term $({s}/{2})\delta \Lambda \cos 2\phi$ in (4.12) can be ignored. The steady-state Fokker–Planck equation reads,

(4.13)\begin{equation} \frac{\partial }{\partial \phi}\left[ \sin^2 \phi + \varepsilon^2 \frac{\partial}{\partial \phi}\right] P(\phi,\delta\phi) = - \Lambda \frac{\partial }{\partial \delta\phi}[\sin 2 \phi \delta\phi P(\phi,\delta\phi)], \end{equation}

where we have defined $\varepsilon ^2 = {D}/{2s}$ in analogy with Meibohm et al. (Reference Meibohm, Pistone, Gustavsson and Mehlig2017). Numerics suggest that the joint probability factorises in the stochastic regime, so that $P(\phi , \delta \phi ) = f(\phi ) g(\delta \phi )$. We use separation of variables and obtain $g(\delta \phi ) \propto \delta \phi ^{-1-{\mu }/{\Lambda }}$. The equation for $f(\phi )$ is a generalised eigenvalue problem,

(4.14)\begin{equation} \frac{\partial }{\partial \phi}\left[ \sin^2 \phi + \varepsilon^2 \frac{\partial}{\partial \phi}\right] f(\phi) = \mu \sin 2\phi \ f(\phi). \end{equation}

It is possible to obtain an eigenvector corresponding to the eigenvalue $\mu = 1$ for (4.14). The eigenvector has two undetermined integration constants which must be found by matching to the solution for large $\phi$ and normalisation. Thus we conclude that $g(\delta \phi ) = \delta \phi ^{-1-{1}/{\Lambda }}$ is a solution for the tail of the distribution of $\delta \phi$ in the regime $|\phi | \ll \phi _0, \delta \phi > {\delta \Lambda }/{4 \Lambda \phi _0}$. When $\Lambda \approx 1$ this would lead to a power law for $\delta \phi$ with exponent $\approx -2$. In figure 4(a) we have plotted the distribution for $\alpha$ which shows a power-law distribution for large $\alpha$, with exponent roughly $-2$. We argue next that since $\alpha$ must be closely related to $\delta \phi$ near the channel wall, the observed exponent $-2$ for $\alpha$ can be explained by our calculation for $\delta \phi$.

So far we have analysed the dynamics for the angle $\phi$ in the channel frame for particles with shape parameter $\Lambda$. Since the Lagrangian stretching direction spends long times aligned with $\hat {\boldsymbol {x}}$, and the Lagrangian contracting direction spends long times aligned with $\hat {\boldsymbol {z}}$, $\alpha , \beta$ get large contributions from $\delta \phi , \delta \theta$, respectively. This means that near the channel wall one can use $\delta \phi$ and $\delta \theta$ as a proxy for understanding $\alpha$ and $\beta$, respectively. The precise calculation of the distributions of $\alpha , \beta$ is an open question left for future work.

In the presence of strong mean shear and weak velocity-gradient fluctuations we have shown that the width of the plateau in the distribution of relative angles scales linearly with $\delta \Lambda$ and that the power-law exponent for the tail of the distribution is $-{3}/{2\Lambda }$ and $-1-{1}/{\Lambda }$ in the deterministic and stochastic regimes, respectively. Consider the relative angle in the case of constant shear without fluctuations. In the long time, the Lagrangian stretching direction reaches a steady state $\phi =0$, and thus the distribution of $\delta \phi$ has a plateau whose width scales as $\delta \Lambda ^{{1}/{2}}$ and an exponent $-2$ for the power-law tail of the distribution. Thus, the presence of weak fluctuations affects the relative angular dynamics sensitively and the observations cannot be explained in terms of just the strong mean shear.