Effective aspect ratio of helices in shear flow

TL;DR

This study investigates how rigid colloidal helices behave in shear flow, revealing that their effective aspect ratio influences their rotation and alignment, with implications for understanding their dynamics in fluid environments.

Contribution

The paper introduces geometric and dynamical effective aspect ratios for helices in shear flow and demonstrates their approximate equivalence across various conditions.

Findings

Elongated helices align more with the flow direction.

Transient axis deflections depend on helix rotation.

Effective aspect ratios are consistent across parameters.

Abstract

We report the results of simulations of rigid colloidal helices suspended in a shear flow, using dissipative particle dynamics for a coarse-grained representation of the suspending fluid, as well as deterministic trajectories of non-Brownian helices calculated from the resistance tensor derived under the slender-body approximation. The shear flow produces nonuniform rotation of the helices, similarly to other high aspect ratio particles, such that more elongated helices spend more time aligned with the fluid velocity. We introduce a geometric effective aspect ratio calculated directly from the helix geometry and a dynamical effective aspect ratio derived from the trajectories of the particles and find that the two effective aspect ratios are approximately equal over the entire parameter range tested. We also describe observed transient deflections of the helical axis into the vorticity direction that can occur when the helix is rotating through the gradient direction and that depend on the rotation of the helix about its axis.

Equations

$${\mathbf{F}}_i=\sum _{j\ne i}\left(F_{ij}^C+F_{ij}^R+F_{ij}^D\right){\hat{r}}_{ij},$$

$$F_{ij}^C=Aw(r_{ij}),$$

$$F_{ij}^D=-\gamma w^2(r_{ij})({\hat{r}}_{ij}\cdot v_{ij}),$$

$$F_{ij}^R=\sigma w(r_{ij}){\alpha }_{ij}/\sqrt{\Delta t},$$

$$\displaystyle w(r)=\left\{\begin{array}[]{lr}1-r/r_{c}&r\leq r_{c}\\ 0&r>r_{c}.\end{array}\right.$$

$$\begin{pmatrix}x(u)\\ y(u)\\ z(u)\end{pmatrix}=\begin{pmatrix}r\cos(2\pi nu/\ell_{c})\\ \ell u/\ell_{c}\\ hr\sin(2\pi nu/\ell_{c})\end{pmatrix}\!\!\!\mbox{\quad with\quad}\!\!\!\begin{array}[]{l}0\leq u\leq\ell_{c}\\ h=\pm 1,\end{array}$$

$$\mathbf{F}={\int }_0^{{\ell }_c}\text{differential}u\left({\alpha }_t{v}_r^t+{\alpha }_n{v}_r^n\right),$$

$$\tau ={\int }_0^{{\ell }_c}\text{differential}ur\text{crossproduct}\left({\alpha }_t{v}_r^t+{\alpha }_n{v}_r^n\right),$$

$$\int \text{differential}u\left({\alpha }_t{v}_f^t+{\alpha }_n{v}_f^n\right)=\int \text{differential}u\left({\alpha }_t{v}_h^t+{\alpha }_n{v}_h^n\right),$$

$$\int \text{differential}ur\text{crossproduct}\left({\alpha }_t{v}_f^t+{\alpha }_n{v}_f^n\right)=\int \text{differential}ur\text{crossproduct}\left({\alpha }_t{v}_h^t+{\alpha }_n{v}_h^n\right).$$

$$\left(\begin{array}{c}{\mathbf{F}}_f\\ {\tau }_f\end{array}\right)=\overline{\overline{\mathbf{R}}}\cdot \left(\begin{array}{c}{v}_{cm}\\ \omega \end{array}\right),$$

$$f_d={\alpha }_t{v}_f^t+{\alpha }_n{v}_f^n={\alpha }_n{v}_f+\left({\alpha }_t-{\alpha }_n\right)v_f^t,$$

$${\mathbf{F}}_f={\int }_0^{{\ell }_c}\text{differential}u{f}_d\text{mbox}and{\tau }_f={\int }_0^{{\ell }_c}\text{differential}u\mathbf{H}\text{crossproduct}{f}_d.$$

$$r_{\psi }=\sqrt{\frac{{\tau }_1}{{\tau }_2}}=\sqrt{\frac{2}{3}}\frac{\ell }{r}\sqrt{\frac{3r^2({\pi }^2{n}^2+\cos (2\psi ))+{\ell }^2}{(4\pi nr{)}^2+{\ell }^2(2+\cos (2\psi ))}},$$

$$r_G=r_{\psi =\pi /4}=\frac{a}{\sqrt{3}}\sqrt{\frac{a^2+3{\pi }^2{n}^2}{a^2+8{\pi }^2{n}^2}}.$$

$$r_e=\frac{1}{⟨{\cos }^2\varphi ⟩}-1,$$

$${\tau }_{\theta }^0=\frac{\dot{\gamma }{\alpha }_{t}r\ell \cos (\psi )\left(h\sin (2\varphi )\left(9{\ell }^2-6{\ell }_c^2-4{\pi }^2{n}^2{r}^2{\cos }^2(\psi )\right)+3\pi nr\ell \sin (\psi )(\cos (2\varphi )+3)\right)}{12\pi n{\ell }_c}.$$

Full text

I Introduction

The behavior of flowing suspensions of fibers and other high aspect ratio particles plays an important role in a wide range of commercially important processes and consequently has been intensively investigated Lundell et al. (2011); Wang et al. (2014). Fibers with intrinsic curvature, including those with helical shapes, appear in many contexts Wang et al. (2014); Zhang et al. (1998); Kagawa et al. (1985), and their behavior under flow represents an interesting and challenging fundamental problem Leal (1980).

A number of studies have specifically investigated the behavior of rigid helices in shear flow experimentally and computationally Kim and Rae (1991); Makino and Doi (2005); Kostur et al. (2006); Marcos et al. (2009); Chen and Zhang (2011); Zhang (2017). Like other high aspect ratio particles such as rods or ellipsoids, shear flow will cause helical filaments to rotate about an axis perpendicular to the velocity gradient and the flow directions (the vorticity axis) with a non-uniform rotation rate, with the particles spending more time with their long axis parallel to the flow than parallel to the shear gradient. For rotationally symmetric ellipsoids at low Reynolds number and in the absence of Brownian motion, the trajectories represent closed orbits (Jeffery Orbits) with analytic solutions that depend only on the particle’s aspect ratio and its orientation with respect to the gradient axis, with no net motion in the vorticity direction Dhont and Briels (2007). More generally, most axisymmetric bodies with fore-aft symmetry will also follow Jeffery Orbits in shear flow, with an effective aspect ratio that is determined by the square root of the ratio of the torque exerted on the body when it is held at rest with its axis along the gradient direction to the torque with its axis along the flow direction Cox (1971), although there are interesting exceptions in special cases Singh et al. (2013).

Helical particles, by contrast, do not follow closed orbits, even in the absence of thermal fluctuations Kim and Rae (1991); Kostur et al. (2006). Furthermore, the particles experience a net drift in the vorticity direction with a sign dependent on the helicity, a phenomenon which has been exploited to use shear flows to separate chiral objects Kim and Rae (1991); Kostur et al. (2006); Marcos et al. (2009); Chen and Zhang (2011); Kramel et al. (2016); Ro et al. (2016). While considerable progress has been made understanding the average long time behavior of helices in shear flow in the presence of thermal fluctuations, we currently lack the ability to predict the short term dynamics of helical filaments, information that is critical for understanding the role of helical particles in suspension rheology or the behavior of particles in complex flows such as turbulence Kramel et al. (2016).

In this work, we report the results of Dissipative Particle Dynamics (DPD) computer simulations of rigid helices in the presence of a shear flow. The DPD technique produces stochastic forces similar to thermal fluctuations and we observe that the orbits of the helices are qualitatively similar to noisy Jeffery Orbits. We compare the simulated trajectories with the deterministic trajectories calculated from the equation of motion for non-Brownian helices in the slender body approximation and report analytic expressions for the forces, torques, resistance tensor, center of mass velocity and angular velocity of a general helix in arbitrary orientation. We derive an analytic expression for a geometric effective aspect ratio calculated directly from the helix geometry and compare that to a dynamical effective aspect ratio calculated from the trajectories of the helices in the simulations and the deterministic trajectories. Over the entire parameter range tested, the geometric aspect ratio matches the measured dynamical aspect ratio in the simulations, within the statistical uncertainty, and accurately predicts the dynamical aspect ratio calculated from the deterministic trajectories. Finally, we discuss the origins of transient deflections of the helical axis into the vorticity direction that occur while the helix is rotating through the gradient direction in some, but not all, of the trajectories.

II Methods

II.1 Simulations

There have been many simulations of fibers in fluid flow using various approximations of hydrodynamic and contact interactions, employing a variety of techniques available with varying degrees of complexity and accuracy Bolintineanu et al. (2014). For these studies, we have used Dissipative Particle Dynamics (DPD) (Hoogerbrugge and Koelman, 1992; Groot and Warren, 1997), an efficient coarse-grained fluid representation that can capture many aspects of the complex hydrodynamic interactions between the helical filament and the surrounding fluid and the effects of thermal fluctuations (Fan, 2006; Duong-Hong, 2006) and is relatively simple to implement. The DPD implementation is similar to one we have used previously to study shear induced aggregation of straight rods Stimatze et al. (2016) and is briefly summarized below.

In DPD, the coarse-grained fluid is represented by soft particles interacting via three pairwise forces: a repulsive force that determines the compressibility of the fluid, a dissipative force that models viscous dissipation, and a random force that determines the steady state temperature of the system.

Thus the total force on particle $i$ is

[TABLE]

where

[TABLE]

is the conservative, soft repulsion contribution to the force exerted by particle $j$ on particle $i$, where ${{\bf\it r}}_{ij}={{\bf\it r}}_{i}-{{\bf\it r}}_{j}$ and $r_{ij}=|{{\bf\it r}}_{ij}|$. The dissipative force is

[TABLE]

where ${\bf\it v}_{ij}=\differential{{\bf\it r}}_{ij}/\differential t$, and the random force is

[TABLE]

where $\alpha_{ij}$ is a random variable with unit variance Gaussian statistics, and $w(r)$ is a weighting function given by

[TABLE]

Following Groot & Warren Groot and Warren (1997), we set the repulsive parameter $A=18.75$, the density $\rho=4$, the random force coefficient $\sigma=3$, the force cutoff radius $r_{c}=1$, and the dissipation coefficient $\gamma=4.5$. The combination of $\rho$ and $A$ determines the compressibility, and values chosen are consistent with the compressibility of water. The combination of the strength of the random force $\sigma$ and the dissipation coefficient $\gamma$ determines the steady state kinetic energy (effective temperature) of the DPD system. For the values used here, $T=1$ in simulation units, which determined the appropriate timestep for the calculations (we use $\Delta t=0.01$)Groot and Warren (1997). All numerical values are given in simulation units, with the relevant length scale being the particle size (unit diameter) and the time scale set by the applied shear.

With parameters in this range, DPD has been shown to reproduce correct hydrodynamics at long length scales Espanol and Warren (1995). The helix length scales simulated in this work are not large compared to the DPD particle size, however, so quantitative agreement with Navier-Stokes hydrodynamics is not expected.

Shear flow is generated by directly simulating moving boundaries at the top and bottom of the simulated fluid.

Specifically, referring to Fig.1, fixed walls one particle diameter thick in the horizontal $(x,z)$ plane are moved at constant, opposite speeds ($v_{\text{wall}}=1.5$) in the $x$ direction, thus producing simple shear with $(x,y,-z)$ = (flow, gradient, vorticity) directions. Periodic boundaries are employed in the $x$ and $z$ directions. For the simulations below, the simulated domain is cubic with a size that is adjusted according to the parameters of the helix being simulated to minimize self-interactions across periodic boundaries as well as contact with walls. The side length of the box was normally set as 1.5 times the helix length, rounded up to the nearest increment of 5 units. No contact with walls was observed for any simulations reported here.

Helices are simulated by rigid strands of spherical particles with centers separated by a fixed spacing of half-unit length in simulation units. The particles interact with the fluid particles by the same DPD interactions described above. The initial helix configuration was set to be

[TABLE]

which gives a left(right)-handed helix when $h=1$($h=-1$) of filament length $\ell_{c}=\sqrt{\ell^{2}+4\pi^{2}n^{2}r^{2}}$ with $n$ turns, radius $r$, end-to-end length $l$ and pitch $p=l/r$, initially oriented parallel to the gradient direction such that a perpendicular from the helical axis to the first bead on the $+y$ end of the helix was in the positive flow direction ($(\phi,\theta,\psi)=(0,0,0)$, see Fig. 2). All parameter sets were run for 50,000 timesteps, then extended until at least 2 “flips” were observed, representing at least one full orbit.

One hundred and five simulations of isolated helices with varying pitch ($p$), radius ($r$) and length ($\ell$) were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Plimpton (1995) environment, including the LAMMPS implementation of DPD and the LAMMPS rigid-body integrator for calculating net forces and torques on the helices. The resulting dynamical equations were solved using the velocity-Verlet integration scheme.

II.2 Analytic Methods

We also present results from deterministic trajectories calculated with the slender body approximation, appropriate for thin filaments in Stokes (low Reynolds number) flow. The equations of motion are determined from the resistance tensor for the helix under the conditions of zero net force and torque.

The force on the helix from slender body theory is given by

[TABLE]

and the torque by

[TABLE]

where ${\bf\it v}_{r}$ is the velocity of the fluid relative to the helical fiber, with superscript $t\ (n)$ representing the component tangent (normal) to the fiber, $\alpha_{t}\ (\alpha_{n})$ is the drag coefficient for flow tangent (normal) to the fiber, and ${\bf\it r}$ is the vector from the helix center to a point on the fiber. Being in a Stokes flow regime immediately implies that $\mathbf{F}={\bf\it\tau}=\mathbf{0}$.

We can break ${\bf\it v}_{r}$ into contributions from the flow and from the helix’s motion as ${\bf\it v}_{r}={\bf\it v}_{f}-{\bf\it v}_{cm}-{\bf\it\omega}\crossproduct{\bf\it r}$ where ${\bf\it v}_{f}$ is the fluid velocity, ${\bf\it v}_{cm}$ is the helix center of mass velocity, and ${\bf\it\omega}$ is the angular velocity of the helix. This allows us to recast Eq. 3 and 4, defining ${\bf\it v}_{h}\equiv{\bf\it v}_{cm}+{\bf\it\omega}\crossproduct{\bf\it r}$ and dropping the integration bounds, as

[TABLE]

and

[TABLE]

Together these two equations may be written as a single matrix equation

[TABLE]

where $\mathbf{F}_{f}$ and ${\bf\it\tau}_{f}$ are the force and torque on a motionless (i.e. anchored in place) helix due to the unperturbed fluid flow, and $\overline{\overline{\mathbf{R}}}$ is the $6\times 6$ symmetric resistance tensor for the helix. We are able to analytically find $\mathbf{F}_{f}$, ${\bf\it\tau}_{f}$, and $\overline{\overline{\mathbf{R}}}$, then solve Eq. 7 to find ${\bf\it v}_{cm}$ and ${\bf\it\omega}$ for an arbitrary helix in an arbitrary orientation.

A helix in an arbitrary orientation can be described by the vector $\mathbf{H}$, obtained by applying rotation matrices for the three angles $(\phi,\theta,\psi)$ to the vector described by Eq. 2.

The drag on the segments of the filament that compose the helix is anisotropic. From slender-body theory at low Reynolds number, the drag coefficient for flow normal to the filament is twice the drag coefficient for flow tangent to the filament, i.e. $\alpha_{n}/\alpha_{t}=2$ Batchelor (1970). Allowing for a nonzero filament thickness results in $\alpha_{n}/\alpha_{t}=2$ to within a logarithmic correction term Marcos and Stocker (2010); Childress (1981). We take $\alpha_{n}=2\alpha_{t}$ following slender body theory throughout the rest of the paper unless explicitly stated otherwise.

We decompose the unperturbed fluid flow, given by ${\bf\it v}_{f}=\dot{\gamma}(\mathbf{H}\cdot\hat{y})\hat{x}$, into tangent, ${\bf\it v}_{f}^{t}=({\bf\it v}_{f}\cdot\hat{t})\hat{t}$, and normal, ${\bf\it v}_{f}^{n}={\bf\it v}_{f}-{\bf\it v}_{f}^{t}$, components where $\hat{t}$ is a unit vector parallel to the filament (i.e. in the direction of $\derivative{\mathbf{H}}{u}$). The drag force per unit length is given by

[TABLE]

and so the force and torque on the fixed helix are

[TABLE]

The resulting expressions are very messy, but can be easily evaluated for any helix parameters and orientation, and are provided in the Supplemental Material sup .

Using the same method as described above, we decompose the helix velocity, ${\bf\it v}={\bf\it v}_{cm}+\mathbf{H}\crossproduct{\bf\it\omega}$, into normal and tangent components. This allows us to evaluate the RHS of Eqs. 5 and 6 from which $\overline{\overline{\mathbf{R}}}$ may be directly read off from the coefficients of ${\bf\it v}_{cm}$ and ${\bf\it\omega}$. The full expression for $\overline{\overline{\mathbf{R}}}$ is provided in the Supplemental Material sup . To our knowledge, this is the first explicit expression for the resistance tensor of an arbitrary helix in an arbitrary orientation.

We then solve Eq. 7 for ${\bf\it v}_{cm}$ and ${\bf\it\omega}$. The resulting expressions are too cumbersome to report but are available in the Supplemental Material as a Mathematica binary sup . Having analytic expressions for the center of mass and angular velocities enabled us to numerically integrate the coupled differential equations governing the dynamics of the helix at a significantly lower computational cost than was previously possible. We are thus able to compare the noisy DPD simulations with the dynamics determined from numerically integrating the solution to the deterministic equations of motion for Stokes flow in the slender body approximation.

III Geometric effective aspect ratio

As described in the introduction, the trajectories of an axisymmetric body with fore-aft symmetry in viscous shear flow are given by Jeffery Orbits that are fully determined by the ratio of the torque exerted on the body when it is held at rest with its axis along the gradient direction to the torque with its axis along the flow direction Cox (1971). Inspired by this result, we use this quantity for helices to define a geometric effective aspect ratio for helices, despite the fact that they are not axisymmetric, nor do they have fore-aft symmetry.

Following Cox Cox (1971), we define a ($\psi$-dependent) geometrical aspect ratio for the helix as

[TABLE]

where $\tau_{1}$ is ${\bf\it\tau}_{f}\cdot\hat{z}$ (as defined in Eq. 4) when the helical axis is parallel to the gradient direction $\hat{y}$ ($\phi=0,\ \theta=0$) and $\tau_{2}$ is ${\bf\it\tau}_{f}\cdot\hat{z}$ when the helical axis is parallel to the flow direction $\hat{x}$ ($\phi=\pi/2,\ \theta=0$).

Unlike the axisymmetric bodies studied by Cox, this definition does not uniquely specify the aspect ratio for a given geometry due to the dependence on the angle $\psi$. Moreover, given that $\psi$ will change during an orbit, the situation is clearly more complicated, and the simple ratio of torques at fixed position will not be sufficient to precisely determine the trajectory. Nonetheless, we can still use Eq. 10 as an aspect ratio that depends only on the geometry and $\psi$ and investigate the extent to which that quantity is a useful predictor of the actual trajectories. Figure 3 shows how $r_{\psi}$ (Eq. 10) varies with $a=\ell/r$, $n$ and $\psi$.

Our goal is to find an aspect ratio based solely on the geometry of the helix that is a good predictor of the fraction of time helices spend aligned with the shear flow for a wide range of helical parameters. We note that for large $n$ and $a$, $r_{\psi}\approx{a}/({2\sqrt{2}})\approx 0.35a$, with significant deviations only for small $n$ (large pitch). Indeed $\lim_{n\to\infty}r_{\psi}=a/(2\sqrt{2})$ for all $a$ and $\psi$. This suggest an approximate geometrical aspect ratio based solely on $a$, $r_{a}=a/(2\sqrt{2})$. As can be seen in Fig. 3, $r_{\psi}\approx r_{a}$ for $n\geq 5$. However, as discussed below, we find that $r_{a}$ is not a good predictor of the fraction of time helices spend aligned with the shear flow for much of the parameter range consider here.

It is perhaps not surprising that helix parameters beyond the bounding cylinder aspect ratio need to be taken into account to accurately describe the tumbling trajectories. Thus we introduce an alternative purely geometrical aspect ratio, which we denote as $r_{G}$, produced by taking $\psi=\pi/4$ in Eq. 10. $r_{G}$ can be calculated from Eq. 10 and is given by

[TABLE]

Below we show that $r_{G}$ works quite well as a predictor of the fraction of time helices spend aligned with the shear flow for all of the parameters studied here, despite the variation of $\psi$ during a trajectory.

IV Results

IV.1 Jeffery-like Orbits

Trajectories for isolated helices were initialized with their helical axis oriented in the gradient direction. In this configuration the shear flow exerts a torque on the helix parallel to the vorticity axis, resulting in a rapid rotation into the flow direction. The torque is reduced as the helix aligns with the flow, so the rotation rate decreases, reaching a minimum when the axis of the helix is perpendicular to the gradient ($y$) axis. This behavior is qualitatively similar to Jeffery Orbits of non-Brownian axisymmetric ellipsoids, which are deterministic, closed orbits with a fixed angle, $\theta$, relative to the vorticity direction. As an example, Fig. 4 displays an orbit for an ellipsoid of revolution with aspect ratio $r_{e}=a/2=4$ in a series of orientation snapshots and plots of $\phi$ and $\theta$ vs. time.

An alternative way to visualize the trajectories is to track the evolution of components of the orientation vector, $\hat{u}$, a unit vector aligned with the axis of the helix. We define $\tilde{u}$, the normalized projection of $\hat{u}$ onto the flow-gradient ($xy$) plane, to more easily visualize the rotation of the helix about the vorticity ($-z$) axis. It is then easy to characterize the orbiting behavior of the helix using $\tilde{u}_{y}\equiv\tilde{u}\cdot\hat{y}=\hat{u}\cdot\hat{y}/\cos{\theta}=\cos{\phi}$ and the deflections into the vorticity direction by $u_{z}\equiv\hat{u}\cdot\hat{z}=\sin{\theta}$. The dotted blue curve in Fig. 4 shows $\tilde{u}_{y}^{2}$, vs. strain ($\dot{\gamma}t)$, where the relatively slow rotation rate when the particle is aligned in the flow direction produces an extended period of time when $u_{y}^{2}$ is small.

Figure 5 displays results from two representative helical geometries, each showing a deterministic calculation (open circles) and a DPD simulation (filled), revealing both the non-uniform rotation rates and the stochastic variations that arise in the DPD simulations. As expected, we observe that squat helices, with smaller length to radius ($a=\ell/r$) values, display relatively small variation in their rotation rates (top panel, $a=5$), whereas helices with high values of $a$ show rotation rates that slow down dramatically when aligned in the flow direction (bottom panel, $a=15$).

For the parameter regimes investigated here, the stochastic component of the motion in the simulations is relatively small, suggesting that the effective rotational Peclet number, the ratio of the shear rate to the rotation diffusion coefficient, $D_{r}$, is large. We reported previously that for a rigid rod of 21 particles under identical conditions, $D_{r}=4\times 10^{-6}$ in simulation units Stimatze et al. (2016), and in general for rods, $D_{r}\propto L^{-3}$.

IV.2 Effective Aspect Ratio

For ideal Jeffery Orbits, the angle $\phi$ of the particle relative to velocity direction is related to the aspect ratio $r_{e}$ according to

[TABLE]

where brackets represent the long time average. Experimentally or computationally accessing the long time average is clearly difficult, but we observe that from Jeffery’s equations (and as visualized in Fig. 4), $\cos^{2}\phi$ is periodic with a period of half an orbit. Furthermore, each half orbit is mirror symmetric about its midpoint. This means that the long time average, for particles obeying Jeffery’s equations, of $\cos^{2}\phi$ is equal to the average over an integral number of quarter orbits.

A spherical particle $(r_{e}=1)$ has a uniform rotation rate, producing $\langle\cos^{2}\phi\rangle=1/2$. As $r_{e}$ increases, the particle spends a longer fraction of its orbit aligned in the flow direction, so $\langle\cos^{2}\phi\rangle$ (and $\langle\tilde{u}_{y}^{2}\rangle$) decreases. Note that this result is independent of $\theta$, the angle with respect to the gradient-velocity plane (which is constant for a Jeffery Orbit and therefore determined uniquely by the initial conditions).

Equation 12 can be easily generalized to calculate an effective aspect ratio from our simulated trajectories (as in Stimatze et al. (2016)), with a couple of caveats. One is that $\theta$ is not constant, varying due to both thermal fluctuations and torques with components in the flow-gradient plane. When $\theta$ approaches $\pi/2$ (helix aligned in the vorticity direction), thermal noise causes $\phi$ to fluctuate erratically. This issue does not create difficulties in this work, where we focus on relatively short trajectories with initial conditions of $\theta=0$. A second caveat is that, as discussed above, the effective aspect ratio calculation from measured trajectories requires an integral number of quarter orbits, which can create selection bias for finite length trajectories. In order to minimize this effect, we have included only parameter ranges where all simulated trajectories included at least one full rotation. We further discard the first half rotation to mitigate any possible transient effects in the simulation.

Here we seek to determine if the empirical effective aspect ratio $r_{e}$, calculated from Eq. 12, can be simply related to geometric parameters of the helix. Figure 6A shows a scatter plot of $r_{e}$, determined from the simulations, versus $r_{G}$, defined by Eq. 11, for a range of helix parameters. Although there is considerable scatter in the data, there is a strong correlation between the two quantities, with no evident dependence on either $\ell$ or $p$ (indicated by the size and color of the points, respectively). A least squares fit to $\log{r_{e}}$ vs. $\log{r_{G}}$ produces $r_{e}=1.01r_{G}^{1.02}$ and is consistent with the two quantities being equivalent ($r_{e}=Ar_{G}^{\alpha}$, with $\alpha=\{0.93,1.10\}$ and $A=\{0.89,1.15\}$ at a 95% confidence level (Cl).

While the geometric aspect ratio defined by Eq. 11 does a remarkably good job of predicting $r_{e}$, it is worth noting that, as discussed above, over much of the parameter range, $r_{G}$ is close to $r_{a}={a}/({2\sqrt{2}})$ (see Fig. 3). A plot of $r_{e}$ vs $r_{a}$ looks qualitatively similar to Fig. 6A, but shows systematic deviations for small $n$. Our results are inconsistent with $r_{e}=r_{a}$ ($r_{e}=A_{1}r_{a}^{\alpha_{1}}$, with $\alpha_{1}=\{1.04,1.22\}$ and $A_{1}=\{0.84,1.08\}$ at a 95% Cl.

We find that $r_{e}$ and $r_{G}$ are similarly closely related for the deterministic helix trajectories, as shown in Fig. 6B, for the same helix parameters used in Fig. 6A, but with initial values of $\psi$ distributed between 0 and $2\pi$. A least squares fit to $\log{r_{e}}$ vs. $\log{r_{G}}$ here produces $r_{e}=1.01r_{G}^{1.01}$. The scatter in the data clearly shows that $r_{e}$ depends on $\psi$, suggesting that the scatter in Fig. 6A arises from a combination of the effects of Brownian motion and the limitations of the definition of a $\psi$-independent $r_{G}$. Fig. 6C shows the result of plotting the same data for $r_{e}$ versus $r_{\psi}$, the $\psi$-dependent geometric aspect ratio calculated from Eq. 10, using the initial value of $\psi$ for each trajectory. This appears to reduce the scatter and gives a least squares fit to $\log{r_{e}}$ vs. $\log{r_{\psi}}$ as $r_{e}=1.01r_{\psi}^{1.00}$. Note that $\psi$ can change significantly over the course of the orbit, as shown below.

These results show that $r_{G}$ provides a reasonably robust estimate of $r_{e}$ for both the DPD simulations and the deterministic trajectories calculated using slender body theory, but that the absence of axisymmetry manifests itself in a dependence on the angle $\psi$ that cannot be easily captured. To our knowledge, this is the first test of an effective aspect ratio of helices in shear flow. Marcos et al. Marcos et al. (2009) numerically calculate effective aspect ratios for helical bacteria based on the ratio of rotation rates for $\phi=0$ and $\pi/2$, which produces very similar results to the $r_{G}$ defined here, but is a significantly messier expression sup .

IV.3 Deflections in Vorticity Direction

Although the primary focus of this study is the rotation about the vorticity axis, we close by noting an interesting behavior that occurs in many, but not all, of the simulated and deterministic trajectories. Figure 7 shows the evolution of $\tilde{u}_{y}^{2}$ and $u_{z}^{2}$ for a representative trajectory from the simulation, where we find a transient deflection in the vorticity direction that peaks when $\tilde{u}_{y}^{2}$ is maximum, i.e when $\phi=n\pi$ and the rotation rate of the Jeffery-like Orbit is at its maximum. Similar, albeit smaller, deflections are observed in the deterministic trajectory for the same initial condition. For both simulated and deterministic trajectories, the deflections vary in magnitude and sign, with no obvious pattern, and are sometimes absent altogether.

We can gain some insight into the origin of these transient deflections by considering the torque on a fixed helix (LHS of Eq. 6). The torque causing deflections into the vorticity direction is given by $\tau_{\theta}={\bf\it\tau}\cdot\hat{\theta}\equiv{\bf\it\tau}\cdot\left(\hat{x}\cos(\phi)+\hat{y}\sin(\phi)\right)$. This quantity ($\tau_{\theta}$) is plotted in the lower panels of Fig. 7, using the instantaneous values of $\phi$, $\theta$ and $\psi$, alongside plots of $\theta(t)$ and $\psi(t)$. It is clear that $\tau_{\theta}$ is highly correlated with $\dot{\theta}$, although it should be noted that the exact behavior of the vorticity deflections is characterized by $\omega_{\theta}\equiv{\bf\it\omega}\cdot\hat{\theta}$ which is an extremely complicated expression.

The full expression for $\tau_{\theta}$ itself is quite unruly, but we can capture some aspects of the relevant behavior if we restrict our attention to a helix in the gradient-velocity plane ($\tau_{\theta}$ at $\theta=0$). This quantity, which we label $\tau_{\theta}^{0}$, is found to be

[TABLE]

While the full expression for $\tau_{\theta}^{0}$ is still complex, in the limit of small $r$ it reduces to simply $\tau_{\theta}^{0}\approx\dot{\gamma}\alpha_{t}h\ell^{2}r\cos(\psi)\sin(2\phi)/(4\pi n)$. The $\sin(2\phi)$ dependence produces a contribution that has extrema at $\phi=\pi/4$ and $3\pi/4$, with a sign change in between, and thus explains the transient deflection. The origin of this dependence can be understood by recalling that simple shear can be represented as a superposition of pure shear (or extensional flow) and pure rotational flow. The contribution from the rotational flow is $\phi$-independent (the very last term in Eq. 13, of order $r^{2}$). The extensional flow is outward along the extensional axis ($\phi=\pi/4$) and inward along the compressional axis ($\phi=3\pi/4$), and thus is responsible for the $\sin(2\phi)$ term in $\tau_{\theta}^{0}$. At $\phi=\pi/4$ and $3\pi/4$, the axial flow produces drag forces that are primarily parallel to the helical axis (although tilted towards the filament normal because of the anisotropic drag), but will contribute to a $\psi$-dependent $\tau_{\theta}^{0}$ because the finite radius of the filament means that each filament element generally has an $r_{z}$ that is non-zero, so $\vec{r}\crossproduct\vec{F}$ can contain elements that contribute to $\tau_{\theta}^{0}$. Interestingly, we find that the expression for $\tau_{\theta}^{0}$ in the limit of small $r$ given above is also valid for isotropic drag ($\alpha_{t}=\alpha_{n}$), indicating that the general behavior can be understood from this geometric argument, without consideration of the tilting of the drag force relative to the flow direction.

Fig. 8, which shows a surface plot of $\tau_{\theta}^{0}(\phi,\psi)$ for a representative set of helix parameters, has a dominant $\sin(2\phi)$ dependence that is modulated by a sinusoidal function of $\psi$. This likely underlies the variability in the sign and magnitude of the observed vorticity deflections. Furthermore, $\psi$ evolves continuously during the orbit, as can be seen in Fig. 7, and this evolution presumably leads to the complexity of the helical trajectories even in the absence of Brownian motion.

V Discussion

As described in the Introduction, rigid helical filaments initially oriented parallel to a shear gradient will rotate about the vorticity axis with a rotation rate that decreases as the helix aligns in the flow direction. Our simulations and deterministic trajectories show, as expected, that the reduction of the rotation rate increases with the aspect ratio of the bounding cylinder ($a/2=l/2r$, Fig. 2). We find that $a$ is the primary determinant of the degree of alignment and therefore that it is relatively insensitive to the other dimensionless ratios that describe a particular helix (such as the number of turns, $n=\ell/p$, and the pitch angle (related to $r/p$)), at least in the range simulated here ($4<a<20$, $1<n<29$, $0.07<r/p<2$). (The thickness of helical filament itself, one particle diameter, was not varied.)

Qualitatively, this behavior can be understood as follows: when the helix is aligned in the gradient direction, the effect of the fluid drag will cause rotation about the vorticity axis with a rate that is comparable to the shear rate. When the helix is aligned in the velocity direction, the torque due to the shear flow, and therefore the rotation rate, is reduced, as with other high aspect ratio particles. The ratio of the torques in those two orientations is mostly determined by $a=\ell/r$. Decreasing the pitch $p$ (or, equivalently, increasing $n$) for a given $\ell$ will increase the torque overall, but not the asymmetry.

We do find, however, that there are significant deviations from the degree of alignment that would be predicted by only considering $\ell/r$, particularly for small values of $n$. The geometric aspect ratio defined by Eq. 11, based on the ratio of the torque exerted on the helix held at rest with its axis along the shear to the torque with its axis in the flow direction using slender body theory appropriate for low Reynolds number flow, accurately accounts for those deviations. In fact the data displayed in Fig. 6A shows that $r_{G}$ is approximately equal to $r_{e}$, to within statistical uncertainty, over the entire parameter range simulated. However, there is considerable scatter in the data for large aspect ratio helices because of the relatively small number of complete orbits observed, leaving the possibility that the behavior in some regimes is more complex. Using deterministic trajectories calculated in the slender body approximation, we are able to show that there is an appreciable dependence of $r_{e}$ on the angle $\psi$ that is not accounted for in Eq. 11.

While the extent of the flow alignment of the helices appears to be relatively insensitive to the pitch, it would likely impact other important physical quantities. As the pitch gets very small, the tightly wound helix will approach a rigid cylindrical shell, which we expect would rotate with nearly complete fluid entrainment except at the ends. By contrast, if $p$ is large (compared to $r$), the amount of fluid displaced by the helix will be determined by filament length, with a logarithmic dependence on the filament thickness, and will be relatively insensitive to the helix radius. Finally, in this study we have only considered isolated helices, but interactions between helices, such as the nature of entanglements, will likely depend on $p/r$.

Acknowledgements.

This work was supported by the AFOSR (FA9550-10-1-0473 and FA9550-14-1-0171). JSU is supported in part by the Georgetown Interdisciplinary Chair in Science Fund. We thank Peter Olmsted for helpful discussions.